[cig-commits] r15987 - doc/geodynamics.org/benchmarks/trunk/cigma
luis at geodynamics.org
luis at geodynamics.org
Mon Nov 16 14:50:22 PST 2009
Author: luis
Date: 2009-11-16 14:50:21 -0800 (Mon, 16 Nov 2009)
New Revision: 15987
Added:
doc/geodynamics.org/benchmarks/trunk/cigma/index.header
doc/geodynamics.org/benchmarks/trunk/cigma/index.html
Log:
Generated cigma/index.html
Added: doc/geodynamics.org/benchmarks/trunk/cigma/index.header
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/cigma/index.header (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/cigma/index.header 2009-11-16 22:50:21 UTC (rev 15987)
@@ -0,0 +1,17 @@
+id: index_html
+title: Cigma
+subject:
+description: Using Cigma
+contributors:
+creators: luis
+effectiveDate: None
+expirationDate: None
+language:
+rights:
+creation_date: 2009/10/18 02:58:59.058 GMT-7
+modification_date: 2009/10/21 12:08:04.944 GMT-7
+excludeFromNav: False
+relatedItems:
+allowDiscussion: None
+Content-Type: text/x-rst
+
Added: doc/geodynamics.org/benchmarks/trunk/cigma/index.html
===================================================================
--- doc/geodynamics.org/benchmarks/trunk/cigma/index.html (rev 0)
+++ doc/geodynamics.org/benchmarks/trunk/cigma/index.html 2009-11-16 22:50:21 UTC (rev 15987)
@@ -0,0 +1,492 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
+<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
+<meta name="generator" content="Docutils 0.5: http://docutils.sourceforge.net/" />
+<title></title>
+<script type="text/javascript" src="http://geodynamics.org/cig/portal_javascripts/Plone%20Default/textheworld6.user.js"></script>
+<link rel="stylesheet" href="../css/voidspace.css" type="text/css" />
+</head>
+<body>
+<div class="document">
+
+
+<p>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.</p>
+<div class="section" id="poisson-problem">
+<h1>Poisson Problem</h1>
+<p>Cigma can compare two functions and return the <em>L</em><sub>2</sub>-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:</p>
+<ol class="arabic simple">
+<li>[;\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</li>
+<li>[;\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;].</li>
+</ol>
+<p>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
+<a class="reference external" href="http://dealii.org/developer/doxygen/tutorial/">http://dealii.org/developer/doxygen/tutorial/</a>.</p>
+<p>In the <tt class="docutils literal"><span class="pre">examples/</span></tt> subdirectory, we include a version of the
+Step 4 tutorial program that has been slightly modified to suit
+our purposes in this section.</p>
+<p>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×64, 32×32, 16×16, 8×8
+by using linear quadrilateral elements.
+The solution is saved in file <tt class="docutils literal"><span class="pre">phi_64x64x64.vtk</span></tt> for the highest resolution,
+for example. In a similar fashion, we solve our 3D Poisson problem on meshes
+of resolution 64×64×64, 32×32×32, 16×16×16, 8×8×8,
+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</p>
+<pre class="literal-block">
+$ 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
+</pre>
+<p>When we run the above code on the command line, or script, Cigma will generate
+a summary of each comparison, including the <em>L</em><sub>2</sub> 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:</p>
+<!-- table:1 -->
+<table border="1" class="docutils">
+<colgroup>
+<col width="29%" />
+<col width="32%" />
+<col width="39%" />
+</colgroup>
+<tbody valign="top">
+<tr><td>[;n;]</td>
+<td>[;h_n;]</td>
+<td>[;||\phi_n - \phi_{64\times 64}||_{L_2} / \sqrt{V};]</td>
+</tr>
+<tr><td>8×8</td>
+<td>0.35355</td>
+<td>0.012312</td>
+</tr>
+<tr><td>16×16</td>
+<td>0.17677</td>
+<td>0.003157</td>
+</tr>
+<tr><td>32×32</td>
+<td>0.08838</td>
+<td>0.000689</td>
+</tr>
+</tbody>
+</table>
+<!-- table:1 -->
+<table border="1" class="docutils">
+<colgroup>
+<col width="32%" />
+<col width="27%" />
+<col width="41%" />
+</colgroup>
+<tbody valign="top">
+<tr><td>[;n;]</td>
+<td>[;h_n;]</td>
+<td>[;||\phi_n - \phi_{64\times 64\times 64}||_{L_2} / \sqrt{V};]</td>
+</tr>
+<tr><td>8×8×8</td>
+<td>0.43301</td>
+<td>0.019205</td>
+</tr>
+<tr><td>16×16×16</td>
+<td>0.21650</td>
+<td>0.004889</td>
+</tr>
+<tr><td>32×32×32</td>
+<td>0.10825</td>
+<td>0.001075</td>
+</tr>
+</tbody>
+</table>
+<p>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:</p>
+<!-- fig:1 -->
+<div align="center" class="figure">
+<img alt="images/bm1_alpha.png" src="images/bm1_alpha.png" style="width: 800px;" />
+<p class="caption"><span class="target" id="figure-1">Figure 1</span>: Convergence of <em>L</em><sub>2</sub> Global Error for Poisson
+problems in 2D and 3D.</p>
+</div>
+<p>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.</p>
+</div>
+<div class="section" id="mantle-convection">
+<h1>Mantle Convection</h1>
+<p>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 <em>L</em><sub>2</sub> 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.</p>
+<p>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<sup>5</sup> 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.</p>
+<p>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×64×64, 32×32×32, 16×16×16,
+and 8×8×8 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.</p>
+<!-- fig:2 -->
+<div align="center" class="figure">
+<img alt="images/bm2_out_fig_fields_64_20000_y099.png" src="images/bm2_out_fig_fields_64_20000_y099.png" style="width: 800px;" />
+<p class="caption"><span class="target" id="figure-2">Figure 2</span>: 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.</p>
+</div>
+<p>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 (<tt class="docutils literal"><span class="pre">*.vtr</span></tt> files),
+and a subsequent <tt class="docutils literal"><span class="pre">*.pvtr</span></tt> file for each solution step.</p>
+<p>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 <em>L</em><sub>2</sub>-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 <cite>Figure 3</cite>,</p>
+<!-- fig:3 -->
+<div align="center" class="figure">
+<img alt="images/bm2_heatflux.png" src="images/bm2_heatflux.png" style="width: 800px;" />
+<p class="caption"><span class="target" id="figure-3">Figure 3</span>: Time variation of average heat flux over top surface.</p>
+</div>
+<p>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.:</p>
+<pre class="literal-block">
+$ cigma compare \
+ -a case64.20000.pvtr:temperature \
+ -b case32.9900.pvtr:temperature \
+ -o steady-state.h5:/error_temperature_64_32
+</pre>
+<p>Since we did not specify an integration mesh through the <tt class="docutils literal"><span class="pre">-m</span></tt> option,
+Cigma will take the mesh from the first field and use it for the integration
+of the <em>L</em><sub>2</sub> error. We have, however, used the <tt class="docutils literal"><span class="pre">-o</span></tt> option to specify
+an output file. This will create the HDF5 file <tt class="docutils literal"><span class="pre">steady-state.h5</span></tt> after
+the comparison is complete, and store the results of that comparison into
+an array called <tt class="docutils literal"><span class="pre">error_temperature_64_32</span></tt> 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,</p>
+<pre class="literal-block">
+$ 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
+</pre>
+<p>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};].</p>
+<pre class="literal-block">
+$ 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
+</pre>
+<p>At this point, you have obtained the local <em>L</em><sub>2</sub>-errors for each
+of the above comparisons, and have stored them in arrays inside the
+file <tt class="docutils literal"><span class="pre">steady-state.h5</span></tt>. 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,</p>
+<pre class="literal-block">
+$ 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
+</pre>
+<p>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 <em>L</em><sub>2</sub> 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.</p>
+<p>In <cite>Figure 4</cite>, we show two plots representing the errors in the temperature
+field on the y=0.99 plane for the two 16×16×16 and 32×32×32
+resolutions, relative to the 64×64×64 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.</p>
+<!-- fig:4 -->
+<div align="center" class="figure">
+<img alt="images/bm2_out_fig_log_error_temperature_64_16_64_32_y099.png" src="images/bm2_out_fig_log_error_temperature_64_16_64_32_y099.png" style="width: 800px;" />
+<p class="caption"><span class="target" id="figure-4">Figure 4</span>: Temperature differences on y=0.99 plane, as compared
+against case with 64×64×64 elements.</p>
+</div>
+<p>Likewise, in <cite>Figure 5</cite> 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.</p>
+<!-- fig:5 -->
+<div align="center" class="figure">
+<img alt="images/bm2_out_fig_log_error_velocity_64_16_64_32_y099.png" src="images/bm2_out_fig_log_error_velocity_64_16_64_32_y099.png" style="width: 800px;" />
+<p class="caption"><span class="target" id="figure-5">Figure 5</span>: Velocity differences on y=0.99 plane, as compared
+against case with 64×64×64 elements.</p>
+</div>
+</div>
+<div class="section" id="circular-inclusion-benchmark">
+<h1>Circular Inclusion Benchmark</h1>
+<p>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};].</p>
+<!-- fig:6 -->
+<div align="center" class="figure">
+<img alt="images/bm3_inclusion_setup.png" src="images/bm3_inclusion_setup.png" />
+<p class="caption"><span class="target" id="figure-6">Figure 6</span>: Schematic for the circular inclusion benchmark.</p>
+</div>
+<p>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.</p>
+<p>From [Schmid Podladchikov 2003], we end up with the following analytic
+formula for the pressure field under the case of simple shear,</p>
+<blockquote>
+[;
+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}
+;]</blockquote>
+<p>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 <tt class="docutils literal"><span class="pre">bm.circular_inclusion.pressure</span></tt>,
+whose definition is summarized in <cite>Listing 1</cite>.</p>
+<!-- listing:1 -->
+<pre class="literal-block">
+#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;
+ }
+};
+</pre>
+<p>One additional step must be taken before we can use the name
+<tt class="docutils literal"><span class="pre">bm.circular_inclusion.pressure</span></tt> 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 <tt class="docutils literal"><span class="pre">FunctionRegistry.cpp</span></tt> source file.</p>
+<!-- listing:2 -->
+<pre class="literal-block">
+FunctionRegistry::FunctionRegistry()
+{
+ //
+ // ... other definitions
+ //
+
+ typedef benchmark::circular_inclusion::Pressure PressureFn1;
+ shared_ptr<PressureFn1> pressure1(new PressureFn1());
+ this->addFunction(“bm.circular_inclusion.pressure”, pressure1);
+}
+</pre>
+<p>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,</p>
+<blockquote>
+[;
+\begin{align*}
+v_x &= -\dot{\epsilon}y, \\
+v_y &= \dot{\epsilon}x,
+\end{align*}
+;]</blockquote>
+<p>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};].</p>
+<!-- fig:7 -->
+<div align="center" class="figure">
+<img alt="images/bm3_out_fig_pressure_512.png" src="images/bm3_out_fig_pressure_512.png" style="width: 800px;" />
+<p class="caption"><span class="target" id="figure-7">Figure 7</span>: Pressure field for circular inclusion problem,
+with resolution of 512×512 elements.</p>
+</div>
+<p>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,</p>
+<pre class="literal-block">
+$ 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
+</pre>
+<p>whose results we summarize in the following Table.</p>
+<!-- table:p512 -->
+<p>from which we can estimate the order of convergence</p>
+<blockquote>
+[;
+\alpha=\log(0.00158/0.00407)/\log(0.0442/0.0884)=1.36.
+;]</blockquote>
+<p>Processing the local error fields stored in <tt class="docutils literal"><span class="pre">circ_inc.h5</span></tt>,
+we obtain the following plots,</p>
+<!-- fig:8 -->
+<div align="center" class="figure">
+<img alt="images/bm3_out_fig_full_log_error_pressure_512_128_512_256.png" src="images/bm3_out_fig_full_log_error_pressure_512_128_512_256.png" style="width: 800px;" />
+</div>
+<!-- fig:9 -->
+<div align="center" class="figure">
+<img alt="images/bm3_out_fig_inc_log_error_pressure_512_128_512_256.png" src="images/bm3_out_fig_inc_log_error_pressure_512_128_512_256.png" style="width: 800px;" />
+</div>
+<p>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,</p>
+<pre class="literal-block">
+$ 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
+</pre>
+<p>whose output we summarize in the following table,</p>
+<!-- table -->
+<p>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.</p>
+<p>After processing the errors in inclusion.h5 with the utility
+program visualize-errors, we obtain the following plots of the
+pressure field differences,</p>
+<!-- fig:10 -->
+<div align="center" class="figure">
+<img alt="images/bm3_out_fig_inc_log_error_pressure_256_512.png" src="images/bm3_out_fig_inc_log_error_pressure_256_512.png" style="width: 800px;" />
+</div>
+<!-- fig:11 -->
+<div align="center" class="figure">
+<img alt="images/bm3_out_fig_inc_log_error_pressure_256_512.png" src="images/bm3_out_fig_inc_log_error_pressure_256_512.png" style="width: 800px;" />
+</div>
+<p>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 <em>L</em><sub>2</sub> 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.</p>
+<p>In summary, we have observed that for this problem the
+<em>L</em><sub>2</sub> 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.</p>
+</div>
+</div>
+</body>
+</html>
More information about the CIG-COMMITS
mailing list