[cig-commits] r15984 - doc/geodynamics.org/benchmarks/trunk/cigma
luis at geodynamics.org
luis at geodynamics.org
Mon Nov 16 14:50:04 PST 2009
Author: luis
Date: 2009-11-16 14:50:04 -0800 (Mon, 16 Nov 2009)
New Revision: 15984
Modified:
doc/geodynamics.org/benchmarks/trunk/cigma/index.rst
Log:
Updated cigma/index.rst
Modified: doc/geodynamics.org/benchmarks/trunk/cigma/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/cigma/index.rst 2009-11-16 22:49:56 UTC (rev 15983)
+++ doc/geodynamics.org/benchmarks/trunk/cigma/index.rst 2009-11-16 22:50:04 UTC (rev 15984)
@@ -1,13 +1,30 @@
-In this chapter, we show how to use Cigma to run specific comparisons on
-included datasets, estimate the order of convergence of the numerical
-methods used in obtaining those solutions, and visualize the resulting
-error fields.
+.. |times| unicode:: U+00D7
+.. |L2| replace:: *L*\ :sub:`2`
+.. |10^5| replace:: 10\ :sup:`5`
+
+.. |8x8| replace:: 8×8
+.. |16x16| replace:: 16×16
+.. |32x32| replace:: 32×32
+.. |64x64| replace:: 64×64
+.. |512x512| replace:: 512×512
+.. |err64x64| replace:: [;||\\phi_n - \\phi_{64\\times 64}||_{L_2} / \\sqrt{V};]
+
+.. |8x8x8| replace:: 8×8×8
+.. |16x16x16| replace:: 16×16×16
+.. |32x32x32| replace:: 32×32×32
+.. |64x64x64| replace:: 64×64×64
+.. |err64x64x64| replace:: [;||\\phi_n - \\phi_{64\\times 64\\times 64}||_{L_2} / \\sqrt{V};]
+
+We can use Cigma to compare our datasets, and estimate the order of
+convergence of the numerical methods used in obtaining those solutions.
+We can also visualize the resulting error fields.
+
Poisson Problem
---------------
-Cigma can compare two functions and return the [;L_{2};]-norm difference
+Cigma can compare two functions and return the |L2|-norm difference
between those two functions. For a given numerical problem, if you
compute the solutions over a range of resolutions, you can use this
global measure to quantify the rate of convergence of your numerical
@@ -27,38 +44,64 @@
finite element software library called Deal.II, which supports Lagrange
finite elements of any order in either 2 or 3 dimensions. In fact, the
two equations we have described are already solved in Step 4 in the list
-of the Deal.II tutorial programs in dealii.org/developer/doxygen/tutorial/.
+of the Deal.II tutorial programs in
+http://dealii.org/developer/doxygen/tutorial/.
-In the ``examples/`` subdirectory, we include a version of the Step 4 tutorial
-program that has been slightly modified to suit our purposes in this section.
+In the ``examples/`` subdirectory, we include a version of the
+Step 4 tutorial program that has been slightly modified to suit
+our purposes in this section.
+
First, we change the output format to use VTK files, in order to make it very
convenient for us to provide our input datasets to Cigma. Next, we solve the
-2D problem over meshes of resolution [;64\times 64;], [;32\times 32;],
-[;16\times 16;], and [;8\times8;] by using linear quadrilateral elements.
+2D problem over meshes of resolution |64x64|, |32x32|, |16x16|, |8x8|
+by using linear quadrilateral elements.
The solution is saved in file ``phi_64x64x64.vtk`` for the highest resolution,
for example. In a similar fashion, we solve our 3D Poisson problem on meshes
-of resolution 64\times64\times64, 32\times32\times32, 16\times16\times16,
-and 8\times8\times8, that use linear hexahedral elements. Note that in
+of resolution |64x64x64|, |32x32x32|, |16x16x16|, |8x8x8|,
+that use linear hexahedral elements. Note that in
either case, we don't really know the exact solution to either problem,
so we use the highest resolution available as the reference point for
our comparisons. The following bash shell command can take care of all
-the comparisons::
+the comparisons
+::
+
$ for r in 32 16 8; do
cigma compare phi_64x64.vtk phi_${r}x${r}.vtk
cigma compare phi_64x64x64.vtk phi_${r}x${r}x${r}.vtk
done
When we run the above code on the command line, or script, Cigma will generate
-a summary of each comparison, including the [;L_{2};] norm of the difference between
+a summary of each comparison, including the |L2| norm of the difference between
the two functions, and other basic mesh information such as its volume V and
maximum cell diagonal h. We summarize that information in the following tables:
-.. fig:
-.. fig:
+.. table:1
++---------+----------+------------+
+| [;n;] | [;h_n;] | |err64x64| |
++---------+----------+------------+
+| |8x8| | 0.35355 | 0.012312 |
++---------+----------+------------+
+| |16x16| | 0.17677 | 0.003157 |
++---------+----------+------------+
+| |32x32| | 0.08838 | 0.000689 |
++---------+----------+------------+
+
+.. table:1
+
++------------+----------+---------------+
+| [;n;] | [;h_n;] | |err64x64x64| |
++------------+----------+---------------+
+| |8x8x8| | 0.43301 | 0.019205 |
++------------+----------+---------------+
+| |16x16x16| | 0.21650 | 0.004889 |
++------------+----------+---------------+
+| |32x32x32| | 0.10825 | 0.001075 |
++------------+----------+---------------+
+
It should become apparent from these results that the global error metric
[;\\varepsilon;] does indeed decrease with increasing resolution at the
expected rate. To verify the exact order of convergence, you can use
@@ -66,8 +109,13 @@
the following plots:
-.. fig:
+.. fig:1
+.. figure:: images/bm1_alpha.png
+ :align: center
+ :width: 800
+ _`Figure 1`: Convergence of |L2| Global Error for Poisson
+ problems in 2D and 3D.
Based on the slope of the lines in the above log-log plots, we have indeed
@@ -75,12 +123,13 @@
from a method that uses linear elements.
+
Mantle Convection
-----------------
In addition to reporting a global error metric, Cigma can also output
local errors between two fields, averaged over each of the integration
-cells used in the calculation of the [;L_{2};] norm (or error). You can
+cells used in the calculation of the |L2| norm (or error). You can
use this scalar error field to ascertain a number of physical insights
by correlating the spatial regions that contribute the most to the
global error.
@@ -88,7 +137,7 @@
For our next example, we use CitcomCU, a CIG code for mantle convection,
to solve a thermal convection problem inside a three-dimensional domain
under base heating, stress-free boundary conditions, constant viscosity,
-and using a Rayleigh number of 10^{5} in the domain \Omega=[0,1]^{3}.
+and using a Rayleigh number of |10^5| in the domain [;\\Omega=[0,1]^{3};].
Under these conditions, we expect the solution to this problem to eventually
converge to a steady state consisting of a single convection cell. We then
use Cigma to compare those steady state solutions in order to recover the
@@ -98,15 +147,26 @@
Again, as in the previous section, we solve the same problem over four
different resolutions, and then use the highest resolution dataset as
the reference point for all the subsequent comparisons. The four cases
-we will use for this example use 64\times64\times64, 32\times32\times32,
-16\times16\times16, and 8\times8\times8 elements. In the following figure,
+we will use for this example use |64x64x64|, |32x32x32|, |16x16x16|,
+and |8x8x8| elements. In the following figure,
we display the temperature and velocity fields on the y=0.99 plane, near
one of the faces of the domain. We can distinctly identify an upwelling
and downwelling cycle, as well as sharp gradients in the temperature
-field. The variations in the velocity field are smoother throughout the slice.
+field. The variations in the velocity field are smoother throughout
+the slice.
-
+.. fig:2
+.. figure:: images/bm2_out_fig_fields_64_20000_y099.png
+ :align: center
+ :width: 800
+ _`Figure 2`: Slice of the steady-state temperature and velocity
+ fields for our thermal convection problem, on the y=0.99 plane.
+ In the temperature plot (a), we show ten equally spaced contour
+ values. In the velocity plot (b), the vectors represent the xz-plane
+ components of the velocity, while the colors represent the
+ velocity magnitude.
+
Once we obtain the solutions for our four resolutions over 20000
solution steps, we process the ASCII output data from CitcomCU
into a number of rectilinear grids in VTK format (``*.vtr`` files),
@@ -115,11 +175,17 @@
As a preliminary step, we also need to ensure that each of our cases
do indeed reach a steady state. We can verify this fact numerically
by using Cigma to compare the same field at two different times, and
-checking that the L_{2}-error is negligible. However, for this problem,
+checking that the |L2|-error is negligible. However, for this problem,
we can also visually verify that the solution has reached steady state
by checking that the oscillations in the values of the average heat
-flux on the top surface have dampened out, as shown by Figure 5.3 (XXX),
+flux on the top surface have dampened out, as shown by `Figure 3`,
+.. fig:3
+.. figure:: images/bm2_heatflux.png
+ :align: center
+ :width: 800
+
+ _`Figure 3`: Time variation of average heat flux over top surface.
Our first comparison of two steady states is simple enough. Although
@@ -136,7 +202,7 @@
Since we did not specify an integration mesh through the ``-m`` option,
Cigma will take the mesh from the first field and use it for the integration
-of the [;L_{2};] error. We have, however, used the ``-o`` option to specify
+of the |L2| error. We have, however, used the ``-o`` option to specify
an output file. This will create the HDF5 file ``steady-state.h5`` after
the comparison is complete, and store the results of that comparison into
an array called ``error_temperature_64_32`` under the HDF5 / root group.
@@ -178,7 +244,7 @@
-b case8.2100.pvtr:velocity \
-o steady-state.h5:/error_velocity_64_8
-At this point, you have obtained the local L_{2}-errors for each
+At this point, you have obtained the local |L2|-errors for each
of the above comparisons, and have stored them in arrays inside the
file ``steady-state.h5``. In order to visualize these errors, you
will need to reattach the original mesh information by using a
@@ -200,34 +266,48 @@
The first two options we give our post-processing program are important.
The first option eliminates the volume scale factor in the definition
-of the local L_{2} error (larger cells contribute more than smaller cells,
+of the local |L2| error (larger cells contribute more than smaller cells,
which would be more obvious if the error were constant). The second option
tells the post-processing program to use a logarithmic scale for writing
those locally averaged errors. Using a log-10 scale makes it much easier
to spot the spatial range of variation in our error fields.
-In Figure 5.4, we show two plots representing the errors in the temperature
-field on the y=0.99 plane for the two 16\times16\times16 and 32\times32\times32
-resolutions, relative to the 64\times64\times64 temperature solution. Examining
-this plot and Figure 5.2a (XXX) reveals that the largest errors are concentrated
-in areas where the temperature gradient is also large, near the boundary layers
-created by activity from the upwelling and downwelling. Similarly, the region
-where the errors are smallest corresponds to a region where both the temperature
+In `Figure 4`, we show two plots representing the errors in the temperature
+field on the y=0.99 plane for the two |16x16x16| and |32x32x32|
+resolutions, relative to the |64x64x64| temperature solution.
+Examining this plot and Figure 5.2a (XXX) reveals that the
+largest errors are concentrated in areas where the temperature gradient
+is also large, near the boundary layers created by activity from
+the upwelling and downwelling. Similarly, the region where the errors
+are smallest corresponds to a region where both the temperature
gradients and velocity magnitude are at their lowest.
-.. fig:
+.. fig:4
+.. figure:: images/bm2_out_fig_log_error_temperature_64_16_64_32_y099.png
+ :align: center
+ :width: 800
+ _`Figure 4`: Temperature differences on y=0.99 plane, as compared
+ against case with |64x64x64| elements.
-Likewise, in Figure 5.5 we show two plots representing the magnitude of the
+
+Likewise, in `Figure 5` we show two plots representing the magnitude of the
velocity field on the same plane. This time, the distribution of the largest
errors are spread out more evenly over the domain. In this case, the largest
errors in the velocity field also correspond to regions where the magnitude
of the velocity field is the largest.
-.. fig:
+.. fig:5
+.. figure:: images/bm2_out_fig_log_error_velocity_64_16_64_32_y099.png
+ :align: center
+ :width: 800
+ _`Figure 5`: Velocity differences on y=0.99 plane, as compared
+ against case with |64x64x64| elements.
+
+
Circular Inclusion Benchmark
----------------------------
@@ -239,37 +319,109 @@
exact analytical solution is known. In a 2003 paper [Schmid Podladchikov 2003],
Schmid and Podladchikov derived an analytic solution for the pressure and
velocity fields of a circular inclusion under simple shear, depicted in
-Figure 5.6. The viscosity parameters introduced are \mu_{m} for the
-viscosity of the matrix, and \mu_{i} for the viscosity of the inclusion.
+Figure 5.6. The viscosity parameters introduced are [;\\mu_{m};] for the
+viscosity of the matrix, and [;\\mu_{i};] for the viscosity of the inclusion.
The kinematic boundary conditions are given generally in terms of the
-simple shear strain rate \dot{\epsilon}.
+simple shear strain rate [;\\dot{\\epsilon};].
+.. fig:6
+.. figure:: images/bm3_inclusion_setup.png
+ :align: center
+ _`Figure 6`: Schematic for the circular inclusion benchmark.
-In particular, we place the inclusion of radius r_{i}=0.1 at the origin,
+
+In particular, we place the inclusion of radius [;r_{i}=0.1;] at the origin,
and exploit the symmetry in this problem by only solving the field on the
top right quarter of the domain.
From [Schmid Podladchikov 2003], we end up with the following analytic
formula for the pressure field under the case of simple shear,
-p=\begin{cases}
-4\dot{\epsilon}\frac{\mu_{m}(\mu_{i}-\mu_{m})}{\mu_{i}+\mu_{m}}\left(\frac{r_{i}^{2}}{r^{2}}\right)\cos(2\theta) & r>r_{i}\\
-0 & r<r_{i}\end{cases}where we use \mu_{i}=2 for the viscosity of the inclusion,
-\mu_{m}=1 for the viscosity of the background media, and \dot{\epsilon}=1
-for the magnitude of the shear. In Cigma, you can refer to this specific
-analytic solution by the somewhat verbose name bm.circular_inclusion.pressure,
-whose definition is summarized in Listing 5.1.
+ [;
+ p=\\begin{cases}
+ 4\\dot{\\epsilon}
+ \\frac{\\mu_{m}(\\mu_{i}-\\mu_{m})}{\\mu_{i}+\\mu_{m}}
+ \\left(\\frac{r_{i}^{2}}{r^{2}}\\right)\\cos(2\\theta) & r>r_{i}\\\\
+ 0 & r<r_{i}\\end{cases}
+ ;]
+where we use [;\\mu_{i}=2;] for the viscosity of the inclusion,
+[;\\mu_{m}=1;] for the viscosity of the background media,
+and [;\\dot{\\epsilon}=1;] for the magnitude of the shear.
+In Cigma, you can refer to this specific
+analytic solution by the somewhat verbose name ``bm.circular_inclusion.pressure``,
+whose definition is summarized in `Listing 1`.
+.. listing:1
+
+::
+
+ #include “Function.h”
+
+ namespace benchmark {
+ namespace circular_inclusion {
+ class Pressure;
+ }
+ }
+
+ class benchmark::circular_inclusion::Pressure : public cigma::Function {
+ public:
+ Pressure() {}
+ ~Pressure() {}
+
+ virtual void init() { }
+ virtual FunctionType getType() const { return Function::GeneralType; }
+
+ virtual int n_dim() { return 2; }
+ virtual int n_rank() { return 1; }
+
+ virtual bool eval(double *x, double *y) {
+ const double R = 0.1;
+ const double mu_m = 1;
+ const double mu_i = 2;
+ const double mag_shear = 1.0;
+ const double C = 4*mag_shear*(mu_m*(mu_i-mu_m)/(mu_i+mu_m))*R*R;
+ const double xc = 0.0;
+ const double yc = 0.0;
+ const double dx = x[0] - xc;
+ const double dy = y[0] - yc;
+ const double r2 = dx*dx + dy*dy;
+
+ if (r2 < R*R) {
+ value[0] = 0.0;
+ } else {
+ double theta = atan2(dy, dx);
+ value[0] = (C/r2) * cos(2*theta);
+ }
+
+ return true;
+ }
+ };
+
+
+
One additional step must be taken before we can use the name
-bm.circular_inclusion.pressure on the Cigma command line. We must
-instantiate our new C++ class, and register it under an assigned
-name in the constructor defined in the FunctionRegistry.cpp source file.
+``bm.circular_inclusion.pressure`` on the Cigma command line.
+We must instantiate our new C++ class, and register it under an assigned
+name in the constructor defined in the ``FunctionRegistry.cpp`` source file.
+.. listing:2
+::
+ FunctionRegistry::FunctionRegistry()
+ {
+ //
+ // ... other definitions
+ //
+ typedef benchmark::circular_inclusion::Pressure PressureFn1;
+ shared_ptr<PressureFn1> pressure1(new PressureFn1());
+ this->addFunction(“bm.circular_inclusion.pressure”, pressure1);
+ }
+
+
Now that we have an analytic expression for the pressure, we would like
to obtain an approximate numerical solution to this circular inclusion problem.
For this purpose we can use Gale, a CIG code for long-term crustal dynamics.
@@ -278,13 +430,26 @@
use the far-field approximation of the velocity field in order to impose
the kinematic boundary conditions in Gale,
-v_{x}Since we expect our solution to vary as p\sim1/r^{2}, we target an
-error of 0.01\% by enlarging the domain under consideration to about 80
-times the radius of the inclusion. This results in the final domain
-\Omega=[0,8]^{2}.
+ [;
+ \\begin{align*}
+ v_x &= -\\dot{\\epsilon}y, \\\\
+ v_y &= \\dot{\\epsilon}x,
+ \\end{align*}
+ ;]
+Since we expect our solution to vary as [;p \\sim 1/r^{2};],
+we target an error of 0.01\% by enlarging the domain under
+consideration to about 80 times the radius of the inclusion.
+This results in the final domain [;\\Omega=[0,8]^{2};].
+.. fig:7
+.. figure:: images/bm3_out_fig_pressure_512.png
+ :align: center
+ :width: 800
+ _`Figure 7`: Pressure field for circular inclusion problem,
+ with resolution of |512x512| elements.
+
First, we'd like to see how well the Gale solutions converge to a common
answer by comparing each other against the highest resolution field available.
Since we only have three solutions, that leaves us with only two meaningful
@@ -306,48 +471,88 @@
whose results we summarize in the following Table.
+.. table:p512
-
from which we can estimate the order of convergence
-\alpha=\log(0.00158/0.00407)/\log(0.0442/0.0884)=1.36.
+ [;
+ \\alpha=\\log(0.00158/0.00407)/\\log(0.0442/0.0884)=1.36.
+ ;]
-Processing the local error fields stored in circ_inc.h5, we obtain the following plots,
+Processing the local error fields stored in ``circ_inc.h5``,
+we obtain the following plots,
+.. fig:8
+.. figure:: images/bm3_out_fig_full_log_error_pressure_512_128_512_256.png
+ :align: center
+ :width: 800
+.. fig:9
+.. figure:: images/bm3_out_fig_inc_log_error_pressure_512_128_512_256.png
+ :align: center
+ :width: 800
-Next, we compare each of the three pressure fields obtained with Gale against the analytical formula for the pressure given earlier. We accomplish that by invoking the following bash command line,
-$ for r in 512 256 128; do
+Next, we compare each of the three pressure fields obtained with Gale
+against the analytical formula for the pressure given earlier.
+We accomplish that by invoking the following bash command line,
- cigma compare \
+::
- -a ${r}_8/fields.00000.pvts:PressureField \
+ $ for r in 512 256 128; do
+ cigma compare \
+ -a ${r}_8/fields.00000.pvts:PressureField \
+ -b bm.circular_inclusion.pressure \
+ -o inclusion.h5:/errror_pressure_${r}
+ done
- -b bm.circular_inclusion.pressure \
+whose output we summarize in the following table,
- -o inclusion.h5:/errror_pressure_${r}
+.. table
- done
+from which we can obtain the order of convergence [;\\alpha=0.5;],
+via a power-law regression analysis. Such a low value of [;\alpha;]
+is expected since the pressure function we are solving has a
+discontinuity across the inclusion boundary.
-whose output we summarize in the following table,
+After processing the errors in inclusion.h5 with the utility
+program visualize-errors, we obtain the following plots of the
+pressure field differences,
+.. fig:10
+.. figure:: images/bm3_out_fig_inc_log_error_pressure_256_512.png
+ :align: center
+ :width: 800
-from which we can obtain the order of convergence \alpha=0.5, via a power-law regression analysis. Such a low value of \alpha is expected since the pressure function we are solving has a discontinuity across the inclusion boundary.
+.. fig:11
+.. figure:: images/bm3_out_fig_inc_log_error_pressure_256_512.png
+ :align: center
+ :width: 800
-After processing the errors in inclusion.h5 with the utility program visualize-errors, we obtain the following plots of the pressure field differences,
+The only appreciable change between these two plots happens in
+the region near the inclusion boundary, which becomes slightly
+more resolved spatially. Also, in contrast to the previous
+two sections, we do not detect much of a color shift between
+two plots, even though the |L2| global error is indeed decreasing.
+This is due to the fact that the error remains large near the
+inclusion with approximately the same maximum value, making
+it harder to see any variation in the local error field.
+In summary, we have observed that for this problem the
+|L2| global error does indeed decrease as we increase
+the mesh resolution, but it does not do so at the optimal rate.
+This indicates that the numerical scheme used to obtain these
+solutions is not completely accurate. Note that many numerical
+codes that solve Stokes flow, Gale included, assume that the
+pressure, velocity, and viscosity fields are continuous, an
+assumption that is violated for this particular problem.
-
-The only appreciable change between these two plots happens in the region near the inclusion boundary, which becomes slightly more resolved spatially. Also, in contrast to the previous two sections, we do not detect much of a color shift between two plots, even though the L_{2} global error is indeed decreasing. This is due to the fact that the error remains large near the inclusion with approximately the same maximum value, making it harder to see any variation in the local error field.
-
-In summary, we have observed that for this problem the L_{2} global error does indeed decrease as we increase the mesh resolution, but it does not do so at the optimal rate. This indicates that the numerical scheme used to obtain these solutions is not completely accurate. Note that many numerical codes that solve Stokes flow, Gale included, assume that the pressure, velocity, and viscosity fields are continuous, an assumption that is violated for this particular problem.
More information about the CIG-COMMITS
mailing list