[cig-commits] r15673 - in doc/geodynamics.org/benchmarks/trunk: . cs geodyn long magma mc mc/2d-cartesian mc/3d-spherical seismo short short/benchmark-landers short/benchmark-rs short/benchmark-rs/results short/benchmark-rs-nog short/benchmark-rs-nog/geofest-input short/benchmark-rs-nog/plots short/benchmark-rs-nog/pylith-0.8-input short/benchmark-rs-nog/results short/benchmark-strikeslip short/benchmark-strikeslip/geofest-input short/benchmark-strikeslip/plots short/benchmark-strikeslip/pylith-0.8-input short/benchmark-strikeslip/results short/savage-prescott short/utilities
luis at geodynamics.org
luis at geodynamics.org
Mon Sep 14 16:45:38 PDT 2009
Author: luis
Date: 2009-09-14 16:45:35 -0700 (Mon, 14 Sep 2009)
New Revision: 15673
Added:
doc/geodynamics.org/benchmarks/trunk/cs/
doc/geodynamics.org/benchmarks/trunk/cs/index.rst
doc/geodynamics.org/benchmarks/trunk/geodyn/
doc/geodynamics.org/benchmarks/trunk/geodyn/index.rst
doc/geodynamics.org/benchmarks/trunk/index.rst
doc/geodynamics.org/benchmarks/trunk/long/
doc/geodynamics.org/benchmarks/trunk/long/circular-inclusion.rst
doc/geodynamics.org/benchmarks/trunk/long/divergence.rst
doc/geodynamics.org/benchmarks/trunk/long/drucker-prager.rst
doc/geodynamics.org/benchmarks/trunk/long/falling-sphere.rst
doc/geodynamics.org/benchmarks/trunk/long/geomod2004.rst
doc/geodynamics.org/benchmarks/trunk/long/geomod2008.rst
doc/geodynamics.org/benchmarks/trunk/long/index.rst
doc/geodynamics.org/benchmarks/trunk/long/relaxation-topography.rst
doc/geodynamics.org/benchmarks/trunk/magma/
doc/geodynamics.org/benchmarks/trunk/magma/index.txt
doc/geodynamics.org/benchmarks/trunk/magma/mckenzie-equations.txt
doc/geodynamics.org/benchmarks/trunk/magma/milestone1.txt
doc/geodynamics.org/benchmarks/trunk/magma/milestone2.txt
doc/geodynamics.org/benchmarks/trunk/magma/milestone3.txt
doc/geodynamics.org/benchmarks/trunk/magma/milestone4.txt
doc/geodynamics.org/benchmarks/trunk/magma/milestone5.txt
doc/geodynamics.org/benchmarks/trunk/magma/running-stgmadds.txt
doc/geodynamics.org/benchmarks/trunk/mc/
doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/
doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/index.rst
doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite1.rst
doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite2.rst
doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite3.rst
doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite4.rst
doc/geodynamics.org/benchmarks/trunk/mc/3d-spherical/
doc/geodynamics.org/benchmarks/trunk/mc/3d-spherical/index.rst
doc/geodynamics.org/benchmarks/trunk/mc/index.rst
doc/geodynamics.org/benchmarks/trunk/mc/notes-on-mantle-convection-benchmarks.rst
doc/geodynamics.org/benchmarks/trunk/seismo/
doc/geodynamics.org/benchmarks/trunk/seismo/index.rst
doc/geodynamics.org/benchmarks/trunk/short/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-landers/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-landers/description-landers.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/description-rs-nog.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/geofest-input/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/geofest-input/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/plots/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/plots/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/pylith-0.8-input/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/pylith-0.8-input/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/results/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/results/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/description-rs.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/results/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/results/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/description-ss.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/geofest-input/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/geofest-input/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/plots/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/plots/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/pylith-0.8-input/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/pylith-0.8-input/index.txt
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/results/
doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/results/index.txt
doc/geodynamics.org/benchmarks/trunk/short/index.txt
doc/geodynamics.org/benchmarks/trunk/short/overview.txt
doc/geodynamics.org/benchmarks/trunk/short/savage-prescott/
doc/geodynamics.org/benchmarks/trunk/short/savage-prescott/index.rst
doc/geodynamics.org/benchmarks/trunk/short/utilities/
doc/geodynamics.org/benchmarks/trunk/short/utilities/CUBITex.txt
doc/geodynamics.org/benchmarks/trunk/short/utilities/index.txt
Log:
Source files for http://geodynamics.org/cig/software/benchmarks/
The basic idea is to take advantage of the new documentation
generation tools that have been developed for Python, which
has recently been entirely migrated into reStructured Text
format (and using Sphinx to generate *.html files from these
*.rst and *.txt files):
* http://docutils.sourceforge.net/FAQ.html
* http://docutils.sourceforge.net/rst.html
* http://sphinx.pocoo.org/
Now, plone can already take reStructured Text as input, but
unfortunately it won't understand embedded latex equations
so we'll have to process those separately.
At any rate, once the bulk of the pages have been generated
and uploaded, the website will become the definitive version Corrections made in Plone (which are subsequently recorded
in the monolithic Zope Database) won't be reflected in these
text source files obviously, so we'll have to pay close attention
to the date-page-modified and date-page-generated timestamps.
Added: doc/geodynamics.org/benchmarks/trunk/cs/index.rst
===================================================================
Added: doc/geodynamics.org/benchmarks/trunk/geodyn/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/geodyn/index.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/geodyn/index.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,78 @@
+
+Historical Benchmark Cases
+==========================
+
+ Historically, there are two cases defined in the benchmark study published
+in the 2001 paper by Christensen et al. [6]. Case 0 is a benchmark of rotating
+non-magnetic convection. Case 1 is a dynamo with an insulating inner core
+co-rotating with the outer boundary. The regions outside the fluid shell are
+electrical insulators and the magnetic field on the boundaries matches with
+appropriate potential fields in the exterior that imply no external sources
+of the field.
+
+ In both cases the Ekman number is $E = 10^{3}$ and the Prandtl number is
+$Pr = 1$. The Rayleigh number is set to $Ra = 100000$. Note that the
+definition of the Rayleigh number differs from the one in the published
+cases [6] by a factor of Ekman number, i.e., $Ra=\frac{Ra}{E}$. The
+magnetic Prandtl number is zero in the non-magnetic convection case 0, and
+is $Pm = 5$ in case 1. The spherical harmonic expansion is truncated at
+degree $\ell_{max}=32$ and a four-fold symmetry is assumed in the
+longitudinal direction (`param.f` should be linked to `param32s4.f` when
+you compile MAG). The input parameter files are `par.bench0` for case 0,
+and `par.bench1` for case 1; both files reside in the `~/src` directory.
+
+ The output files of the benchmark cases are stored n the directory
+`~/bench-data/data_bench0` and `~/bench-data/data-bench1` respectively.
+In the following table we see the solutions from MAG agree with the
+benchmark suggested value with a small margin of difference. In both case 0
+and case 1, the values listed were obtained with low resolution and a
+relatively short run of MAG
+
+
++--------------+------------------------+------------+------------------------+-------------+
+| | Case 0 Suggested value | Mag Case 0 | Case 1 Suggested Value | Mag Case 1 |
++--------------+------------------------+------------+------------------------+-------------+
+| $E_{kin}$ | $58.348 \pm 0.050$ | 58.35 | $30.733 \pm 0.020$ | 30.72 |
+| $E_{mag}$ | | | $626.41 \pm 0.40$ | 627.15 |
+| $T$ | $0.42812 \pm 0.00012$ | | $0.37338 \pm 0.00040$ | |
+| $\mu_{\phi}$ | $-10.1571 \pm 0.0020$ | -10.80 | $-7.6250 \pm 0.0060$ | -7.84 |
+| $B_{\theta}$ | | | $-4.9289 \pm 0.0060$ | |
+| $\omega$ | $0.1824 \pm 0.0050$ | | $-3.1017 \pm 0.0040$ | |
++--------------+------------------------+------------+------------------------+-------------+
+
+
+Reversal Dynamo Case
+====================
+
+In this benchmark, we produce a magnetic field reversal using MAG. The
+input parameter in the source directory for this case is `~/src/par.Rev`.
+There is no longitudinal symmetry in this case, so when you compile MAG,
+use `param32s1.f` linking to `param.f`. The Ekman number is $E=0.02$, the
+Prandtl number is $Pr=1$ and the magnetic Prandtl number is $Pm=10$. The
+Rayleigh number is $Ra=12000$.
+
+Results and Discussions
+-----------------------
+
+This case was run on 32-bit and 64-bit Intel processors. Figure
+[fig:Field-Plot1] shows a plot of mean velocity Vrms, mean magnetic
+field Brms, the axial dipole and the dipole tilt on the outer boundary.
+It indicated a magnetic field reversal between times 25 and 30.
+Figure [fig:Field-Plot2] shows a longer run of MAG, where we see the
+magnetic field reversed again. At this time, the magnetic field had
+weakened substantially. In Figure [fig:The-pole], the top is the pole plot
+before the second field reversal and the bottom is the pole plot after the
+second field reversal.
+
+ image:: images/field-64.ps
+ Figure [fig:Field-Plot1]: Field Plot for Reversal Dynamo Case
+
+ image:: images/field-64-revR.ps
+ Figure [fig:Field-Plot2]: Field Plot for Reversal Dynamo Case (longer run)
+
+ image:: images/g1revR.ps
+ image:: images/g7revR.ps
+ Figure [fig:The-pole]: Magnetic Field Pole Plot. The top s the pole plot
+ at the beginning of the second field reversal; the bottom is the pole
+ plot at the end of the second field reversal.
+
Added: doc/geodynamics.org/benchmarks/trunk/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/index.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/index.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,17 @@
+
+Software Verification
+=====================
+
+Benchmark Efforts by Working Group
+----------------------------------
+
+ * "Short-Term Crustal Dynamics":http://geodynamics.org/cig/software/benchmarks/short/
+
+ * "Long-Term Tectonics":http://geodynamics.org/cig/software/benchmarks/long/
+
+ * "Mantle Convection":http://geodynamics.org/cig/software/benchmarks/mc/
+
+ * "Magma Migration":http://geodynamics.org/cig/software/benchmarks/magma/
+
+ * "Geodynamo":http://geodynamics.org/cig/software/benchmarks/geodyn/
+
Added: doc/geodynamics.org/benchmarks/trunk/long/circular-inclusion.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/circular-inclusion.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/circular-inclusion.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,75 @@
+
+Circular Inclusion
+==================
+
+ Schmid and Podladchikov [Clast] derived a simple analytic solution for
+ the pressure and velocity fields for a circular inclusion under simple
+ shear as in Figure [fig:inclusion-setup].
+
+ image:: images/inclusion_setup.eps
+ Figure [fig:inclusion-setup]
+ Schematic for the circular inclusion benchmark
+
+ The file `input/benchmarks/circular_inclusion/README` has instructions
+ on how to run this benchmark.
+
+ Because of the symmetry of the problem, we only have to solve over the
+ top-right quarter of the domain. For the velocity boundary conditions,
+ the analytic solution is a bit complicated. So we used the simple
+ relation $$v_x = -\dot{\epsilon}y, v_y=\dot{\epsilon}x,$$ for the
+ boundaries, where $\dot{\epsilon}$ is the magnitude of the shear and $x$
+ and $y$ are the coordinates. This induces an error of order $r_i^2/r^2$,
+ where $r_i=0.1$ is the radius of the inclusion, and $r$ is the radius. We
+ have the boundaries at 80 times the radius of the inclusion, giving an
+ error of about $0.01\%$, which is much smaller than the other errors we
+ were looking at. Just to make sure, we did runs with boundaries at 40
+ times the radius of the inclusion and got very similar results.
+
+ A characteristic of the analytic solution is that the pressure is zero
+ inside the inclusion, while outside it follows the relation
+ $$p_m=4\dot{\epsilon}\frac{\mu_m(\mu_i-\mu_m)}{\mu_i+\mu_m}\frac{r_i^2}{r^2}\cos(2\theta),$$
+ where $\mu_i=2$ is the viscosity of the inclusion and $\mu_m=1$ is the
+ viscosity of the background media. Many numerical codes that solve Stokes
+ flow (Eq. [eq:simple momentum conservation] and [eq:continuity]),
+ including Gale, assume that pressure, velocity, and viscosity are
+ continuous. The pressure discontinuity at the surface of the inclusion
+ violates that assumption, so the error tends to concentrate near the
+ surface of the inclusion.
+
+ Figure [fig:Pressure-inclusion] plots the error in the pressure along the
+ line $y=x/2$ for different resolutions. Inside the inclusion near the
+ surface, the pressure is consistently wrong. The pressure does not
+ converge with higher resolution, giving us a clue that the default
+ numerical scheme is not accurate.
+
+ image:: images/inclusion_r8_p.png
+ Figure [fig:Pressure-inclusion]
+ Pressure along the line $y=x/2$ for resolutions of $128 \times 128$
+ (blue), $256 \times 256$ (red), and $512 \times 512$ (black). The
+ inclusion has radius $r_i=0.1$. Note that the pressure should be zero
+ inside the inclusion, but the numerical solutions consistently
+ underestimate the pressure.
+
+ Outside the inclusion, the error is better behaved. Figure
+ [fig:Pressure-error] plots the error in the pressure along the line
+ $y=x/2$ outside the inclusion for different resolutions. While there are
+ still problems near the surface, away from the surface the solutions are
+ quite good. Figure [fig:Scaled-pressure-error] plots the error scaled
+ with resolution, and we can see that the error scales linearly with
+ resolution. This gives us confidence that, at least away from the
+ inclusion, the code is giving the right answer. This kind of result,
+ where the solution is bad close to the surface, but good otherwise, is
+ typical for numerical solutions of this problem [FD Stokes].
+
+ image:: images/inclusion_r8_p_error.png
+ Figure [fig:Pressure-error]
+ Error in the pressure outside the inclusion along the line $y=x/2$ for
+ resolutions of $128 \times 128$ (blue), $256 \times 256$ (red), and
+ $512 \times 512$ (black). The inclusion has radius $r_i=0.1$.
+
+ image:: images/inclusion_r8_p_scaled_error.png
+ Figure [fig:Scaled-pressure-error]
+ As in Figure [fig:Pressure-error], but with the error scaled with $h$.
+ So the medium-resolution error is multiplied by 2 and the
+ high-resolution error is multiplied by 4.
+
Added: doc/geodynamics.org/benchmarks/trunk/long/divergence.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/divergence.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/divergence.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,40 @@
+
+Divergence
+==========
+
+ This benchmark tests the implementation of the divergence term in the
+ equation [eq:divergence]. In 2D, a constant divergence is applied to a
+ square domain, and the velocity on the corners is set to enforce a
+ spreading from the center of the square. Figure [fig:Divergence_v_sri]
+ shows the velocity and strain rate invariant for a numerical solution.
+ For a constant divergence $d$, the analytic solution for this setup is
+ $$v_x = x \cdot d/2, v_y = y \cdot d/2$$
+
+ In 3D, the analytic solution is
+ $$v_x = x \cdot d/3, v_y = y \cdot d/3, v_z = z \cdot d/3$$
+
+ In both cases, the strain rate invariant equals $\sqrt{d/2}$. As shown in
+ Figure [fig:Divergence_2D_error], the main source of error in 2D comes
+ from inaccuracies in the solver. Figure [fig:Divergence_3D_error] paints
+ a different picture in 3D, where the main source of error comes from
+ having a finite number of particles.
+
+ image:: images/divergence_v.png
+ Figure [fig:Divergence_v_sri]
+ Velocity and Strain Rate Invariant solution for the 2D Divergence
+ benchmark. The variation in the strain rate invariant is uniformly
+ small.
+
+ image:: images/divergence_2D_error.eps
+ Figure [fig:Divergence_2D_error]
+ Maximum error in the strain rate invariant for the 2D Divergence
+ benchmark vs. tolerance in the linear solver. The resolution is kept at
+ $32 \times 32$, and the number of particles per cell is kept at 30.
+
+ image:: images/divergence_3D_error.eps
+ Figure [fig:Divergence_3D_error]
+ Maximum error in the strain rate invariant for the 3D Divergence
+ benchmark vs. the number of particles in each cell. The resolution is
+ kept at $16 \times 16 \times 16$, and the tolerance in the linear
+ solver is kept at $10^{-7}$.
+
Added: doc/geodynamics.org/benchmarks/trunk/long/drucker-prager.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/drucker-prager.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/drucker-prager.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,61 @@
+
+Drucker-Prager
+==============
+
+Analytic Treatment
+------------------
+
+ For the Drucker-Prager rehology in 2D, we can write the yielding relation
+ as $$\sigma_{ns}=\sigma_{nn}\tan\varphi+C,$$ where $\sigma_{ns}$ is the
+ shear stress perpendicular to the fault plane, $\sigma_{nn}$ is the shear
+ stress parallel to the fault plane, $\varphi$ is the internal angle of
+ friction, and $C$ is the cohesion. Decomposing this into principal
+ stresses $\sigma_{I}$, $\sigma_{II}$, and $\sigma_{III}$ gives
+ $$\sin(2\Theta)(\sigma_{I}-\sigma_{III})/2=\tan\varphi\left((\sigma_{I}+\sigma_{III})/2+\cos(2\Theta)(\sigma_{I}-\sigma_{III})/2\right)+C,$$
+ where $\Theta$ is the angle that the fault makes relative to the maximum
+ shear stress. Assuming that the fault forms where the shear stress
+ $\sigma_{I}-\sigma_{III}$ is a minimum, a little algebra gives us
+ $$\Theta=\pm\left(\frac{\pi}{4} + \frac{\varphi}{2}\right).$$
+
+ Using this, we can construct a simple plasticity experiment and make sure
+ that Gale gives the correct faulting angle.
+
+Model Setup
+-----------
+
+ We performed a shortening experiment as shown in Figure
+ [fig:Mohr-Coulomb-setup]. We only solve the Stokes equation and look at
+ the strain rate invariant to find incipient faults. We do not take any
+ time steps, removing any confounding effects they may cause.
+
+ image:: images/Mohr_Coulomb_setup.eps
+ Figure [fig:Mohr-Coulomb-setup]
+ The setup for the shortening experiment. The box is 1 unit on a side,
+ and the low viscosity region has a radius of 0.01 (its size is
+ exaggerated).
+
+Numerical Results
+-----------------
+
+ Figure [fig:Mohr-Coulomb-sri] shows the results for three different
+ resolutions for $\varphi=45^{\deg}$. There is not much difference between
+ the medium ($256 \times 256$) and high ($512 \times 512$) results,
+ suggesting that we have sufficient resolution. Figure
+ [fig:Mohr-Coulomb-comparison] shows a plot of the numerical vs. analytic
+ results for several different angles for medium resolution. This gives us
+ confidence that, at least in compression(sp?) in 2D, our Drucker-Prager
+ implementation gives the correct results.
+
+ image:: images/Mohr_coulomb_resolutions.png
+ Figure [fig:Mohr-Coulomb-sri]
+ Strain rate invariant for the shortening experiment where
+ $\varphi=45^{\deg}$ with three different resolutions:
+ $128 \times 128$, $256 \times 256$, $512 \times 512$.
+ Any differences between the medium and high resolution runs are swamped
+ by uncertainties in determining the overall angle of faulting.
+
+ image:: images/mohr_coulomb_angles.eps
+ Figure [fig:Mohr-Coulomb-comparison]
+ Numerical vs. analytic results for fault angles as a function of
+ internal angle of friction.
+
Added: doc/geodynamics.org/benchmarks/trunk/long/falling-sphere.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/falling-sphere.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/falling-sphere.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,74 @@
+
+Falling Sphere
+==============
+
+ This benchmark simulates a rigid sphere falling through a cylinder filled
+ with a viscous medium as in Figure [fig:Sphere-Cylinder]
+
+ image:: sphere_cylinder.eps
+ Figure [fig:Sphere-Cylinder]
+ Schematic of a Sphere falling through a Cylinder.
+
+ The file `input/benchmarks/falling_sphere/README` has instructions on
+ running this benchmark. In an infinitely large cylinder, the analytic
+ solution for the drag on a sphere is $$F=6\pi\eta r u,$$ where $\eta$ is
+ the viscosity of the medium, $r$ is the radius of the sphere, and $u$ is
+ the velocity of the sphere. Conversely, the buoyancy force is
+ $$F=\frac{4}{3}\pi{r^3}g\delta\rho,$$ where $g$ is the gravitational
+ constant and $\delta\rho$ is the density difference between the sphere
+ and the medium. Balancing these two forces and solving for the velocity
+ gives $$u = \frac{2}{9}{r^2}g\delta\rho / \eta.$$
+
+ Seetting $g=1$, $r=1$, $\delta\rho=0.01$, and $\eta=1$ gives a velocity
+ of $$u=0.00222.$$
+
+ In our case, we simulate a rigid sphere with a high viscosity sphere.
+ This allows some internal circulation within the sphere, and so the
+ expression for the velocity becomes [Landau & Lifschitz],
+ $$u = \frac{1}{3}\frac{r^2{g}\delta\rho}{\eta}\frac{\eta+\eta'}{\eta+\frac{3}{2}\eta'},$$
+ where $\eta'$ is the viscosity of the sphere. Four our case, the
+ background medium's viscosity is 1 and the sphere's viscosity is 100, so
+ the correction is about $1\%$. This turns out to be smaller than other
+ effects for the cases we ran.
+
+ when the boundaries are not infinitely far away, we can expand the
+ solution in terms of the ratio of the radius of the sphere ($r$) to the
+ radius of the cylinder ($R$). One solution by Habermann [Stokes Sphere]
+ gives a drag force of
+ $$F_H=6\pi\eta{ru}\frac{1-0.75857\left(\frac{r}{R}\right)^5}{1+f_H\left(\frac{r}{R}\right)},$$
+ where
+ $$f_H\left(\frac{r}{R}\right)=-2.1050(r/R)+2.0865(r/R)^3-1.7068(r/R)^5+0.72603(r/R)^6.$$
+
+ For our case with $r=1$, $R=4$, this gives a velocity of
+ $$u=1.122747319\cdot{10^{-3}},$$ which agrees closely with the result
+ from Habermann.
+
+ Another possible artifact is that we do not simulate an infinitely long
+ cylinder. This turns out to be a small effect. We use a cylinder with a
+ height of 8, and place the sphere halfway down. We did runs where the
+ cylinder was twice as tall, and the results were essentially unchanged.
+
+ The errors in the computed velocity compared to the Faxen solution are
+ plotted in Figure [fig:Error-in-velocity]. These were done with
+ resolutions of $8 \times 16 \times 8$, $16 \times 32 \times 16$, and
+ $64 \times 128 \times 64$, corresponding to grid sizes ($h$) of $0.5$,
+ $0.25$, $0.125$, and $0.0625$. Because of the symmetries of the problem
+ we only have to simulate a quarter of the domain. As we increase the
+ resolution (decrease $h$), the error decreases. Since we are simulating a
+ high viscosity sphere rather than a completely rigid sphere, the velocity
+ inside the sphere is not uniform. The error bars indicate the variation
+ in velocity across the sphere.
+
+ image:: images/Sphere_Error.eps
+ Figure [fig:Error-in-velocity]
+ Error in computed velocity vs. resolution.
+
+ Scaling the error with resolution gives Figure [fig:Scaled-error-velocity].
+ The error scales linearly with resolution, giving us confidence that we
+ can accurately solve this problem.
+
+ image:: images/Sphere_Scaled_Error.eps
+ Figure [fig:Scaled-error-velocity]
+ As in figure [fig:Error-in-velocity], but with the error scaled with $h$.
+ So the higher resolution errors are multiplied by 2, 4, and 8.
+
Added: doc/geodynamics.org/benchmarks/trunk/long/geomod2004.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/geomod2004.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/geomod2004.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,107 @@
+
+===========
+Geomod 2004
+===========
+
+ Two benchmarks were created to validate numerical codes against analog
+ sandbox experiments [Buiter et al Numerical Sandbox]: one benchmark
+ simulates extension, and the other simulates shortening. A number of
+ investigators with different codes ran these benchmarks, giving us a good
+ standard against which to compare.
+
+Extension
+=========
+
+ This benchmark simulates a sandbox being extended as in Figure
+ [fig:Extension-model-setup]. The right side and half of the bottom are
+ translated to the right. This creates a velocity discontinuity at the
+ center which is the initial seed for localization. Gale's implementation
+ of this benchmark is in `input/benchmarks/extension.xml`.
+
+ Like half of the codes in the benchmark, boundary friction was not
+ included. Rather, the material is held fixed to the bottom boundary, and
+ the velocity discontinuity is smoothed over 0.2 cm. In addition, the
+ exact background viscosity is not prescribed by the benchmark. We have
+ used $10^{12}\ Pa \cdot s$, the same as used in the I2ELVIS calculations.
+ This value is somewhere in the middle of the range of values used in the
+ calculations for other codes.
+
+ image:: images/Extension_setup.png
+ Figure [fig:Extension-model-setup]
+ Extension model setup. Reproduced, with permission, from Buiter et al.
+ [Buiter et al Numerical Sandbox]
+
+ A comparison against the other codes is in Figure
+ [fig:Comparison-extension]. While it is difficult to perform exact
+ comparisons, Gale produces similar fault patterns. It is interesting to
+ note that this qualitative comparison holds true even though the code is
+ not exactly convergent. For example, Figure [fig:extension-convergence]
+ shows a run with varying resolution. Changing the resolution alters the
+ details, but, after a certain minimum resolution, does not change the
+ character of the solution.
+
+ image:: images/Extension_comparision.png (sp?)
+ Figure [fig:Comparison-extension]
+ Strain rate invariant for the numerical extension models after 5 cm of
+ extension. The resolutions of the various models are:
+ I2ELVIS: $400 \times 75$
+ LAPEX-2D: $301 \times 71$
+ Microfem: $201 \times 61$
+ SloMo: $401 \times 71$
+ Sopale: $401 \times 71$
+ Gale: $1024 \times 128$
+ Uper images reproduced, with permission, from Buiter et al.
+ [Buiter et al Numerical Sandbox].
+
+ image:: images/extension_sensitivity.pdf
+ Figure [fig:extension-convergence]
+ Strain rate invariant for the extension model after 5 cm of extension
+ for four different resolutions: (a) 128x16, (b) 256x32, (c) 512x64,
+ and (d) 1024x128.
+
+
+Shortening
+==========
+
+ This benchmark simulates a sandbox being shortened as in Figure
+ [fig:Shortening-model-setup]. The right side is moved to the left,
+ creating a velocity discontinuity at the bottom right corner. Gale's
+ implementation of this benchmark is in `input/benchmarks/shortening.xml`.
+
+ Like most of the codes in the benchmark, we applied diffusion (diffusion
+ coefficient $10^{-6}\ m^2\ s^{-1}$) to the top surface to smooth steep
+ slope angles. This lies within the range used in calculations used by the
+ other codes (LAPEX-2D: $10^{-6}\ m^2\ s^{-1}$, Microfem: unspecified,
+ Sopale: $10^{-9}\ m^2\ s^{-1}$). Again, the exact background viscosity is
+ not prescribed by the benchmark, so we have used $10^{12}\ Pa \cdot s$.
+
+ image:: images/Shortening_setup.png
+ Figure [fig:Shortening-model-setup]
+ Shortening model setup. Reproduced, with permission, from Buiter et al.
+ [Buiter et al Numerical Sandbox].
+
+ A comparison against the other codes' calculations at 14 cm of cumulative
+ shortening is in Figure [fig:shortening-compare]. There is more variance
+ among the different codes, so it is difficult to tell whether Gale's
+ behavior agrees with the other codes. Figure [fig:shortening-convergence]
+ shows a run with a few different resolutions, and even there we see
+ marked differences in behavior as we increase resolution.
+
+ image:: images/Shortening_comparison.png
+ Figure [fig:shortening-compare]
+ Strain rate invariant for the numerical shortening models after 14 cm
+ of shortening. The resolutions of the various models are:
+ I2ELVIS: 900 x 75,
+ LAPEX-2D: 351 x 71,
+ Microfem: 201 x 36,
+ Sopale: 401 x 71,
+ Gale: 512 x 128.
+ The upper portion of the figure is reproduced, with permission, from
+ Buiter et al. [Buiter et al Numerical Sandbox].
+
+ image:: images/shortening_sensitivity.pdf
+ Figure [fig:shortening-convergence]
+ Strain rate invariant for the shortening model after 14 cm of
+ shortening for three different resolutions: (a) 128 x 32,
+ (b) 256 x 64, (c) 512 x 128.
+
Added: doc/geodynamics.org/benchmarks/trunk/long/geomod2008.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/geomod2008.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/geomod2008.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,187 @@
+
+===========
+Geomod 2008
+===========
+
+ Using the lessons learned from the Geomod 2004 benchmarks, new
+ benchmarks were created that would make it easier to compare numerical
+ experiments with each other and with analog experiments [Geomod 2008].
+
+Stable Wedge
+============
+
+ This benchmark simulates a wall pushing a wedge as in Figure
+ [fig:Wedge_setup]. There is an analytic solution [Dahlen Wedge] which
+ details what the friction on the bottom and sides should be to provide
+ enough resistance so that the wedge does not collapse under its own
+ weight, but not so much as to cause any internal deformation as it
+ slides. The derivation of the solution assumes that the friction along
+ the sides has no cohesion. So the force at the tip will go to zero as the
+ thickness of the material goes to zero. However, analog experiments
+ suggest a finite cohesion, so this benchmark specifies a boundary
+ cohesion.
+
+ We modeled the wedge using a relatively low viscosity ($1\ Pa \cdot s$)
+ air layer on top. This low viscosity region does not, for the most part,
+ affect the dynamics.
+
+ We modeled boundary friction by first fixing the sand to the boundary. We
+ then modify the material properties in the element next to the boundary
+ so that it provides the correct resistance. So in the bulk, the sand's
+ internal angle of friction is $36^{\deg}$ weakening to $31^{\deg}$, while
+ in the element at the boundary the internal angle of friction is
+ $16^{\deg}$ weakening to $14^{\deg}$. Similarly, in the bulk, the
+ cohesion is $10\ Pa$, while at the boundary the cohesion is $10\ Pa$
+ weakening to $0.01\ Pa$. If we do not weaken the cohesion, when we try to
+ model an unstable wedge by reducing the internal angle of friction, the
+ wedge never collapses on itself.
+
+ Figure [fig:Stable_sri] shows the strain rate invariant after the wall
+ has moved 4 cm, and Figure [fig:Stable_particles] shows the particles.
+ The bulk translates with almost no deformation, although, as expected,
+ the tip deforms. The odd structures at the tip are below the grid
+ resolution. Figure [fig:Stable_unstable] shows a simulation when we
+ reduce the boundary friction to $1^{\deg}$. The wedge quickly becomes
+ unstable and collapses.
+
+ image:: images/Geomod2008_wedge_setup.eps
+ Figure [fig:Wedge_setup]
+ Setup for the stable wedge benchmark. Image courtesy of Susanne Buiter.
+
+ image:: images/Stable_wedge_sri.png
+ Figure [fig:Stable_sri]
+ Strain rate invariant for the stable wedge benchmark within the wedge.
+ Outside the wedge, the strain rates are large because of the air's low
+ viscosity. The resolution is 512.128, and the wedge has translated
+ 4 cm.
+
+ image:: images/Stable_wedge_particles.png
+ Figure [fig:Stable_particles]
+ Material particles for the stable wedge benchmark. The deformation at
+ the tip is caused by a finite boundary cohesion, although the actual
+ structure is not resolved. The resolution is 512x128, and the wedge has
+ translated 4 cm.
+
+ image:: images/Stable_wedge_unstable.png
+ Figure [fig:Stable_unstable]
+ Strain rate invariant and velocity arrows for the stable wedge
+ benchmark, but with the friction angle reduced to $1^{\deg}$. Note that
+ the strain rates are much higher than in Figure [fig:Stable_sri]. The
+ wedge collapses almost immediately. The resolution is 512x128, and the
+ wedge has translated 0.17 cm.
+
+
+Unstable Shortening
+===================
+
+ This benchmark simulates a wall pushing against a wall of sand as in
+ Figure [fig:Unstable-setup]. There are three layers of sand, with the
+ middle layer being a little heavier and sticking a little more to the
+ boundary. Otherwise it is identical. Figures
+ [fig:unstable_sri_128],
+ [fig:unstable_sri_256],
+ [fig:unstable_sri_512],
+ [fig:unstable_particles_128],
+ [fig:unstable_particles_256], and
+ [fig:unstable_particles_512] show results at different times and
+ different resolutions.
+
+ image:: images/Geomod2008_unstable_setup.eps
+ Figure [fig:Unstable-setup]
+ Setup for the unstable shortening benchmark.
+ Image courtesy of Susanne Buiter.
+
+ image:: images/Geomod2008_unstable_sri128x32.png
+ Figure [fig:unstable_sri_128]
+ Strain rate invariant for the unstable shortening benchmark with a
+ resolution of 128x32. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_unstable_sri256x64.png
+ Figure [fig:unstable_sri_256]
+ Strain rate invariant for the unstable shortening benchmark with a
+ resolution of 256x64. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_unstable_sri512x128.png
+ Figure [fig:unstable_sri_512]
+ Strain rate invariant for the unstable shortening benchmark with a
+ resolution of 512x128. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_unstable_particles128x32.png
+ Figure [fig:unstable_particles_128]
+ Material particles for the unstable shortening benchmark with a
+ resolution of 128x32. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_unstable_particles256x64.png
+ Figure [fig:unstable_particles_256]
+ Material particles for the unstable shortening benchmark with a
+ resolution of 256x64. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_unstable_particles512x128.png
+ Figure [fig:unstable_particles_512]
+ Material particles for the unstable shortening benchmark with a
+ resolution of 512x128. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+
+Brittle Shortening
+==================
+
+ This benchmark is very similar to unstable shortening. The only
+ difference is that part of the bottom is also moving along as shown in
+ Figure [fig:Brittle_setup]. This causes the deformation to start from
+ inside the sand box, rather than along the walls. F
+ Figures
+ [fig:brittle_sri_128],
+ [fig:brittle_sri_256],
+ [fig:brittle_sri_512],
+ [fig:brittle_particles_128],
+ [fig:brittle_particles_256], and
+ [fig:brittle_particles_512] show results at different times and
+ different resolutions.
+
+ image:: images/Geomod2008_brittle_setup.eps
+ Figure [fig:Brittle_setup]
+ Setup for the brittle shortening benchmark.
+ Image courtesy of Susanne Buiter.
+
+ image:: images/Geomod2008_brittle_sri128x32.png
+ Figure [fig:brittle_sri_128]
+ Strain rate invariant for the brittle shortening benchmark with a
+ resolution of 128x32. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_brittle_sri256x64.png
+ Figure [fig:brittle_sri_256]
+ Strain rate invariant for the brittle shortening benchmark with a
+ resolution of 256x64. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_brittle_sri512x128.png
+ Figure [fig:brittle_sri_512]
+ Strain rate invariant for the brittle shortening benchmark with a
+ resolution of 512x128. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_brittle_particles128x32.png
+ Figure [fig:brittle_particles_128]
+ Material particles for the brittle shortening benchmark with a
+ resolution of 128x32. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_brittle_particles256x64.png
+ Figure [fig:brittle_particles_256]
+ Material particles for the brittle shortening benchmark with a
+ resolution of 256x64. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
+ image:: images/Geomod2008_brittle_particles512x128.png
+ Figure [fig:brittle_particles_512]
+ Material particles for the brittle shortening benchmark with a
+ resolution of 512x128. The snapshots are taken at 0, 2.5, 5, 7.5, and
+ 10 cm of shortening.
+
Added: doc/geodynamics.org/benchmarks/trunk/long/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/index.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/index.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,38 @@
+
+Benchmarks
+==========
+
+* Falling Sphere
+
+* Circular Inclusion
+
+* Relaxation of Topography
+
+* Divergence
+
+* Drucker-Prager
+
+* Geomod 2004
+
+ * Extension
+
+ * Shortening
+
+* Geomod 2008
+
+ * Stable Wedge
+
+ * Unstable Shortening
+
+ * Brittle Shortening
+
+
+
+Links
+
+* "Four Gale Tutorials":http://geodynamics.org/cig/software/packages/long/gale/tutorials
+
+* "Original Work Area (currently empty)":http://geodynamics.org/cig/workinggroups/long/workarea
+
+
+
Added: doc/geodynamics.org/benchmarks/trunk/long/relaxation-topography.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/long/relaxation-topography.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/long/relaxation-topography.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,44 @@
+
+Relaxation of Topography
+========================
+
+ Given an infinitely deep purely viscous medium with an infinitesimal
+ sinusoidal height profile, the topography will decay exponentially
+ with the timescale [Folds] $$t_r = \frac{4\pi\eta}{gL},$$ where
+ $\eta$ is the viscosity, $g$ is the gravitational constant, and $L$ is
+ the wavelength of the initial sinusoid.
+
+ In our case, we simulate a medium with finite depth and finite height.
+ The internal fields decay exponentially with depth with a length scale of
+ $L/{2\pi}$. The error in the solution due to a finite height is of order
+ $(2\pi{A}/L)^2$, where $A$ is the amplitude of the sinusoid. We use $L=1$
+ and $A=0.01$, giving errors of order $0.02\%$ and $0.4\%$.
+
+ The file `input/benchmarks/sinusoid/README` explains how to run this
+ benchmark. Figure [fig:Strain-topo] shows the results of a low-resolution
+ run. Even this run is not particularly small ($128 \times 256$), because
+ we need fairly high resolution to be able to accurately resolve the small
+ ($1\%$) height difference. Also note that we use symmetry to only
+ simulate half of the wavelength.
+
+ image:: images/Paraview_topography.png
+ Figure [Strain-topo]
+ Strain rate and velocities for a sinusoidal topography relaxing under
+ gravity.
+
+ Running the code with multiple resolutions and measuring the error in the
+ height in the trough gives Figure [fig:topo-error]. Scaling the error
+ with resolution gives Figure [fig:scaled-topo-error]. The error decreases
+ linearly with increasing resolution, giving us confidence in our ability
+ to accurately track topography.
+
+ image:: images/topo_error.eps
+ Figure [fig:topo-error]
+ Error in the height at the trough
+
+ image:: images/topo_scaled_error.eps
+ Figure [fig:scaled-topo-error]
+ As in Figure [fig:topo-error], but with the error scaled with $h$.
+ So the medium-resolution error is multiplied by 2 and the
+ high-resolution error is multiplied by 4.
+
Added: doc/geodynamics.org/benchmarks/trunk/magma/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,171 @@
+Benchmarks
+
+* An Introduction and Tutorial to the "McKenzie Equations" for Magma Migration
+ [last modified 2007-01-15] --> mckenzie-equations.txt
+
+ A new formulation for the equations of magma migration in viscous materials
+ as originally derived by McKenzie is presented, as well as a set of
+ well-understood special case problems that form a useful benchmark-suite
+ for developing and testing new codes.
+
+
+* Running stgMADDs Benchmarks
+ [last modified 2009-04-02] --> running-stgmadds.txt
+
+ The Magma Development team has finished the alpha release of the
+ Magma Dynamics Demonstration Suite (MADDs). The initial code implements
+ the zero porosity / no melting magma benchmark for mid-ocean ridge
+ solid flows in 2D and 3D built on the Underworld framework. The purpose
+ of this code is principally to validate accurate pressure solvers for
+ Stokes flow in current CIG supported software. The stgMADDs source
+ code is available in CIG's Mercurial Repository (geodynamics.org/hg).
+
+-------------------------------------------------------------------------------
+
+* Milestone1 Results and Analysis
+ [last modified 2008-02-08] --> milestone1.txt
+
+ Details how to run the first milestone of the MADDs project in 2D and 3D
+ and provides some results of those simulations. It also gives the rates
+ of convergence of the pressure gradient solutions as the resolution
+ is increased.
+
+ 2D Ridge Model
+ Velocity, pressure, and pressure gradients solutions and L2 errors
+ for a 2D ridge model with 120 x 60 elements.
+
+ 3D Ridge Model
+ Velocity, pressure, and pressure gradients solutions and L2 error fields
+ for 3D ridge model.
+
+ Global Pressure Gradient Errors for 2D Ridge Model
+ Normalized global L2 errors.
+
+ Global Pressure Gradient Errors for 3D Ridge Model
+ Global normalized L2 pressure gradient errors at varying resolutions.
+
+-------------------------------------------------------------------------------
+
+* Milestone2 Results and Analysis
+ [last modified 2008-02-08] --> milestone2.txt
+
+ Details the results of the Milestone2 simulations and analyzes the accuracy
+ of the advection scheme.
+
+ Gaussian Porosity Field Advection
+ Advection of Gaussian porosity field as a Stokes equation force term.
+ The lower density porosity region rises due to gravity.
+
+ Ridge Model with Gaussian Porosity Field
+ Stokes flow with 2D ridge model boundary conditions and Gaussian
+ porosity initial distribution, driven by a porosity dependent
+ force term.
+
+ Semi Lagrangian Advection Scheme Test - Step Function
+ Diagonal step function initial distribution subjected to a
+ shearing velocity field.
+
+ Semi Lagrangian Advection Scheme Test - Gaussian Distribution
+ Gaussian initial distribution subjected to a shearing velocity field.
+
+ Error Convergence for Advection Scheme - Step Function IC
+ Normalized global L2 errors for semi Lagrangian advection scheme
+ with a diagonal step function initial condition as a function
+ of resolution.
+
+ Error Convergence for Advection Scheme - Gaussian IC
+ Normalized global L2 errors for semi Lagrangian advection scheme
+ with Gaussian initial distribution as a function of resolution.
+
+-------------------------------------------------------------------------------
+
+* Milestone3 Results
+ [last modified 2008-02-08] --> milestone3.txt
+
+ Details the results for the third milestone, in which the melt velocity
+ was determined given the existing solid velocity and pressure fields.
+
+ Melt Model - 2D Ridge with Constant Porosity
+ Solid and melt velocity, pressure, and pressure gradient fields
+ for 2D ridge model with constant porosity. Melt velocity magnitudes
+ are significantly larger near the point of discontinuity due to
+ their proportionality to the pressure gradients, which are largest
+ at these points.
+
+ Melt Model - Gaussian Porosity Driven Flow
+ Solid and melt velocity, pressure, and pressure gradient fields
+ for Stokes flow driven by a Gaussian initial porosity distribution.
+
+-------------------------------------------------------------------------------
+
+* Milestone4 Results and Analysis
+ [last modified 2008-09-28] --> milestone4.txt
+
+ Discussion of the system being modeled, and details of how to run the
+ model with different initial conditions in 2D and 3D.
+
+ 2D Solitary Wave
+ A 2D solitary wave with a wave speed of 7 rising through a solid
+ with a constant speed of -2. The wave shows no visible diffusive
+ behavior.
+
+ Noisy 1D Solitary Wave Initial Condition
+ Initial condition of a vertically changing 1D solitary wave with
+ a certain amount of introduced noise, which allows 2D solitary
+ waves to emerge over time.
+
+ Emerging 2D Solitary Waves
+ Solitary waves emerging from a noisy 1D solitary wave initial condition.
+
+ Emergent 2D Solitary Waves
+ Solitary waves having emerged from a noisy 1D solitary wave initial distribution.
+
+-------------------------------------------------------------------------------
+
+* Milestone5 Results and Analysis
+ [last modified 2009-03-26] --> milestone5.txt
+
+ Results and analysis for the isoviscous McKenzie equations (with melting)
+ driven by a corner flow velocity BC.
+
+ Isoviscous McKenzie System with Corner Flow BC - 1
+ After 1 time step
+
+ Isoviscous McKenzie System with Corner Flow BC - 50
+ After 50 time steps
+
+ Isoviscous McKenzie System with Corner Flow BC - 3200
+ After 3200 time steps
+
+ Velocity - x component
+ x-component of the velocity field for the 3D isoviscous McKenzie model
+ with ridge BCs at time step 150.
+
+ Velocity - y component
+ y-component of the velocity field for the 3D isoviscous McKenzie model
+ with ridge BCs at time step 150.
+
+ Porosity
+ Porosity field for the 3D isoviscous McKenzie model with ridge BCs
+ at time step 150.
+
+ Compaction pressure
+ Compaction pressure due to compressibility of the solid phase for the
+ 3D isoviscous McKenzie model with ridge BCs at time step 150.
+
+ Melt fraction
+ Melt fraction field representing the melt to solid phase of the
+ 3D isoviscous McKenzie model with ridge BCs at time step 150.
+
+ Melt velocity - x component
+ x-component of the melt velocity field for the 3D isoviscous McKenzie model
+ with ridge BCs at time step 150.
+
+ Melt velocity - y component
+ y-component of the melt velocity field for the 3D isoviscous McKenzie model
+ with ridge BCs at time step 150.
+
+-------------------------------------------------------------------------------
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/
Added: doc/geodynamics.org/benchmarks/trunk/magma/mckenzie-equations.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/mckenzie-equations.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/mckenzie-equations.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,11 @@
+
+An Introduction and Tutorial to the "McKenzie Equations" for Magma Migration
+ [last modified 2007-01-15]
+
+ A new formulation for the equations of magma migration in viscous materials
+ as originally derived by McKenzie is presented, as well as a set of well-understood
+ special case problems that form a useful benchmark-suite for developing and
+ testing new codes.
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/McKenzieIntroBenchmarks.pdf/view
Added: doc/geodynamics.org/benchmarks/trunk/magma/milestone1.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/milestone1.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/milestone1.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,60 @@
+Plone Metadata
+
+ milestone1results
+
+ Milestone 1 Results and Analysis
+
+ Details how to run the first milestone of the MADDs project in 2D and 3D
+ and provides some results of these simulations. It also gives the rates
+ of convergence of the pressure gradient solutions as the resolution is
+ increased.
+
+Running the code
+
+ In order to run the simulations for milestone 1 of the MADDs project
+ (in 2D), first::
+
+ cd Magma/Models/Milestone1/Ridge2D_Quadratic
+
+ Then make a symbolic link to the executable binary (assuming the code
+ has been successfully built)::
+
+ ln -s ../../../../build/bin/StGermain .
+
+ The simulation may then be run (in parallel), passing the respective XML
+ file as input::
+
+ mpiexec -np <#-of-procs> ./StGermain Ridge2D.xml
+
+ Alternatively, the 3D simulation may be run as::
+
+ cd Magma/Models/Milestone1/Ridge3D_Quadratic
+ ln -s ../../../../build/bin/StGermain .
+ mpiexec -np <#-of-procs> ./StGermain Ridge3D.xml
+
+
+Simulation results and error convergence
+
+ These simulations will produce graphical output of the velocity,
+ pressure, and pressure gradient solutions, as well as the analytic
+ reference solutions and the element-wise normalized L2 error fields for
+ the pressure and pressure gradients, as shown below. It will also
+ generate text files to the output directory giving the node-wise results
+ for the respective fields.
+
+ As the resolution is increased, the normalized global L2 errors are
+ observed to decrease. This decrease is approximately linear for the 2D
+ ridge mode and slightly poorer for the 3D model. Graphs detailing the
+ global errors as a function of resolution are given below.
+
+
+Related Item(s)
+
+ 2D Ridge Model
+ 3D Ridge Model
+ Global Pressure Gradient Errors for 2D Ridge Model
+ Global Pressure Gradient Errors for 3D Ridge Model
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/milestone1results/
Added: doc/geodynamics.org/benchmarks/trunk/magma/milestone2.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/milestone2.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/milestone2.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,87 @@
+Plone Metadata
+
+ milestone2results
+
+ Milestone 2 Results and Analysis
+
+ Details the results of the Milestone 2 simulations and analyzes the
+ accuracy of the advection scheme.
+
+Model Descriptions
+
+ The second milestone contains several different models. In each a
+ porosity distribution has been supplied, which acts as a force term on
+ the stokes equation. The porosity is advected according to a semi
+ Lagrangian scheme, with no natural diffusion. Details of these models are
+ given below.
+
+Simple Gaussian porosity field
+
+ The first model is a simple demonstration of how the porosity dependent
+ force term drives the porosity distribution up through the domain. The
+ domain is subjected to free slip boundary conditions, which also distort
+ the porosity distribution as it is advected. This model is housed in
+ the directory::
+
+ Magma/Models/Milestone2/ZPNM_SemiLagrangianPorosity/
+
+ with the XML file to be passed at run time being the 'Porridge.xml' file.
+
+Ridge model with Gaussian porosity distribution
+
+ The second model is an extension of the first milestone, that is a ridge
+ model with the intial boundary conditions loaded from a reference
+ solution. However in this case a Gaussian porosity distribution has been
+ added. The porosity distribution is distorted as it moves up through the
+ domain in accordance with the velocity field generated from the boundary
+ conditions. This model can be found in::
+
+ Magma/Models/Milestone2/RidgeModelWithGaussianPorosity/
+
+ with the XML file to be passed being 'Ridge2D.xml'
+
+
+Validation of the advection scheme
+
+ As well as the models, a test is also supplied for validating the
+ accuracy of the semi Lagrangian advection scheme. This involves an
+ initial porosity distribution (either Gaussian or diagonal line step
+ function, as determined from the XML), which is subjected to a static
+ shearing velocity field. The porosity distribution is subjected to the
+ velocity field for a finite number of time steps (as determined from the
+ XML input file), before the velocity field is reversed and the
+ distribution advected back for the same number of time steps. The
+ normalized global L2 error between the initial and final distributions is
+ then calculated. This test is housed in the directory::
+
+ Magma/Models/Milestone2/tests/
+
+ and may be run using the input file 'testSemiLagrangianAdvection.xml'.
+
+ Running this simulation at varying resolutions, the convergence of the
+ errors for the advection scheme were determined using both the Gaussian
+ distribution and diagonal step function as initial conditions. The errors
+ were recorded for schemes which used a cubid method as well as a
+ quadratic method based on the element shape functions for interpolating
+ the value at the end point of the characteristic. As can be observed from
+ the convergence error plots below, the quadratic and cubic interpolation
+ method schemes converged at comparable (less than linear) rates using the
+ step function initial condition, with an improvement in those results
+ obtained using the cubic interpolation method. When the Gaussian initial
+ distribution was applied (which was easier to solve accurately on account
+ of the smoother gradients involved), both interpolation methods converged
+ at a rate much better than linear, with the cubic interpolation method
+ proving far superior to the quadratic method.
+
+Related Item(s)
+
+ Gaussian Porosity Field Advection
+ Ridge Model With Gaussian Porosity Field
+ Semi Lagrangian Advection Scheme Test - Step Function
+ Semi Lagrangian Advection Scheme Test - Gaussian Distribution
+ Error Convergence for Advection Scheme - Step Function IC
+ Error Convergence for Advection Scheme - Gaussian IC
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/milestone2results/
Added: doc/geodynamics.org/benchmarks/trunk/magma/milestone3.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/milestone3.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/milestone3.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,39 @@
+Plone Metadata
+
+ milestone3results
+
+ Milestone 3 Results
+
+ Details the results for the third milestone, in which melt velocity was
+ determined given the existing solid velocity and pressure fields.
+
+Model descriptions
+
+ Two different models were implemented for which the melt velocity was
+ determined. The first of these was an extension of the ridge model
+ implemented in Milestone 1, but with a constant porosity field, such that
+ the solution was static in time. This model can be found in::
+
+ /Magma/Models/Milestone3/Ridge2D_Field_BasedConstantPorosity
+
+ The second model was an extension of the porosity driven Stokes flow with
+ a Gaussian intial distribution implemented in Milestone 2. This model
+ resides at::
+
+ /Magma/Models/Milestone3/FieldBasedPorosityDrivenFlow2D
+
+ Since the melt velocity is decoupled from the McKenzie equations, it is
+ relatively simple to calculate, provided that the pressure and solid
+ velocity fields have already been accurately determined. As such no tests
+ were applied to validate its accuracy, however qualitatively their
+ behavior is observed to be correct.
+
+
+Related Item(s)
+
+ Melt Model - Gaussian Porosity Driven Flow
+ Melt Model - 2D Ridge with Constant Porosity
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/milestone3results/
Added: doc/geodynamics.org/benchmarks/trunk/magma/milestone4.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/milestone4.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/milestone4.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,71 @@
+Plone Metadata
+
+ milestone4results
+
+ Milestone 4 - Results and Analysis
+
+ Discussion of the system being modeled, and details of how to run the
+ model with different initial conditions in 2D and 3D.
+
+Problem Description
+
+ This milestone solves a porosity pressure system which involves the
+ coupling of a Darcy flow to describe the compressibility of the permeable
+ solid matrix and a time dependent advection equation for the porosity
+ field. Together these equations allow for non-linear dispersive porosity
+ waves. Given an initial porosity distribution which is itself a porosity
+ wave, this wave should advect at a speed determined by the amplitude and
+ the power used to determine the permeability from the porosity (as well
+ as the velocity of the background solid), without any diffusion. If the
+ initial porosity distribution is not itself a solitary porosity wave,
+ with time these should emerge from the porosity field.
+
+Running the Simulations
+
+ In order to run the solitary waves model, (in 2D) first::
+
+ cd Magma/Models/Milestone4/SolitaryWaves2D
+
+ Then make a symbolic link to the executable binary as::
+
+ ln -s ../../../../build/bin/StGermain .
+
+ The simulation may then be run in parallel as::
+
+ mpirun -np <#-of-procs> ./StGermain SolitaryWaves.xml
+
+ This will run the code with the default initial porosity distribution of
+ a solitary wave with a wave speed of 7 and a porosity exponent of 3. The
+ background solid velocity has been set in the file 'VelocityField.xml' as
+ -2, such that the wave should rise with a speed of 5. In order to verify
+ that the wave is advecting at the correct speed, this may be changed to
+ -7, which should then show the wave to be stationary.
+
+ An alternative initial porosity distribution of a vertically changing
+ noisy 1D solitary wave may be set from the file 'sWaveSetup.xml' by
+ modifying the 'referenceSolutionFileName' parameter to
+ './input/solitaryWaves1DGlobal.dat'. This should then show a set of 2D
+ solitary porosity waves emerging from the 1D distribution with time.
+
+ A 3D model may also be run by changing directories to::
+
+ cd Magma/Models/Milestone4/SolitaryWaves3D
+
+ and then repeating the procedures detailed above for the 2D system. The
+ initial condition, read in from the file
+ './input/solitaryWaves3DGlobal.dat', is that of a single 1D solitary wave
+ in the vertical direction set against a noisy background distribution
+ which evolves with time into a set of 3D solitary waves.
+
+
+Related item(s)
+
+ Emerging 2D Solitary Waves
+ 2D Solitary Wave
+ Noisy 1D Solitary Wave Initial Condition
+ Emergent 2D Solitary Waves
+ Emergent 3D Solitary Waves
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/milestone4results/
Added: doc/geodynamics.org/benchmarks/trunk/magma/milestone5.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/milestone5.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/milestone5.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,66 @@
+Plone Metadata
+
+ milestone5results
+
+ Milestone 5 - Results and Analysis
+
+ Results and analysis for the isoviscous McKenzie equations (with melting)
+ driven by a corner flow velocity BC.
+
+
+Problem Description
+
+ This model describes the full isoviscous McKenzie equations (with
+ melting), for a system driven by a corner flow velocity boundary
+ condition for the solid phase. These equations couple a Stokes system for
+ the solid phase with a Darcy flow for the melt moving through the
+ permeable solid, and an advection term for the porosity field. The flow
+ is driven by a corner flow boundary condition for the solid, which
+ creates a region of low dynamic pressure about the area of discontinuity,
+ and a linear ramp in the melting function.
+
+Running the Simulation
+
+ The model is run from the directory::
+
+ Magma/Models/Milestone4/IsoviscousMcKenzieRidge2D/
+
+ and executing as::
+
+ ./StGermain IsoviscousMcKenzieRidge2D.xml
+
+ taking care to create a soft link to the StGermain binary in the build
+ directory as before.
+
+3D Model with Ridge Velocity Boundary Conditions
+
+ A 3D model was also implemented, which is driven by Direchlet BCs on the
+ velocity field, which are interpolated onto the prescribed domain from an
+ input file (the same one as for Milestone 1). The directory and execution
+ command for running this model are given as::
+
+ Magma/Models/Milestone5/IsoviscousMcKenzieRidge3D
+ ./StGermain IsoviscousMcKenzieRidge3D.xml
+
+ Some images for the x- and y- velocity components, the dynamic and
+ compaction pressures, the porosity, the melt fraction, and the x- and y-
+ melt velocity components are attached below.
+
+
+Related Item(s)
+
+ Isoviscous McKenzie System with Corner Flow BC - 1
+ Isoviscous McKenzie System with Corner Flow BC - 50
+ Isoviscous McKenzie System with Corner Flow BC - 3200
+ Velocity - x component
+ Velocity - y component
+ Dynamic (Stokes) pressure
+ Porosity
+ Compaction pressure
+ Melt fraction
+ Melt velocity - x component
+ Melt velocity - y component
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/milestone5results/
Added: doc/geodynamics.org/benchmarks/trunk/magma/running-stgmadds.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/magma/running-stgmadds.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/magma/running-stgmadds.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,38 @@
+Plone Metadata
+
+ stgmadds
+
+ Running stgMADDs Benchmarks
+
+ The Magma Development team has finished the alpha release of the Magma
+ Dynamics Demonstration Suite (MADDs). The initial code implements the
+ zero porosity/no melting magma benchmark for mid-ocean ridge solid flows
+ in 2D and 3D built on the Underworld framework. The purpose of this code
+ is principally to validate accurate pressure solvers for Stokes flow in
+ current CIG supported software. The stgMADDs source code is available in
+ CIG's Mercurial Repository (geodynamics.org/hg).
+
+Download and Install stgMADDs
+
+ For a first time download of the stgMADDs repository, do the following:
+
+ 1 Create the topmost repository with::
+
+ hg clone http://geodynamics.org/hg/magma/3D/stgMADDs
+
+ 2 Then obtain all the other repositories using::
+
+ ./obtainRepositories.py
+
+ 4 To push, you may have to use the ssh syntax, e.g.::
+
+ hg push ssh://hg@geodynamics.org/hg/magma/3D/stgMADDs
+
+ Caveat Emptor: This is very much an alpha release code for
+ experimentation with the accuracy of different mixed FEM pressure
+ solvers. Questions, complaints and bug reports should be directed to
+ "cig-magma at geodynamics.org":mailto:cig-magma at geodynamics.org
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/magma/workarea/benchmark/stgmadds/
Added: doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/index.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/index.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,104 @@
+
+General Description of the Benchmark Problem
+============================================
+
+ The benchmark will be similar to the benchmark of Blankenbach *et al.*
+ (1989) in methodology.
+
+ The benchmark problem is 2-D thermal convection of a non-rotating
+ anelastic liquid of infinite Prandtl number in a Cartesian, closed, unit
+ cell. The governing equations is based on Truncated An-Elastic
+ Approximation (TALA).
+
+ (attach equations as image here...)
+
+ We will have several cases, steady or unsteady, constant or variable
+ viscosity, bottom or internal heated, heat or mechanically driven.
+
+
+Grids
+-----
+
+ Each case will be run at **3 different resolutions** (grid resolution
+ 32x32, 64x64, 128x128, or higher if needed) to quantify the convergence
+ asymptotically. By comparing the asymptotically converged result, we
+ probably can negate the need of mesh refinement near the boundary and
+ reduce the uncertainty associated with various
+ interpolation/extrapolation schemes in calculating derived information
+ (e.g. geoid).
+
+Velocity BC's
+-------------
+
+ All boundaries (top, bottom, left, right) are **impermeable** (i.e., zero
+ normal velocity) and **free-slip** (i.e., zero tangential stress), except
+ for the mechanically driven case, where the top boundary is impermeble
+ and zero-slip (i.e. fixed horizontal velocity).
+
+Temperature BC's
+----------------
+
+ All non-dimensional numbers are defined at the top surface. There are
+ five non-dimensional numbers:
+
+ * **Ra**: Rayleigh number
+
+ * **H**: volumentric heat production number, **H = 0**, except for
+ internal heated cases.
+
+ * **Di**: Dissipation number
+
+ * **Gamma**: Gruneisen parameter
+
+ * **T_0**: Surface temperature, **T_0 = 0.1** for all cases.
+
+Reference State
+---------------
+
+ The reference density profile is $rho_ref(z) = exp((1-z)*Di/Gamma)$
+
+ The reference temperature profile is $T_ref(z) = T_0 * exp((1-z) * Di) T_0$
+
+ These physical properties are constant: thermal diffusivity, coefficient
+ of thermal expansion, gravitational acceleration.
+
+
+Required Information (all quantities are non-dimensional unless specified)
+--------------------------------------------------------------------------
+
+ * Nusselt number
+
+ * Mean Temperature
+
+ * Total Kinetic Energy
+
+ * RMS(V_x at top surface)
+
+ * Max(V_x at top surface)
+
+ * Total Dissipation Heating
+
+ * Total Adiabatic Cooling
+
+ * Dynamic Topography: Values required at 4 corners.
+
+ * Geoid Anomaly (dimensional): values required at the upper corners.
+ The following dimensional constants (in SI units) are used for the
+ calculation of geoid:
+
+ * Gravitational constant **G** = 6.673x10^-11
+
+ * depth of the box **R** = 3x10^6
+
+ * density at surface **rho_0** = 4000
+
+ * coefficient of thermal expansion **alpha_0** = 3x10^-5
+
+ * temperature contrast **Delta_T** = 3000
+
+ * viscosity **visc_0** = 10^21
+
+ * thermal diffusivity **kappa_0** = 10^-6
+
+ The density is 0 above the top of the box and is 2 below the bottom of
+ the box.
Added: doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite1.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite1.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite1.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,50 @@
+
+Suite 1
+=======
+
+ Testing the implementation of driving forces, isoviscous, comparing with
+ analytical solutions.
+
+ The non-dimensional numbers are moderately low in this case (Ra=10^5,
+ Di=0.25, gamma=1.5). The viscosity is constant. The purpose of this suite
+ is to ensure the driving forces are implemented correctly. This set of
+ non-dimensional numbers will give low compressibility and slow
+ convection. Therefore, most of the codes should behave well in this
+ parameter range.
+
+Case 1a: Bottom Heated
+----------------------
+
+ The temperature at the bottom is fixed at 1. The initial temperature
+ condition is
+
+ T = 0.5 everwhere, except at z = 0.5,
+ where T = 0.5 + cos(pi * x) * 0.001 * elz
+
+ where elz is the number of elements in the z-direction.
+
+ The initial temperature perturbation mimics a delta function. However,
+ the amplitude of the perturbation, with a factor of 1/1000, doesn't match
+ with a delta function. A comparison with analytical solution is possible
+ for the 0th-step velocity.
+
+Case 1b: Internal Heated
+------------------------
+
+ The temperature BC at the bottom is no-heatflux. The initial temperature
+ is the same as Case 1a. H = 1.
+
+Case 1c: Mechanically Driven
+----------------------------
+
+ The temperature at the bottom is fixed at 0. The initial temperature is
+ zero everywhere. The horizontal velocity BC at the top boundary is
+
+ V_x = 1000 * x^2 * (x-1)^2
+
+ so that V_x = 0 at x = 0 and 1, and dV_x/dx = 0 at x = 0 and 1.
+ (This case is optional.)
+
+ image:: Vx.png
+ "Horizontal velocity boundary condition"
+
Added: doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite2.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite2.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite2.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,30 @@
+
+=======
+Suite 2
+=======
+
+
+Steady state, basal heated
+==========================
+
+ Gamma = 1.1 in this suite.
+
+ Temperature dependent viscosity: eta = exp(-5 * T)
+
+Case 2a
+-------
+
+ * Ra = 3x10^5, Di = 0.5
+
+
+Case 2b
+-------
+
+ * Ra = 10^6, Di = 0.75 (not sure whether this case can reach steady state)
+
+Case 2c
+-------
+
+ * Ra = 3x10^6, Di = 1.0 (not sure whether this case can reach steady state)
+
+
Added: doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite3.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite3.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite3.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,12 @@
+
+=======
+Suite 3
+=======
+
+Steady state, internal heated
+=============================
+
+ This suite will have several cases taken from Jarvis and McKenzie (1980)
+
+ H = 1, no-heatflux for the bottom boundary.
+
Added: doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite4.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite4.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/2d-cartesian/suite4.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,10 @@
+
+=======
+Suite 4
+=======
+
+Time-dependent, unstead convection
+==================================
+
+ Parameters to be determined ...
+
Added: doc/geodynamics.org/benchmarks/trunk/mc/3d-spherical/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/3d-spherical/index.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/3d-spherical/index.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,285 @@
+
+=========================================
+3D Spherical Mantle Convection Benchmarks
+=========================================
+
+Presents initial benchmark results for CitcomS and comparisons with Stemmer et
+al {2006} and Ratcliff et al {1996} and Yoshida and Kageyama {2004} when the
+solutions are available from these codes. The community's comments and
+suggestions are solicited (use "add comment" feature at the bottom of the
+page).
+
+(1) Introduction
+=================
+
+There have been quite a lot of new developments in 3D spherical mantle
+convection models in the last few years, particularly in Japan and Europe.
+It seems to be a good time now to start doing some benchmarks for these
+codes. Benchmarks for spherical mantle convection were discussed at the
+Boulder mantle convection workshop in June of 2005. After the workshop,
+there were some email discussions on defining benchmark problems. This
+year, Uli Hansen's group published a paper in PEPI in which they included
+some benchmark calculations from their code and several other codes
+{Stemmer et al., 2006, PEPI}. The papers by Stemmer et al {2006}, Ratcliff
+et al. {1996, JGR}, Zhong et al. {2000, JGR}, and Yoshida and Kageyama
+{2004, GRL} and the previous discussions form a basis for benchmark
+calculations presented here and for discussions of our future benchmark
+studies.
+
+Here we present some initial benchmark results for CitcomS (20 cases so
+far) and comparisons with Stemmer et al {2006} and Ratcliff et al {1996}
+and Yoshida and Kageyama {2004} when the solutions are available (9 of
+them) from these codes. We define benchmark cases and present averaged
+global quantities as well as local properties for benchmark cases. At the
+moment, all the cases are for entirely basal heating and steady state or
+weakly time-dependent results, and viscosity variations from
+temperature-dependence vary from 1 (mobile-lid regime) to 1e6 (stagnant lid
+regime). And Ra<sub>0.5</sub> defined by viscosity at T=0.5 varies from 7e3
+to 1e5 (i.e., to achieve steady state or weakly time-dependent state).
+
+We hope that the community will join our effort and discuss what other
+benchmark cases should be included in our future benchmarks (e.g., internal
+heating, time-dependent cases at relatively high Ra, ...).
+
+(2) Definition of Benchmark Cases
+=================================
+
+(2.1) General features
+----------------------
+
+All the cases assume the Boussinesq approximation and constant material
+properties except for viscosity. There are no phase transitions, no
+compositional anomalies, and no internal heating (i.e., entirely basal
+heating) for all the cases. The top and bottom boundaries have radii
+r<sub>t</sub>=1 and r<sub>b</sub>=0.55, respectively.
+
+(2.2) Viscosity, Rayleigh number and scalings.
+----------------------------------------------
+
+The viscosity can be temperature-dependent and the non-dimensional
+viscosity, eta, is given by
+
+ image:: visc.gif (1)
+
+where E is the activation energy and T is non-dimensional temperature.
+Rayleigh number is defined as
+
+ image:: Ra.gif (2)
+
+where d is the thickness of the mantle (i.e., the length scale) and <font
+face="Symbol">h</font><sub>0.5</sub> is the viscosity at T=0.5. That is, Ra
+is what we sometimes call Ra<sub>0.5</sub>. This also implies that time and
+velocity scalings are defined using the thickness of the mantle, which is
+somewhat different from those in CitcomS which uses planetary radius as the
+length scale. All the CitcomS results in this benchmark study are re-scaled
+to use the mantle thickness as length scaling, to be consistent with those
+used in other studies.
+
+(2.3) Boundary conditions
+-------------------------
+
+Boundary conditions are free-slip at the top and bottom boundaries and
+isothermal with non-dimensional temperatures of 0 and 1 at the top and
+bottom boundaries, respectively.
+
+(2.4) Initial conditions.
+-------------------------
+
+Initial conditions come in two types. While most cases use initial
+temperature that is given as a function of coordinates (i.e., perturbation
+at some given spherical harmonics superimposed on a conductive temperature
+profile), some other cases use temperature field from other calculations as
+initial temperature. For the first type of initial condition, the initial
+temperature is given by
+
+ image:: IC.gif (3)
+
+where <font face="Symbol">e</font><sub>c</sub> and <font
+face="Symbol">e</font><sub>s</sub> are the magnitudes of perturbation for
+cosine and sinine terms respectively, l and m are spherical harmonic degree
+and order, r, <font face="Symbol">q</font>, and <font
+face="Symbol">f</font> are radial, co-latitude and longitude coordinates,
+and p<sub>lm</sub> is a normalized associated Legendre polynomial that is
+related to the associated Legendre polynomial P<sub>lm</sub> as:
+
+ image:: plm.gif (4)
+
+The magnitude of the perturbations (i.e., <font
+face="Symbol">e</font><sub>c</sub> and <font
+face="Symbol">e</font><sub>s</sub>) are 0.01, but for some cases, we only
+use the cosine term with <font face="Symbol">e</font><sub>s</sub>=0.0. For
+some cases, perturbations can be given at two different sets of harmonics
+(e.g., l=4, m=0 and l=4, m=4 for all the cubic symmetry cases). Detailed
+initial conditions for each case will be discussed later.
+
+We do not use random perturbations for benchmark studies here.
+
+(3) Quantifying Outputs of Benchmarks
+=====================================
+
+All the cases are computed to a steady state for global properties. Only
+steady state results are quantified. For some cases, we also give an output
+file that lists time-dependence of the results for better comparisons with
+other codes.
+
+(3.1) Global properties
+-----------------------
+
+For every 20 timesteps, we compute Nusselt numbers for both the top and
+bottom boundaries, averaged temperature for the whole mantle or shell <font
+face="Symbol">á</font>T<font face="Symbol">ñ</font>, and averaged
+RMS velocity <font face="Symbol">á</font>v<sub>rms</sub><font
+face="Symbol">ñ</font>.
+
+ image:: T.gif (5)
+
+ image:: vrms.gif (6)
+
+When a case reaches a steady state, we then compute the time-averaged
+values and standard deviation of Nusselt numbers, averaged temperature, and
+averaged RMS velocity over a certain period of time.
+
+(3.2) Local properties
+----------------------
+
+For every 20 timesteps, we also compute the maximum and minimum radial
+velocity and temperature at the middle depth of the mantle (i.e., at
+r=0.775). Again, we compute their time-averaged values and standard
+deviations over a certain period of time after a steady state is reached.
+
+These global and local properties are also used in Stemmer et al {2006} in
+their benchmark study.
+
+(4) Results
+===========
+
+Two sequences of cases are presented here, one associated with tetrahedral
+symmetry and the other with cubic symmetry. Ra<sub>0.5</sub> ranges from
+7e3 to 1e5 and viscosity variations due to temperature-dependence ranges
+from 1 (i.e., isoviscous) to 1e6 in stagnant lid regime. These cases are
+similar to what Ratcliff et al {1996} and Stemmer et al. {2006} had
+presented. However, we added more cases with larger viscosity contrasts.
+For cases with Ra<sub>0.5</sub>=7e3, we use 12x32x32x32 elements, and for
+cases with Ra<sub>0.5</sub>=1e5, 12x48x48x48 elements are used. Redial
+resolution is refined near the top and bottom boundary layers.
+
+In addition to those outputs in Tables, we also provide files that list the
+full time-dependence of our standard outputs. For some cases, thermal
+structure images are also provided.
+
+(4.1) The cases with tetrahedral symmetry and its variations
+------------------------------------------------------------
+
+For all the cases in this category, we use Ra<sub>0.5</sub>=7e3, but these
+cases have different temperature dependent viscosity, <font
+face="Symbol">Dh</font> of 1, 10, 20, 100, 1000, 1e4, 1e5, and 1e6 (Cases
+BM1A-BM1H in <a href=table1>Table 1</a> for definitions of all cases. The
+initial condition is the same for these cases with <font
+face="Symbol">e</font><sub>c</sub>=<font
+face="Symbol">e</font><sub>s</sub>=0.01 in equation 3. The first three
+cases, BM1A, 1B, and 1C were also presented in Stemmer et al {2006},
+Ratcliff et al. {1996}, and Yoshida and Kageyama {2004}; and Zhong et al.
+{2000} and Bercovici et al. {1989} computed case BM1A. Only the first five
+cases display tetrahedral symmetry. The last two cases, BM1G and 1H are in
+stagnant-lid regime, showing a large number of small plumes. Case BM1F with
+<font face="Symbol">Dh</font> of 1e4 is a transitional case.
+
+<a href=table2>Table 2</a> shows the output results and comparisons with
+previous studies when available. Note that in addition to average values,
+Table 2 also gives the time duration over which the averages are computed
+and standard deviations. It is clear that for Cases BM1A, 1B and 1C,
+CitcomS' results compare well with Stemmer et al. {2006}, Ratcliff et al.
+{1996} and Yoshida and Kageyama {2004}.
+
+Files for time-dependence of the outputs from t=0 to steady state and
+representative snapshots of residual temperature field for selective cases
+can be downloaded at <a
+href="http://anquetil.colorado.edu/szhong/CitcomS.html">Thermal Convection
+Benchmarks for CitcomS</a>. The format of these output files is also
+explained there.
+
+(4.2) The cases with cubic symmetry and its variations
+------------------------------------------------------
+
+For cubic symmetry cases with Ra<sub>0.5</sub>=7e3, we use the initial
+condition similar to Ratcliff et al. {1996}:
+
+ image:: IC40.gif (7)
+
+where sin(m <font face="Symbol">f</font>) is excluded. We have computed
+cases with temperature dependent viscosity, <font face="Symbol">Dh</font>
+of 1, 20, 30, 100, 1000, 1e4, 1e5 and 1e6 (Cases BM2A-BM2H in <a
+href=table1>Table 1</a> for definitions of all cases). The first three
+cases, BM2A, 2B, and 2C, can again be compared with previous studies. The
+first six cases (BM2A-BM2F) display cubic symmetry with 6 upwelling plumes.
+Different from tetrahedral cases, now BM2F with <font
+face="Symbol">Dh</font> of 1e4 also shows flow cubic symmetry, the same as
+cases with smaller viscosity contrast. The last two cases BM2G and BM2H are
+again in stagnant lid regime, with flow pattern breaking away from the
+cubic symmetry. Notice that the slight difference in Nu between BM1 and BM2
+cases for the same Ra.
+
+<a href=table3>Table 3</a> shows the output results and comparisons with
+previous studies when available. Again, CitcomS compares well with all the
+previous studies for BM2A, BM2B and BM2C.
+
+For cases with Ra<sub>0.5</sub>=1e5, we have computed cases with
+temperature dependent viscosity, <font face="Symbol">Dh</font> of 1, 10,
+30, and 100, so far (cases BM3A-3D, in <a href=table1>Table 1</a>). We use
+the same initial condition as in equation 7 only for the case with
+uniform viscosity that leads to cubic symmetry flow. For other
+Ra<sub>0.5</sub>=1e5 cases with variable viscosity, we use the final
+temperature field from either the uniform viscosity case (i.e., BM3A)
+or other cases as initial condition to help preserve the cubic symmetry
+flow. <a href=table1>Table 1</a> lists the detailed initial conditions.
+We found that if we use equation 7 as initial conditions for these high
+Ra cases with variable viscosity, we can no long achieve cubic symmetry
+flow (at least after 25000 time steps for BM3B). More plumes develop
+from the CMB than expected from cubic symmetry and the plumes are
+stable for a long time.
+
+<a href=table3>Table 3</a> shows the results and comparisons with Ratcliff
+et al {1996} for cases BM3A, BM3B and BM3C for which Ratcliff et al {1996}
+also published RMS velocity and Nu results. While the results from CitcomS
+compares reasonably with Ratcliff et al {1996}, the deviations are
+significantly larger than those for small Ra cases presented earlier. We
+think that the difference is caused by the resolution. Ratcliff et al.
+{1996} used the same resolution (i.e., 32, 64, and 128 cells in radial,
+co-latitude and longitude directions, respectively, for cubic symmetry
+cases) for different Ra cases. We think that cases with
+Ra<sub>0.5</sub>=1e5 should require higher resolution.
+
+Again, files for time-dependence of the outputs from t=0 to steady state
+and representative snapshots of residual temperature field for selective
+cases can be downloaded at <a
+href="http://anquetil.colorado.edu/szhong/CitcomS.html">Thermal Convection
+Benchmarks for CitcomS</a>.
+
+(5) Conclusions
+===============
+
+Here we present results of Nussult numbers, RMS velocity, averaged
+temperature, and maximum and minimum flow velocity and temperature at the
+mid-mantle depth for 20 cases from CitcomS. For these cases,
+Ra<sub>0.5</sub> is either 7e3 or 1e5, and viscosity contrast varies from 1
+to 1e6. The style of convection varies from mobile lid to stagnant lid
+regimes. For nine of the 20 cases, we could compare with previously
+published results {Stemmer et al., 2006; Ratcliff et al., 1996; Yoshida and
+Kageyama, 2004}. Comparisons show that CitcomS's results are generally
+consistent with these previous studies, although at high Ra, the relatively
+low resolution from previous studies may degrade the agreement. Time
+dependent results from CitcomS are compiled into a file for each case, and
+all of them can be downloaded for more detailed comparisons. Compared with
+benchmarks in Zhong et al. {2000}, here we added significantly more thermal
+convection cases.
+
+We encourage our colleagues to send in their results for these cases. We
+think that more benchmark cases are needed to test time-dependent
+convection at even higher Ra, and internal heating convection. We welcome
+suggestions.
+
+Acknowledgement: Mr. Nan Zhang, a graduate student at the University of
+Colorado, helped with calculations and analyses.
+
+
+
Added: doc/geodynamics.org/benchmarks/trunk/mc/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/index.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/index.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,41 @@
+
+
+============================
+Mantle Convection Benchmarks
+============================
+
+Benchmarks
+==========
+
+2D Cartesian Compressible Convection Benchmarks
+-----------------------------------------------
+
+ * Suite 1
+
+ * Suite 2
+
+ * Suite 3
+
+ * Suite 4
+
+
+3D Spherical Mantle Convection Benchmarks
+-----------------------------------------
+
+ * BM1{A-H}
+
+ * BM2{A-H}
+
+ * Bm3{A-D}
+
+
+
+Links
+-----
+
+ * "Mantle Convection Benchmarks":http://geodynamics.org/cig/workinggroups/mc/workarea/benchmark/
+
+ * "Benchmarks for 2D Cartesian compressible convection":http://geodynamics.org/cig/Members/tan2/benchmarks/
+
+ * "3D Spherical Mantle Convection Benchmarks":http://geodynamics.org/cig/workinggroups/mc/workarea/benchmark/3dconvention/
+
Added: doc/geodynamics.org/benchmarks/trunk/mc/notes-on-mantle-convection-benchmarks.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/mc/notes-on-mantle-convection-benchmarks.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/mc/notes-on-mantle-convection-benchmarks.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,78 @@
+
+Notes On Mantle Convection Benchmarks
+=====================================
+
+ (1) Ra = 3x10**6 (I think this is high enough to give some real time
+ dependence without pushing available resolution very much).
+
+ (2) Constaint properties (thermal expansivity, thermal diffusivity,
+ density, gravity, viscosity, internal heat generation) to keep things
+ very simple.
+
+ (3) Free slip upper and lower boundaries.
+
+ (4) Radius ratio = 0.546 (cmb/surface radius)
+
+ (5) purely internally heated
+
+ (6) insulating at cmb, constant temperature at surface
+
+ (7) model resolution: 65 nodes (64 layers) radially, with some packing of
+ nodes near the top and bottom boundaries. (We'll send you the actual
+ radii we use, assuming you can vary them at will.)
+
+ (8) initial diagnostics: (basically, these are just to get started and
+ see if we're in the same universe)
+
+ * (a) Nu vs. time (this should square with the internal heating in a
+ time-average sense)
+
+ * (b) Radial temperature profile vs. time - this is effectively a
+ measure of the efficiency of heat transfer, or equivalent of Nu for
+ bottom heated cases.
+
+ * (c) Spherical harmonic expansion of temperature field at all radial
+ levels at beginning and ending time (see below).
+
+ * (d) peak velocity and peak temperature in each radial layer vs. time
+
+ * (e) for now, let's ignore dynamic topography, since it's derived from
+ primitive results
+
+ (9) Initial conditions and run time: This is a bit thorny, so here's a
+ proposal. We can run TERRA to equilibrium under the specified model
+ conditions. Equilibrium is where Nu has settled down to fluctuations
+ about a steady mean value. At some point, call it time = 0.0, we'll stop
+ the code and output the full temperature field in the form of a spherical
+ harmonic expansion up to degree 128, which corresponds to the highest
+ model resolution. We can then restart both TERRA and CitcomS using this
+ spherical harmonic expansion (NOT the full temperature field at each
+ node, since this would prejudice things with regard to the particular
+ horizontal discretization.) Then both codes can run for a defined amount
+ of model time, keeping track of Nu, peak T, and peak V as a function of
+ time as indicated above. At the end of this time, or at several times
+ along the way, we can output spherical harmonic representations of T at
+ each layer for comparison.
+
+ I added the following comments:
+
+ (1) We use some analytic expressions for initial conditions (e.g., some
+ radial profile superimposed with a small perturbation of a given harmonic
+ function). In this way, others, if they want to benchmark their codes, do
+ not need to get the Terra output. Also in case some summary report comes
+ out of this effort, we can simply write down the initial conditions.
+
+ (2) We aim to reproduce four benchmark cases in steady of just one. The
+ four cases at the moment in my mind can be: three constant property cases
+ with purely basal heating at Ra=1e5 (case 1), and Ra=1e6 (case 2), and
+ purely internal heating at Ra=1e6 (case 3), and one temperature dependent
+ viscosity and purely basal heating at Ra=1e6 (case 4).
+
+ Case 1 will likely reach to a steady state, which is always a good thing
+ for a benchmark. Cases 2 and 3 are almost identical to what you have
+ suggested recently, and they are most likely time-dependent. The 1e6 Ra
+ is smaller than what you suggested today but is consistent with your
+ earlier suggestion. With Ra=1e6, we may not need grid refinement, which
+ is also good for benchmark purposes (again, others can do it later).
+ Case 4 is obviously of interest too.
+
Added: doc/geodynamics.org/benchmarks/trunk/seismo/index.rst
===================================================================
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-landers/description-landers.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-landers/description-landers.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-landers/description-landers.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,77 @@
+Plone Metadata
+ description-landers
+ Benchmark Description
+ Benchmark problem description
+
+Summary
+
+ Viscoelastic (Maxwell) relaxation of stresses from the 1992 M7.3 Landers earthquake,
+ focusing on the deformation in the area of the 1999 M7.1 Hector Mine earthquake.
+
+Problem Specification
+
+ Model size -- [NEED SPECS FOR CARL'S MESH]
+
+ Material properties
+
+ Elastic -- The material properties are a simplified 1-D version of the
+ 3-D SCEC Community Velocity Model. The 1-D model contains 11 layers with
+ uniform material properties within each layer and a minimum layer thickness
+ of 2 km. The elastic properties are given in an "ASCII":materials_layers2km.txt
+ file.
+
+ Viscoelastic -- Maxwell linear viscoelasticity (based on values in Pollitz, EPSL, 2003)
+
+ Upper crust (-19 km ≤ z) -- η = 1.0e+25 Pa-s (essentially elastic)
+
+ Lower crust (-30 km ≤ z < - 19 km) -- η = 32.2e+18 Pa-s
+
+ Mantle (z < -30 km) -- η = 4.6e+18 Pa-s
+
+ Fault geometry and slip distribution
+
+ The Landers and Hector Mine fault geometries and slip distribution
+ for Landers are incorporated into the LaGriT mesh.
+
+ Boundary conditions
+
+ Bottom and side displacements are pinned. Top of the model is a free surface.
+
+ Discretization
+
+ [GET SPECS FROM CARL'S MESH]
+
+ Element types
+
+ Linear and/or quadratic tetrahedral elements
+
+Requested Output
+
+ Solution
+
+ Displacements at all nodes at times of 0, 0.5, 1, 2, 4, and 7 years
+ as well as the mesh topology (i.e., element connectivity arrays and
+ coordinates of vertices) and basis functions. Also compute the traction
+ vector computed at the quadrature points of the faces making up the
+ Hector Mine faults.
+
+ June 30, 2006 -- Use ASCII output for now. In the future we will
+ switch to using HDF5 files.
+
+ Performance
+
+ * CPU time
+
+ * Wallclock time
+
+ * Memory usage
+
+ * Compiler and platform info
+
+"Truth"
+
+ You can't handle the truth
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-landers/description-landers
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/description-rs.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/description-rs.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/description-rs.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,104 @@
+Plone Metadata
+ description-rs
+ Benchmark Description
+ Benchmark problem description. Formerly known as benchmark 6b.
+
+Summary
+
+ Viscoelastic (Maxwell) relaxation of stresses from a single, finite, reverse-slip
+ earthquake in 3D with gravity. Evaluate results with imposed displacement boundary
+ conditions on a cube with sides of length 24 km. The displacements imposed are
+ the analytic elastic solutions. Symmetry boundary conditions are imposed at y = 0,
+ so the solution is equivalent to that for a domain with a 48 km length in the
+ y direction.
+
+ The effects of gravitational loading should be relaxed before the fault slip is
+ imposed. Alternatively, Winkler nodes could be used to calculate the gravitational
+ restoring forces resulting from the deformed upper surface.
+
+Problem Specificaqtion
+
+ Model size -- 0 km ≤ x ≤ 24 km; 0 km ≤ y ≤ 24 km; -24 km ≤ z ≤ 0 km
+
+ Top layer -- -12 km ≤ z ≤ 0 km
+
+ Bottom layer -- -24 km ≤ z ≤ -12 km
+
+ Material properties -- The top layer is nearly elastic whereas the bottom layer
+ is viscoelastic.
+
+ Elastic -- Poisson solid, G = 30 GPa, ρ = 3000 kg/m^3; g = 9.80665 m/s^2
+
+ Maxwell viscoelastic material properties
+
+ Top layer -- η = 1.0e+25 Pa-s (essentially elastic)
+
+ Bottom layer -- η = 1.0e+18 Pa-s
+
+ Boundary conditions
+
+ Bottom and side displacements set to analytic solution. (Note: the side
+ at y = 0 km has zero y-displacements because of the symmetry.) Top of the
+ model is a free surface.
+
+ Discretization
+
+ The model should be discretized with a nominal spatial resolution of 1000m,
+ 500m, and 250m. If possible, also run the models with a nominal spatial
+ resolution of 125 m. Optionally, use meshes with variable (optimal)
+ spatial resolution with the same number of nodes as the uniform resolution
+ meshes.
+
+ Element types
+
+ Linear and/or quadratic and tetrahedral and/or hexahedral
+
+ Fault specifications
+
+ Type -- 45 degree dipping reverse fault.
+
+ Location -- Strike parallel to y-direction with top edge at x = 4 km
+ and bottom edge at x = -12 km. 0 km ≤ y ≤ 16 km; -16 km ≤ z ≤ 0 km
+
+ Slip distribution -- 1 m of uniform thrust slip motion for 0 km ≤ y ≤ 12 km
+ and -12 km ≤ z ≤ 0 km with a linear taper to 0 slip at y = 16 km and z = -16 km.
+ In the region where the two tapers overlap, each slip value is the minimum
+ of the two tapers (so that the taper remains linear).
+
+ Boundary conditions
+
+ Lateral and bottom displacements are set to analytic elastic solution.
+ Note that the side at y = 0 km has zero y-displacements because of the
+ imposed symmetry at y = 0 km.
+
+Requested Output
+
+ Solution
+
+ Displacements at all nodes at times of 0, 1, 5, and 10 years
+ as well as the mesh topology (i.e., element connectivity arrays and
+ coordinates of vertices) and basis functions.
+
+ June 30, 2006 -- Use ASCII output for now. In the future we will switch
+ to using HDF5 files.
+
+ Performance
+
+ * CPU time
+
+ * Wallclock time
+
+ * Memory usage
+
+ * Compiler and platform info
+
+"Truth"
+
+ Okada routines are available to generate an elastic solution. The 'best'
+ viscoelastic answer will be derived via mesh refinement. Analytical
+ solutions to the viscoelastic problem are being sought if anyone has
+ any information.
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-rs/description-rs
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/results/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/results/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs/results/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,11 @@
+Plone Metadata
+ results
+
+ Results
+
+ Results from benchmark runs. Place tarballs containing the requested results
+ in this folder and describe the run in the `description` field.
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-rs/results
+
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/description-rs-nog.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/description-rs-nog.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/description-rs-nog.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,93 @@
+Plone Metadata
+ description-rs-nog
+ Benchmark Description
+ Benchmark problem description. Formerly known as benchmark 5b.
+
+Summary
+
+ Viscoelastic (Maxwell) relaxation of stresses from a single finite, reverse-slip
+earthquake in 3-D without gravity. Evaluate results with imposed displacement boundary
+conditions on a cube with sides of length 24 km. The displacements imposed are the
+analytic elastic solutions. Symmetry boundary conditions are imposed at y = 0,
+so the solution is equivalent to that for a domain with a 48 km length in the
+y direction.
+
+Problem Specification
+
+ Model size -- 0 km ≤ x ≤ 24 km; 0 km ≤ y ≤ 24 km; -24 km ≤ z ≤ 0 km
+
+ Top layer -- -12 km ≤ z ≤ 0 km
+
+ Bottom layer -- -24 km ≤ z ≤ -12 km
+
+ Material properties -- The top layer is nearly elastic whereas the bottom layer is viscoelastic.
+
+ Elastic -- Poisson solid, G = 30 GPa
+
+ Viscoelasticity -- Maxwell linear viscoelasticity
+
+ Top layer -- η = 1.0e+25 Pa-s (essentially elastic)
+
+ Bottom layer -- η = 1.0e+18 Pa-s
+
+ Fault specifications
+
+ Type -- 45 degree dipping reverse fault.
+
+ Location -- Strike parallel to y-direction with top edge at x = 4 km,
+ and bottom edge at x = 20 km. 0 km ≤ y ≤ 16 km; -16 km ≤ z ≤ 0 km
+
+ Slip distribution -- 1 m of uniform thrust slip motion for 0 km ≤ y ≤ 12 km
+ and -12 km ≤ z ≤ 0 km with a linear taper to 0 slip at y = 16 km and
+ z = -16 km. In the region where the two tapers overlap, each slip value
+ is the minimum of the two tapers (so that the taper remains linear).
+
+ Boundary conditions
+
+ Bottom and side displacements set to analytic solution. (Note: the side
+ at y = 0 km has zero y-displacements because of symmetry). Top of the
+ model is a free surface.
+
+ Discretization
+
+ The model should be discretized with nominal spatial resolutions of
+ 1000 m, 500 m, 250 m. If possible, also run the models with a nominal
+ spatial resolution of 125 m. Optionally, use meshes with variable (optimal)
+ spatial resolution with the same number of nodes as the uniform resolution
+ meshes.
+
+ Element types
+
+ Linear and/or quadratic and tetrahedral and/or hexahedral.
+
+
+Requested Output
+
+ Solution
+
+ Displacements at all nodes at times of 0, 1, 5, and 10 years as well as
+ the mesh topology (i.e., element connectivity arrays and coordinates of
+ vertices) and basis functions.
+
+ June 30, 2006 -- Use ASCII output for now. In the future we will switch
+ to using HDF5 files.
+
+ Performance
+
+ * CPU time
+
+ * Wallclock time
+
+ * Memory usage
+
+ * Compiler and platform info
+
+"Truth"
+
+ Okada routines are available to generate an elastic solution. The 'best'
+ viscoelastic answer will be derived via mesh refinement. Analytical solutions
+ to the viscoelastic solution are being sought if anyone has information.
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-rs-nog/description-rs-nog
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/geofest-input/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/geofest-input/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/geofest-input/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,23 @@
+GeoFEST Input
+
+ Input files for GeoFEST
+
+ * bmrsnog_tet4_1000m.gft.gz (2006-08-31)
+ Gzipped GeoFEST input file for benchmark using linear tetrahedral elements
+ with a 1000m nominal node spacing.
+
+ * bmrsnog_tet4_0500m.gft.gz (2006-08-31)
+ Gzipped GeoFEST input file for benchmark using linear tetrahedral elements
+ with a 500m nominal node spacing.
+
+ * reverse slip (no grav), refined grid 01, no smoothing (GeoFEST 4.5)
+ (2006-09-06) Carl Gable's mesh,
+ see http://meshing.lanl.gov/proj/crustal_dyn_reverse_fault_bm/catalog.html
+
+ * reverse slip (no grav), refined grid 02, no smoothing (GeoFEST 4.5)
+ (2006-09-06) Carl Gable's mesh #02,
+ see http://meshing.lanl.gov/proj/crustal_dyn_reverse_fault_bm/catalog.html
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-rs-nog/geofest-input
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/plots/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/plots/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/plots/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,45 @@
+Plone Metadata
+ Plots of Reverse-Slip No Gravity Benchmark Results
+ Plots of global and local errors for reverse-slip no gravity benchmark
+
+Displacement Field
+
+ "PyLith soln":img:tet4_1000m_pylith_disp_t00.png
+
+ "GeoFEST soln":img:tet4_1000m_geofest_disp_t00.png
+
+Global Error
+
+ "Plot of global error":img:globalerror.png
+
+Local Error
+
+ Elastic solution: Code versus Analytic
+
+ 1000m resolution
+
+ "PyLith error":img:tet4_1000m_pylith_analytic_t00.png
+
+ "GeoFEST error":img:tet4_1000m_geofest_analytic_t00.png
+
+ "COMSOL error":img:tet10_2000m_femlab_analytic_t00.png
+
+ 500m resolution
+
+ "PyLith error":img:tet4_0500m_pylith_analytic_t00.png
+
+ "GeoFEST error":img:tet4_0500m_geofest_analytic_t00.png
+
+ Viscoelastic solution: PyLith versus GeoFEST
+
+ "t0yr":img:tet4_0500m_pylith_geofest_t00.png
+
+ "t1yr":img:tet4_0500m_pylith_geofest_t01.png
+
+ "t5yr":img:tet4_0500m_pylith_geofest_t05.png
+
+ "t10yr":img:tet4_0500m_pylith_geofest_t10.png
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-rs-nog/plots
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/pylith-0.8-input/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/pylith-0.8-input/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/pylith-0.8-input/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,30 @@
+PyLith-0.8 Input
+
+ Input files for PyLith-0.8
+
+ * bmrsnog_hex_1000m.tgz (2006-07-20)
+ Tarball containing PyLith-0.8 input files for benchmark using linear hexahedral
+ elements with a 1000m nominal node spacing.
+
+ * bmrsnog_tet4_1000m.tgz (2006-07-20)
+ Tarball containing PyLith-0.8 input files for benchmark using linear tetrahedral
+ elements with a 1000m nominal node spacing.
+
+ * bmrsnog_tet4_0500m.tgz (2006-07-20)
+ Tarball containing PyLith-0.8 input files for benchmark using linear tetrahedral
+ elements with a 500m nominal node spacing.
+
+ * bmrsnog_tet4_0250m.tgz (2006-07-20)
+ Tarball containing PyLith-0.8 input files for benchmark using linear tetrahedral
+ elements with a 250m nominal node spacing.
+
+ * reverse slip (no grav), refined grid 01, no smoothing (2006-09-06)
+ Carl Gable's mesh,
+ see http://meshing.lanl.gov/proj/crustal_dyn_reverse_fault_bm/catalog.html
+
+ * reverse slip (no grav), refined grid 02, no smoothing (2006-09-06)
+ Carl Gable's mesh #2,
+ see http://meshing.lanl.gov/proj/crustal_dyn_reverse_fault_bm/catalog.html
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-rs-nog/pylith-0.8-input
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/results/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/results/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-rs-nog/results/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,69 @@
+Results
+
+ Results from benchmark runs. Place tarballs containing the requested results
+ in this folder and describe the run in the description field.
+
+ * GeoFEST reverse fault results - 1 km (2006-08-17)
+ Tarball contains input and output files as well as text file
+ containing run-time information
+
+ * GeoFEST reverse fault results - 500 m (2006-08-17)
+ Tarball contains input and output files as well as text file
+ containing run-time information
+
+ * Geofest reverse slip var_res_mesh_01_soln (2006-09-05)
+ fixed the BCs, Geofest 4.5, dt=0.1 constant
+
+ * PyLith, 1 proc, linear hex, 1 km resolution, dt=0.1yr (2006-08-29)
+ PyLith results run on 1 processor of a Power Mac G5.
+ Linear hexahedral mesh at 1 km resolution.
+ Constant time step size of 0.1 years.
+
+ * PyLith, 1 proc, linear tet, 1 km resolution, dt=0.1yr (2006-08-29)
+ PyLith results run on 1 processor of a Power Mac G5.
+ Linear tetrahedral mesh at 1 km resolution.
+ Constant time step size of 0.1 years.
+
+ * PyLith, 1 proc, linear tet, 500 m resolution, dt=0.1yr (2006-08-29)
+ PyLith results run on 1 processor of a Power Mac G5.
+ Linear tetrahedral mesh at 500 m resolution.
+ Constant time step size of 0.1 years.
+
+ * PyLith, 1 proc, linear tet, 250 m resolution, dt=0.1yr (2006-08-29)
+ PyLith results run on 1 processor of a Power Mac G5.
+ Linear tetrahedral mesh at 250 m resolution.
+ Constant time step size of 0.1 years.
+
+ * tet_var_res_01_pylith_soln.tgz (2006-09-04)
+ PyLith-0.8 results for Carl Gable's variable resolution (no smoothing)
+ mesh 01 for the reverse slip benchmark - constant dt=0.1yr
+
+ * Femlab 2 km resolution, elastic (2006-10-16)
+ This model has 19544 quadratic tetrahedral elements and is twice the
+ size in y of the model description, since there is no symmetric boundary.
+ This yields a resolution close to 2 km. The model and solver require
+ about 800 MB and is solved in about 10 minutes on a 1.8 GHz AMD Opteron.
+
+ * Femlab 1 km resolution, t = 0 years (2006-10-18)
+ This model has ~162000 linear tetrahedral elements and is twice the size
+ in y of the model description, since there is no symmetric boundary.
+ This yields a resolution close to 1 km. The model and the solver require
+ about 800MB and is solved in about 3 minutes on a 1.8 GHz AMD Opteron.
+ An iterative solver was used, which uses the Incomplete LU preconditioner
+ with a drop tolerance of 0.01
+
+ * Femlab 1 km resolution, t = 1 year
+ Viscoelastic problem requires ~3.5GB and takes about 4.5 hrs to run.
+ Drop tolerance is 0.01
+
+ * Femlab 1 km resolution, t = 5 years
+ Viscoelastic problem requires ~3.5GB and takes about 4.5 hrs to run.
+ Drop tolerance is 0.01
+
+ * Femlab 1 km resolution, t = 10 years
+ Viscoelastic problem requires ~3.5GB and takes about 4.5 hrs to run.
+ Drop tolerance is 0.01
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-rs-nog/results
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/description-ss.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/description-ss.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/description-ss.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,103 @@
+Plone Metadata
+
+ description-ss
+ Benchmark Description
+ Benchmark problem description. Formerly known as benchmark 4b.
+
+Summary
+
+ Viscoelastic (Maxwell) relaxation of stresses from a single, finite,
+strike-slip earthquake in 3-D without gravity. Evaluate results with imposed
+displacement boundary conditions on a cube with sides of length 24 km. The
+displacements imposed are the analytic elastic solutions. Anti-plane strain
+boundary conditions are imposed at y = 0, so the solution is equivalent
+to that for a domain with a 48 km length in the y direction.
+
+Problem Specification
+
+ "Problem geometry":img:benchmark_geometry.png
+
+ Model size -- 0 km ≤ x ≤ 24 km; 0 km ≤ y ≤ 24 km; -24 ≤ z ≤ 0 km
+
+ Top layer -- -12 km ≤ z ≤ 0 km
+
+ Bottom layer -- -24 km ≤ z ≤ -12 km
+
+ Material properties -- The top layer is nearly elastic whereas the bottom layer
+ is viscoelastic.
+
+ Elastic -- Poisson solid, G = 30 GPa
+
+ Viscoelasticity -- Maxwell linear viscoelasticity
+
+ Top layer -- η = 1.0e+25 Pa-s (essentially elastic)
+
+ Bottom layer -- η = 1.0e+18 Pa-s
+
+ Fault specifications
+
+ Type -- Vertical right-lateral strike-slip fault.
+
+ Location --
+ Strike parallel to y-direction at center of model (x = 12 km)
+ 0 km ≤ y ≤ 16 km; -16 km ≤ z ≤ 0 km.
+
+ Slip distribution --
+ 1 m of uniform strike slip motion for 0 km ≤ y ≤ 12 km
+ and -12 km ≤ z ≤ 0 km with a linear taper to 0 slip
+ at y = 16 km and z = -16 km. In the region where the two
+ tapers overlap, each slip value is the minimum of the
+ two tapers (so that the taper remains linear).
+
+ Boundary conditions
+
+ Bottom and side displacements are set to the elastic analytical solution,
+ and the top of the model is a free surface. There are two exceptions to
+ these applied boundary conditions. The first is on the y = 0 plane, where
+ y-displacements are left free to preserve symmetry, and the x- and
+ z-displacements are set to zero. The second is along the line segment
+ between (12, 0, -24) and (12, 24, -24), where the analytical solution
+ blows up in some cases. Along this line segment, all 3 displacement
+ components are left free.
+
+ Discretization
+
+ The model should be discretized with nominal spatial resolutions of
+ 1000 m, 500 m, and 250 m. If possible, also run the models with a nomial
+ spatial resolution of 125 m. Optionally, use meshes with variable
+ (optimal) spatial resolution with the same number of nodes as the
+ uniform resolution meshes.
+
+ Element types
+
+ Linear and/or quadratic and tetrahedral and/or hexahedral.
+
+
+Requested Output
+
+ Solution
+
+ Displacement at all nodes at times of 0, 1, 5, and 10 years as well
+ as the mesh topology (i.e., element connectivity arrays and coordinates
+ of vertices) and basis functions.
+
+ June 30, 2006 -- Use ASCII output for now. In the future we will switch
+ to using HDF5 files.
+
+ Performance
+
+ * CPU time
+
+ * Wallclock time
+
+ * Memory usage
+
+ * Compiler and platform info
+
+"Truth"
+
+ Okada routines are available to generate an elastic solution. The 'best'
+ viscoelastic answer will be derived via mesh refinement. Analytical solutions
+ to the viscoelastic solution are being sought if anyone has information.
+
+
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/geofest-input/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/geofest-input/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/geofest-input/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,28 @@
+
+GeoFEST Input
+
+ Input files for GeoFEST
+
+ * GeoFEST linear tet 1km resolution dt = 0.1 year
+
+ * GeoFEST linear tet 500m resolution dt = 0.1 year
+
+ * GeoFEST linear tet 250m resolution input file
+
+ * GeoFEST / PYRAMID 1km
+ PYRAMID input file for parallel 1km run.
+
+ * GeoFEST / PYRAMID 500m
+ PYRAMID input file for parallel GeoFEST run.
+
+ * GeoFEST / PYRAMID 250m
+ PYRAMID input file for parallel GeoFEST run.
+
+ * GeoFEST linear tet 500m dt = 0.1 year (NEW)
+ The taper problem has been fixed
+
+ * GeoFEST linear tet 250m dt = 0.1 year (NEW)
+ The taper problem has been fixed.
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-strikeslip/geofest-input
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,29 @@
+
+Strike-Slip Benchmark (no gravity)
+
+ Benchmark for strike-slip fault without gravity.
+
+ * Benchmark Description
+ Benchmark problem description. Formerly known as benchmark 4b.
+
+ * PyLith-0.8 Input
+ Input files for PyLith-0.8
+
+ * GeoFEST Input
+ Input files for GeoFEST
+
+ * Results
+ Results from benchmark runs. Place tarballs containing the
+ requested results in this folder and describe the run in the
+ description field.
+
+ * Plots of Benchmarking Results
+ Plots of benchmarking results showing global and local errors.
+
+ * Geometry for Strike-Slip Benchmark
+ Domain and fault geometry for the strike-slip benchmark.
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks
+
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/plots/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/plots/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/plots/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,48 @@
+Plots of Strike-Slip No Gravity Benchmark Results
+ Plots of global and local errors for strike-slip no gravity benchmark
+
+Displacement Field
+
+ "PyLith soln":img:tet4_1000m_pylith_disp_t00.png
+
+ "GeoFEST soln":img:tet4_1000m_geofest_disp_t00.png
+
+Global Error
+
+ "Plot of global error":img:globalerror.png
+
+Local Error
+
+ Elastic solution: Code versus Analytic
+
+ 250m resolution
+
+ "PyLith error":img:tet4_0250m_pylith_analytic_t00.png
+
+ "GeoFEST error":img:tet4_0250m_geofest_analytic_t00.png
+
+ 500m resolution
+
+ "PyLith error":img:tet4_0500m_pylith_analytic_t00.png
+
+ "GeoFEST error":img:tet4_0500m_geofest_analytic_t00.png
+
+ Viscoelastic solution: PyLith versus GeoFEST
+
+ 250m resolution
+
+ "t0yr":img:tet4_0250m_pylith_geofest_t00.png
+
+ 500m resolution
+
+ "t0yr":img:tet4_0500m_pylith_geofest_t00.png
+
+ "t1yr":img:tet4_0500m_pylith_geofest_t01.png
+
+ "t5yr":img:tet4_0500m_pylith_geofest_t05.png
+
+ "t10yr":img:tet4_0500m_pylith_geofest_t10.png
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-strikeslip/plots
+
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/pylith-0.8-input/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/pylith-0.8-input/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/pylith-0.8-input/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,19 @@
+PyLith-0.8 Input
+
+ Input files for PyLith-0.8
+
+ * bmssnog_tet4_1000m.tgz
+ Tarball containing PyLith-0.8 input files for benchmark using
+ linear tetrahedral elements with a 1000m nominal node spacing.
+
+ * bmssnog_tet4_0500m.tgz
+ Tarball containing PyLith-0.8 input files for benchmark using
+ linear tetrahedral elements with a 500m nominal node spacing.
+
+ * bmssnog_tet4_0250m.tgz
+ Tarball containing PyLith-0.8 input files for benchmark using
+ linear tetrahedral elements with a 250m nominal node spacing.
+
+Original URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-strikeslip/pylith-0.8-input/
+
Added: doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/results/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/results/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/benchmark-strikeslip/results/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,84 @@
+Results
+
+ Results from benchmark runs. Place tarballs containing the requested results
+ in this folder, and describe the run in the description field.
+
+ * PyLith, 1 proc, linear tet, 1km resolution, dt=0.1yr (2006-08-29)
+ PyLith results run on 1 processor of a Power Mac G5.
+ Linear tetrahedral mesh at 1km resolution.
+ Constant time step size of 0.1 year.
+
+ * PyLith, 1 proc, linear hex, 1km resolution, dt=0.1yr (2006-08-29)
+ PyLith results on 1 processor of a Power Mac G5.
+ Linear hexahedral mesh at 1km resolution.
+ Constant time step size of 0.1 year.
+
+ * PyLith, 1 proc, linear tet, 500m resolution, dt=0.1yr (2006-08-29)
+ PyLith results run on 1 processor of a Power Mac G5.
+ Linear tetrahedral mesh at 500m resolution.
+ Constant time step size of 0.1 year.
+
+ * PyLith Revised Results, 500m, New BC and Split Node Input (2007-01-30)
+ New solution using revised BC and split node inputs. The revised BC
+ take care of the problems of defining BC on the fault plane (or in
+ some cases the projected fault plane). The new split node inputs
+ no longer assume a bilinear slip distribution in the region where
+ the fault tapers overlap, and now assumes a taper consistent with
+ what is used for the analytical solution.
+
+ * PyLith Revised Results, 500m, Altered BC for Viscoelastic Solution (2007-02-06)
+ New version where BC have been altered from those of previous version
+ to make viscoelastic results consistent with those from GeoFEST.
+ The revised BC do not pin y-component on the y=0 plane, and no BC
+ are applied along the intersection of the fault plane (or its projection)
+ along y=0 and z=-24.
+
+ * PyLith, 1 proc, linear tet, 250m resolution, dt=0.1yr (2006-09-07)
+ PyLith results run on 1 processor of an Opteron 2.4GHz Linux machine.
+ Linear tetrahedral mesh at 250m resolution.
+ Constant time step size of 0.1 year.
+
+ * GeoFEST / PYRAMID 1km (2006-09-06)
+ Parallel results using 64 processors of Intel/Linux Cluster
+ with GeoFEST-4.5 and Pyramid-2.1.3
+
+ * GeoFEST / PYRAMID 500m (2006-09-06)
+ Parallel results using 64 processors of Intel/Linux Cluster
+ with GeoFEST-4.5 and Pyramid-2.1.3
+
+ * GeoFEST / PYRAMID 250m (2006-09-06)
+ Parallel results using 128 processors of Intel/Linux Cluster
+ with GeoFEST-4.5 and Pyramid-2.1.3
+
+ * GeoFEST linear tet 1km resolution dt=0.1yr (updated) (2006-09-21)
+ The taper error has been fixed.
+
+ * GeoFEST Linear-Tet 500m Re-Run (2006-11-29)
+
+ * GeoFEST Linear-Tet 250m Re-Run (2006-11-29)
+
+ * Femlab 1km resolution, t = 0 years (2006-10-17)
+ This model has ~162,000 linear tetrahedral elements and is twice the
+ size in y of the model description, since there is no symmetric boundary.
+ This yields a resolution close to 1km. The model and solver require
+ almost 800MB and is solved in about 3 minutes on a 1.8 GHz AMD Opteron.
+ An iterative solver was used, which uses the Incomplete LU preconditioner
+ with a drop tolerance of 0.01. Decreasing this value has very little
+ effect on the error but takes longer to solve.
+
+ * Femlab 1km resolution, t = 1 year (2006-10-17)
+ Viscoelastic problem requires ~3.5GB and takes about 4.5 hours to run.
+ Drop tolerance is 0.01.
+
+ * Femlab 1km resolution, t = 5 years (2006-10-17)
+ Viscoelastic problem requires ~3.5GB and takes about 4.5 hours to run.
+ Drop tolerance is 0.01.
+
+ * Femlab 1km resolution, t = 10 years (2006-10-17)
+ Viscoelastic problem requires ~3.5GB and takes about 4.5 hours to run.
+ Drop tolerance is 0.01.
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-strikeslip/results
+
Added: doc/geodynamics.org/benchmarks/trunk/short/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,52 @@
+
+Short-Term Tectonics Benchmarks
+
+"Overview":overview
+
+Benchmarks
+
+ Strike-slip (no gravity)
+
+ * "Description":benchmark-strikeslip/description-ss
+
+ * "PyLith input":benchmark-strikeslip/pylith-0.8-input
+
+ * "GeoFEST input":benchmark-strikeslip/geofest-input
+
+ * "Submitted results":benchmark-strikeslip/result
+
+ * "Plots of results":benchmark-strikeslip/plots
+
+ Reverse-slip (no gravity)
+
+ * "Description":benchmark-rs-nog/description-rs-nog
+
+ * "PyLith input":benchmark-rs-nog/pylith-0.8-input
+
+ * "GeoFEST input":benchmark-rs-nog/geofest-input
+
+ * "Submitted results":benchmark-rs-nog/results
+
+ * "Plots of results":benchmark-rs-nog/plots
+
+ Reverse-slip (with gravity)
+
+ * "Description":benchmark-rs/description-rs
+
+ * PyLith input (coming soon)
+
+ * GeoFEST input (coming soon)
+
+ * "Submitted results":benchmark-rs/results
+
+ Landers-Hector Mine
+
+ * "Description":benchmark-landers/description-landers
+
+ * Mesh constructed with LaGriT (coming soon)
+
+Utilities
+
+ * "Analytic and Semi-Analytic Codes":utilities
+
+ * "CUBIT examples":utilities/CUBITex
Added: doc/geodynamics.org/benchmarks/trunk/short/overview.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/overview.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/overview.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,69 @@
+Short Name: overview
+Title: Overview
+Description: Overview of benchmark suite
+
+Overview of Goals and Objectives
+
+ In order to test the accuracy and speed of various elastic and viscoelastic
+finite element calculations using different codes on different platforms, we
+have developed the following benchmark comparisons. Information resulting from
+the benchmark comparisons will be used for the following purposes.
+
+ 1 Confirming proper numerical implementation of the physics include
+ rheological laws, fault constitutive laws, etc.
+
+ 1 Testing the accuracy of the numerical implementations as a function
+ of meshing scheme, number of nodes, element type, time-stepping scheme,
+ code, etc.
+
+ 1 Testing the computational efficiency of different codes, solvers, and
+ modeling techniques as a function of meshing scheme, number of nodes,
+ element type, time-stepping scheme, code, etc.
+
+ 1 Comparing and evaluating available finite element codes.
+
+ Based on the comparisons, we would like to be able to (1) identify and correct any
+errors in numerical implementation which currently exist in any of the considered
+codes, (2) quantify differences in numerical solutions as a function of meshing
+scheme, number of nodes, element type, time-stepping scheme, code, etc., and (3)
+quantify and, if possible, minimize model induced uncertainties resulting from
+discretization, model boundaries, unexpected transients in time-dependent materials,
+etc.
+
+General Methodology
+
+ All benchmark descriptions assume a right-handed Cartesian coordinate system
+with the x-direction running east, the y-direction running north, and the
+z-direction running up. If a boundary condition is applied at a depth, d, this
+will correspond to z = -d. The surface is always assumed to be at z = 0.
+Use whatever coordinate system is most convienient for your program, but please
+convert the results to the one defined here.
+
+ Benchmark meshes will be described at the coarsest level to be run. Because meshes
+may be either structured or unstructured, the mesh nodal spacing described refers
+to the average. If memory, time, and patience allow, also run models at 1/2, 1/4,
+1/8, etc. the original coarse node spacing. This will make it possible to see how
+accuracy and speed scale with mesh spacing. If your code permits a variety of
+element types, also run models using various types of elements (linear vs.
+quadrilateral; hexahedral vs. tetrahedral, full vs. reduced integration). This
+will make it possible to see how accuracy and speed change with element type.
+Finally, variable mesh spacing degrades accuracy, but, for economy, we would like
+to employ a variable mesh (e.g., to resolve stress variations at the fault tips,
+etc.) If time permits, investigate the trade-offs involved in using variable
+resolution meshes.
+
+ With regard to output, there are a number of parameters which should be noted
+for each model. For the purposes of determining accuracy, please record displacements
+at all nodes and stresses (all 6 independent components) at all Gauss points
+at each specified time. For the purposes of evaluating performance, please try
+to keep track of memory usage (including the size of stiffness matrix and mean
+execution time, etc.) More details regarding the submission of your results for
+inclusion in the summary analysis can be found at this document. To ease the
+burden on those who are compiling the data, your results will not be accepted
+unless they are in the specified format.
+
+ When available, analytical solutions for the various benchmarks will be given
+on the CIG website at http://XXXXXXXX. Whenever possible, please check to make
+sure your results are essentially correct before submitting them for the summary
+analysis.
+
Added: doc/geodynamics.org/benchmarks/trunk/short/savage-prescott/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/savage-prescott/index.rst (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/savage-prescott/index.rst 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,107 @@
+
+=============================
+Savage and Prescott Benchmark
+=============================
+
+ This benchmark problem computes the viscoelastic (Maxwell) relaxation of
+ stresses from repeated infinite, strike-slip earthquakes in 3D without
+ gravity. The files needed to run the benchmark may be found at
+ geodynamics.org/svn/cig/short/2.5D/benchmarks/savageprescott. An analytical
+ solution to this problem is described by Savage and Prescott
+ [Savage:Prescott:1978], which provides a simple way to check our numerical
+ solution. A python utility code is provided in the utils directory to compute
+ the analytical solution. Although this problem is actually 2.5D (infinite
+ along-strike), we solve it using a 3D finite element model.
+
+
+Problem Description
+===================
+
+ Figure [fig:benchmark:savageprescott:geometry] shows the geometry of the
+ problem, as described by [Savage:Prescott:1978]. The analytical solution
+ describes the surface deformation due to repeated earthquakes on an infinite
+ strike-slip fault embedded in an elastic layer overlying a Maxwell
+ viscoelastic half-space. The upper portion of the fault (red in the figure)
+ is locked between earthquakes, while the lower portion (blue in the figure)
+ creeps at plate velocity. At regular recurrence intervals, the upper portion
+ of the fault abruptly slips by an amount equal to the plate velocity
+ multiplied by the recurrence interval, thus 'catching up' with the lower part
+ of the fault.
+
+ There are some differences between the analytical solution and our numerical
+ representation. First, the analytical solution represents the earthquake
+ cycle as the superposition of uniform fault creep and an elementary
+ earthquake cycle. Uniform fault creep is simply the uniform movement of the
+ two plates past each other at plate velocity. For the elementary earthquake
+ cycle, no slip occurs below the locked portion of the fault (blue portion in
+ the figure). On the locked (red) portion of the fault, backslip equal to
+ plate velocity occurs until the earthquake recurrence interval, at which
+ point abrupt forward slip occurs. In the finite element solution, we perform
+ the simulation as described in the figure. Velocity boundary conditions are
+ applied at the extreme edges of the model to simulate block motion, steady
+ creep is applied along the blue portion of the fault, and regular earthquakes
+ are applied along the upper portion of the fault. It takes several earthquake
+ cycles for the velocity boundary conditions to n as described in the figure.
+ Velocity boundary conditions are applied at the extreme edges of the model to
+ simulate block motion, steady creep is applied along the blue portion of the
+ fault, and regular earthquakes are applied along the upper portion of the
+ fault. It takes several earthquake cycles for the velocity boundary
+ conditions to n as described in the figure. Velocity boundary conditions are
+ applied at the extreme edges of the model to simulate block motion, steady
+ creep is applied along the blue portion of the fault, and r
+
+ Due to the difficulties in representing solutions in an infinite domain,
+ there are several meshes that have been tested for this problem. The simplest
+ meshes have uniform resolution (all cells have equal dimensions); however,
+ such meshes typically do not provide accurate solutions since the resolution
+ is too coarse in the region of interest. For that reason, we also tested
+ several meshes where the mesh resolution decreases away from the center. In
+ the problem description that follows, we will focus on the graded3 mesh (in
+ meshes/spbm_hex8_graded4_20km.exo.gz), which provides reasonable results. It
+ will first be necessary to unzip this mesh so that it may be used by PyLith.
+
+
+ Domain The domain for this mesh spans the region-1000\leq x\leq1000\ km,The top (elastic) layer occupies the region -40\ km\ \leq z\leq0 and the bottom (viscoelastic) layer occupies the region -400\ km\ \leq z\leq-40\ km.
+
+ Material properties The material is a Poisson solid with a shear modulus (\mu) of 30 GPa. The domain is modeled using an elastic isotropic material for the top layer and a Maxwell viscoelastic material for the bottom layer. The bottom layer has a viscosity (\eta) of 4.73364e19 Pa-s, yielding a relaxation time (2\eta/\mu) of 100 years.
+
+ Fault The fault is a vertical, left-lateral strike-slip fault. The strike is parallel to the y-direction at the center of the model:x=0\ km,The locked portion of the fault (red section in Figure [fig:benchmark:savageprescott:geometry]) extends from -20\: km\leq z\leq0, while the creeping section (blue) extends from -40\: km\leq z\leq0. Along the line where the two sections coincide (z=-20\: km), half of the coseismic displacement and half of the steady creep is applied (see parameters/finalslip_rupture.spatialdb and parameters/sliprate_creep.spatialdb).
+
+ Boundary conditions On the bottom boundary, vertical displacements are set to zero, while on the y-boundaries the x-displacements are set to zero. On the x-boundaries, the x-displacements are set to zero, while constant velocities of +/- 1 cm/yr are applied in the y-direction, giving a relative plate motion of 2 cm/year.
+
+ Discretization For the graded3 mesh (meshes/spbm_hex8_graded4_20km.exo), the resolution at the outer boundaries is 20 km. An inner region is then put through one level of refinement, so that near the center of the mesh the resolution is 10 km. All meshes were generated with Cubit.
+
+ Basis functions We use trilinear hexahedral cells.
+
+ Solution We compute the surface displacements and compare these to the
+ analytical solution [fig:benchmark:savageprescott:solution].
+
+Running the Benchmark
+=====================
+
+After checking out the benchmark files from the CIG SVN repository, change to the meshes directory. Decompress the gzipped files in the mesh directory,
+
+ gunzip \*.gz
+
+Alternatively, simply unzip the mesh you want to use (e.g., meshes/spbm_hex8_graded3_20km.exo). Change to the parameters directory. There are a number of .cfg files corresponding to the different meshes, as well as a pylithapp.cfg file defining parameters common to all problems, and a timedep.cfg file defining fault slip behavior and output frequency. Each problem uses three .cfg files: pylithapp.cfg, a mesh-specific file (e.g., spbm_hex8_graded3_20km.cfg), and timedep.cfg. There are also scripts for each mesh. You can then run the problem either as
+
+ pylith spbm_hex8_graded3_20km.cfg timedep.cfg
+
+or as
+
+ sh ./spbm_hex8_graded3_20km.sh
+
+This will run the problem for 10 earthquake cycles of 200 years each, using a time-step size of 10 years, for a total simulation time of 2000 years. Ground surface output occurs every 10 years, while all other outputs occur every 50 years. The graded3 benchmark requires approximately 5 GB of disk space.
+
+Once the problem has run, results will be placed in the results directory in the mesh-specific subdirectory. These results may be viewed directly with a package such as ParaView; however, to compare results to the analytical solution, some postprocessing is required. First, generate the analytical results by going to the utils directory and running the savpres_ss.py script. This will produce comma-delimited files for displacements and velocities (savpres_displ.csv, savpres_vel.csv) that are easy to use with a plotting package. There is an additional script in the utils directory (vtkdiff.py) that computes velocities from PyLith displacement results. To generate velocity results, you can type, for example,
+
+./vtkdiff.py vtkdiff_hex8_graded3_20km.cfg
+
+Results will be placed in the results directory for the appropriate model. This will create a set of VTK files that you can use with ParaView. You can generate profile results using ParaView (e.g., the `Plot Over Line` filter). By loading the entire time series for the ground surface velocity results, it is possible to generate profile results for all time steps. For comparison with the analytical results, it is easiest to use a profile coincident with the x-axis that extends from x=0 to the right-hand side of the mesh. The profiles can also be saved as CSV files from ParaView.
+
+
+Benchmark Results
+=================
+
+Figure [fig:benchmark:savageprescott:solution] shows the computed surface displacements for the 10th earthquake cycle compared with the analytical solution. The profile results were obtained as described above, and then all results (analytical and numerical) were referenced to the displacements immediately following the last earthquake. We find very good agreement between the analytical and numerical solutions, even for meshes with uniform refinement. This is a relatively recent benchmark, and we have not yet explored quantitative fits as a function of mesh resolution. For this benchmark, it is also important to consider the distance of the boundary from the region of interest. Also note that the agreement between analytical and numerical solutions is poor for early earthquake cycles, due to the differences in simulating the problem, as noted above.
+
Added: doc/geodynamics.org/benchmarks/trunk/short/utilities/CUBITex.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/utilities/CUBITex.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/utilities/CUBITex.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,17 @@
+
+Example Journal Files for CUBIT
+
+ From the 2009 NMCDEF meeting in Golden, CO
+
+ * Meshing example (2009-06-23)
+ Very short example of meshing a pyramid and displaying mesh
+
+ * Geometry test (2009-06-23)
+ Example of building geometrical shapes with merging, subtracting, moving, etc.
+
+ * Fault example (2009-06-23)
+ Example of building and meshing (coarsely) a region around a dipping fault patch.
+
+
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/utilities/CUBITex
Added: doc/geodynamics.org/benchmarks/trunk/short/utilities/index.txt
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/short/utilities/index.txt (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/short/utilities/index.txt 2009-09-14 23:45:35 UTC (rev 15673)
@@ -0,0 +1,2 @@
+URL
+ http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/utilities
More information about the CIG-COMMITS
mailing list