[cig-commits] r15976 - doc/geodynamics.org/benchmarks/trunk/cs
luis at geodynamics.org
luis at geodynamics.org
Mon Nov 16 14:48:59 PST 2009
Author: luis
Date: 2009-11-16 14:48:59 -0800 (Mon, 16 Nov 2009)
New Revision: 15976
Modified:
doc/geodynamics.org/benchmarks/trunk/cs/index.rst
Log:
Added cs/index.rst
Modified: doc/geodynamics.org/benchmarks/trunk/cs/index.rst
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/cs/index.rst 2009-11-16 22:48:56 UTC (rev 15975)
+++ doc/geodynamics.org/benchmarks/trunk/cs/index.rst 2009-11-16 22:48:59 UTC (rev 15976)
@@ -0,0 +1,353 @@
+
+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.
+
+Poisson Problem
+---------------
+
+Cigma can compare two functions and return the [;L_{2};]-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
+method. To demonstrate how this works, we will use the two Poisson
+problems:
+
+(1) [;\\nabla^{2}\\phi(x,y)=4x^{4}+4y^{4};] inside [;\\Omega=[-1,1]^{2};],
+ subject to the Dirichlet boundary condition
+ [;\\phi(x_{0},y_{0})=x_{0}^{2}+y_{0}^{2};] on the boundary
+ [;\\partial\\Omega;], and
+
+(2) [;\\nabla^{2}\\phi(x,y,z)=4x^{4}+4y^{4}+4z^{4};] on [;\\Omega=[-1,1]^{3};]
+ subject to [;\\phi(x_{0},y_{0},z_{0})=x_{0}^{2}+y_{0}^{2}+z_{0}^{2};]
+ on the boundary [;\\partial\\Omega;].
+
+These differential equations can be easily solved numerically by using a
+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/.
+
+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.
+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
+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::
+
+ $ 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
+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:
+
+
+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
+regression analysis to fit a power law through the above points, and obtain
+the following plots:
+
+
+.. fig:
+
+
+
+Based on the slope of the lines in the above log-log plots, we have indeed
+verified that the global errors behave as [;\\varepsilon=Ch^{2};], as expected
+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
+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.
+
+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}.
+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
+order of convergence of the CitcomCU numerical code and to examine how
+the error behaves on a representative slice of the original domain.
+
+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 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.
+
+
+
+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),
+and a subsequent ``*.pvtr`` file for each solution step.
+
+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,
+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),
+
+
+
+Our first comparison of two steady states is simple enough. Although
+it is not critical for this example since we have already converged,
+in general time-dependent problems you must ensure you are comparing
+solutions at the same absolute time. Thus, even though it appears
+that the following command is comparing two different solution steps,
+they actually correspond to roughly the same physical time.::
+
+ $ cigma compare \
+ -a case64.20000.pvtr:temperature \
+ -b case32.9900.pvtr:temperature \
+ -o steady-state.h5:/error_temperature_64_32
+
+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
+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.
+You may also refer to the same HDF5 file in subsequent runs of Cigma.
+The target dataset will be appended to the file if it doesn't exist,
+and overwritten if it does. The commands for other two comparisons are,
+
+::
+
+ $ cigma compare \
+ -a case64.20000.pvtr:temperature \
+ -b case16.4700.pvtr:temperature \
+ -o steady-state.h5:/error_temperature_64_16
+
+ $ cigma compare \
+ -a case64.20000.pvtr:temperature \
+ -b case8.2100.pvtr:temperature \
+ -o steady-state.h5:/error_temperature_64_8
+
+Comparisons involving the velocity can be obtained in a similar manner.
+Note that even though we are comparing vectors here, the resulting
+errors will be scalars, since Cigma is integrating the scalar function
+[;||\\vec{v}_{a}(\\vec{x})-\\vec{v}_{b}(\\vec{x})||^{2};].
+
+::
+
+ $ cigma compare \
+ -a case64.20000.pvtr:velocity \
+ -b case32.9900.pvtr:velocity \
+ -o steady-state.h5:/error_velocity_64_32
+
+ $ cigma compare \
+ -a case64.20000.pvtr:velocity \
+ -b case16.4700.pvtr:velocity \
+ -o steady-state.h5:/error_velocity_64_16
+
+ $ cigma compare \
+ -a case64.20000.pvtr:velocity \
+ -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
+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
+post-processing utility program called visualize-errors. We can
+accomplish this for our six comparisons, by using the following
+bash comand line,
+
+::
+
+ $ for variable in temperature velocity; do
+ for b in 32 16 8; do
+ visualize-errors \
+ --use-logarithmic-scale \
+ -m case64.20000.vtr \
+ -i steady-state.h5:/error_${variable}_64_${b} \
+ -o log_error_temperature_64_${b}.vtk:log_error_${variable}
+ done
+ done
+
+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,
+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
+gradients and velocity magnitude are at their lowest.
+
+
+.. fig:
+
+
+Likewise, in Figure 5.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:
+
+
+Circular Inclusion Benchmark
+----------------------------
+
+In some cases, you will be able to obtain an exact solution to your
+differential equations. In these cases, you can extend Cigma by
+implementing your function directly in C++, and registering a convenient
+name under which to call it. In this section, we'll show you exactly how
+to accomplish this by considering a two-dimensional problem for which an
+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.
+The kinematic boundary conditions are given generally in terms of the
+simple shear strain rate \dot{\epsilon}.
+
+
+
+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.
+
+
+
+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.
+
+
+
+
+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.
+Since the analytic formula for the velocity field is more complicated than
+the analytic formula for the pressure, we will simply expand the domain and
+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}.
+
+
+
+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
+comparisons,
+
+::
+
+ $ cigma compare \
+ -a 512_8/fields.00000.pvts:PressureField \
+ -b 128_8/fields.00000.pvts:PressureField \
+ -o inclusion.h5:/error_pressure_512_128
+
+ $ cigma compare \
+ -a 512_8/fields.00000.pvts:PressureField \
+ -b 256_8/fields.00000.pvts:PressureField \
+ -o circ_inc.h5:/error_pressure_512_256
+
+
+whose results we summarize in the following Table.
+
+
+
+
+
+from which we can estimate the order of convergence
+\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,
+
+
+
+
+
+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
+
+ cigma compare \
+
+ -a ${r}_8/fields.00000.pvts:PressureField \
+
+ -b bm.circular_inclusion.pressure \
+
+ -o inclusion.h5:/errror_pressure_${r}
+
+ done
+
+whose output we summarize in the following table,
+
+
+
+
+
+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.
+
+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 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