[cig-commits] r12715 - mc/2D/ConMan/trunk/doc
sdk at geodynamics.org
sdk at geodynamics.org
Tue Aug 26 11:31:18 PDT 2008
Author: sdk
Date: 2008-08-26 11:31:18 -0700 (Tue, 26 Aug 2008)
New Revision: 12715
Modified:
mc/2D/ConMan/trunk/doc/conman.tex
Log:
essagenumerous additions corrections and improvements
Modified: mc/2D/ConMan/trunk/doc/conman.tex
===================================================================
--- mc/2D/ConMan/trunk/doc/conman.tex 2008-08-26 18:30:08 UTC (rev 12714)
+++ mc/2D/ConMan/trunk/doc/conman.tex 2008-08-26 18:31:18 UTC (rev 12715)
@@ -79,7 +79,7 @@
\noindent \begin{center}
\thispagestyle{empty}%
\begin{figure}[H]
- \includegraphics[width=0.75\paperwidth]{images/conman_cover}
+ \includegraphics[width=0.75\paperwidth]{./images/conman_cover.pdf}
\end{figure}
\par\end{center}
@@ -190,14 +190,15 @@
\section{The Finite Element Method}
-I am going to follow a similar form to Hughes (The Finite Element
+This section closely follows Hughes (The Finite Element
Method) Chapter~1 sections 1-4. There are two forms we can write
-the equation, the strong and the weak form. You are used to seeing
-the strong form, but not used to seeing the weak form. The finite
+the equation, the strong and the weak form. More readers are probably
+more familiar with
+the strong form, and less familiar with the weak form. The finite
element method, is cast in the weak form. In elasticity, for example,
the weak form comes from a variational principal, such as the principal
-of virtual displacements. For viscous flow, there is also a variational
-form, but we will not go into that.
+of virtual displacements in elasticity. For viscous flow, there is also a variational
+form, but we will not discuss that here.
In general, the finite element method takes a differential equation
(strong form) and transforms it into an integral equation (weak form).
@@ -223,7 +224,7 @@
\begin{equation}
-u,_{x}\left(0\right)=h\label{eq:}\end{equation}
-
+\noindent
This choice of initial conditions allows us to examine both kinds
of boundary conditions. The solution is trivial, but that does not
matter. For completeness, it is \begin{equation}
@@ -241,19 +242,19 @@
\begin{equation}
\intop_{0}^{1}w,_{x}\left(x\right)u,_{x}\left(x\right)=\intop_{0}^{1}w\left(x\right)f\left(x\right)dx+w\left(0\right)h\label{eq:}\end{equation}
-
+\noindent
$\nu$ is the set of weighting functions defined by
\begin{equation}
\nu=\left\{ w\left(x\right)|w\left(x\right)\epsilon H^{1},\, w\left(1\right)=0\right\} \label{eq:}\end{equation}
-
+\noindent
and $\mathcal{L}$ is a set of trial solutions defined by
\begin{equation}
\mathcal{L=\left\{ \mathrm{\mathit{u\left(x\right)|u\left(x\right)}\epsilon H^{1},\,\mathit{u}\left(1\right)=\mathit{g}}\right\} }\label{eq:}\end{equation}
-
+\noindent
$H^{1}$ is the set of all functions whose first derivatives are square
integrable on {[}0, 1]. The integral equation is then solved by integrating
over each element in the domain and adding the result. The result
@@ -262,14 +263,14 @@
\begin{equation}
\left[K\right]x=b\label{eq:}\end{equation}
-
+\noindent
where \emph{{[}K]} is referred to as the element stiffness matrix.
There will be more to say about the implementation in Section \ref{cha:Input-Guide}.
\subsection{Galerkin's Approximation}
-Now we have a start on the finite element method. I want to continue
+Now we have a start on the finite element method. we continue
to follow Hughes; however his notation becomes quite difficult to
keep up with. Now, lets begin to think about putting a solution on
the computer. Because we will have a finite approximation, related
@@ -299,7 +300,7 @@
where the $d_{A}$'s are unknown constants to be solved for.
In the next section we will make the shape functions more concrete.
-I want you to see how general this is, because in principle there
+It is useful to see how general this is, because in principle there
is a great deal of flexibility in how we choose the shape functions.
We have not said anything more about this function $w^{h}$ and how
@@ -361,16 +362,15 @@
\subsection{Shape Functions}
-At this point, I will narrow the focus to deal with specifically the
-elements in my code, ConMan. It is possible to think very general
+At this point, we narrow the focus to deal with specifically the
+elements in ConMan. It is possible to think very general
shape functions, but in practice, people use triangles or quadralaterals
-(in 2-D). We will extend our finite element formulation to a 2-D equation
-next time. In terms of the level of approximation, there are also
+(in 2-D). In terms of the level of approximation, there are also
a lot of possibilities. We will stick to the simplest form, bi-linear
elements, but you should be aware that higher order elements (bi-quadratic
-or bi-cubic spline elements) are also popular with some people. I
-am condensing a lot of very useful material from Chapter 3 of Hughes'
-book into one lecture. If you want to see more complete derivations,
+or bi-cubic spline elements) are also popular with some people.
+This section is condensing a lot of very useful material from Chapter 3 of Hughes'
+book into s short overview. If you want to see more complete derivations,
proofs of convergence, etc., of how to go about using higher order
elements, look at Hughes book, Chapter 3.
@@ -391,19 +391,25 @@
N_{4} & = & \frac{(a-x)(b+y)}{4ab}\end{eqnarray}
these shape functions satisfy the conditions in Equation~\ref{eq:normalize}
and Equations ~\ref{eq:interpx} and~\ref{eq:interpy}. A good exercise
-would be to show this is true.
+would be to show this is true.
+A visual representation of these shape functions is shown in Figure\ref{fig:The-bilinear-shape}.
-\begin{quote}
-\textbf{Problem 7} Show that the shape functions defined above satisfy
-Equation~\ref{eq:normalize} and Equations~\ref{eq:interpx} and~\ref{eq:interpy}.
-\end{quote}
-Notice that by convention, I start numbering my element nodes in the
-lower left hand corner and work counter-clockwise. {\em This is
-a very important point. It is followed through out all my finite element
-codes.} There is no magic reason, you just have to choose a starting
-place.
+\noindent \begin{center}
+%
+\begin{figure}
+\caption{\label{fig:The-bilinear-shape}The bilinear shape function for a single
+element (top) and the four elements whose shape functions combine
+to form the global shape function for node A (bottom). Figure taken
+from Hughes, Sec 3.2.}
-In ConMan, as in Hughes, we further choose to normalize this by setting
+\noindent \begin{centering}
+\includegraphics[scale=0.4]{images/conman-fig2}
+\par\end{centering}
+\end{figure}
+
+\par\end{center}
+
+In ConMan we further choose to normalize this by setting
$a=1$ and $b=1$. This choice gives us an element whose area is 1,
which is a convenient way to think about things. (This is because
we left the factor of 1/4 in the denominator). In my code, to make
@@ -429,7 +435,23 @@
where $J$ is the Jacobian of the transformation. This is a very
powerful point. When we are thinking of solving a regular Cartesian
domain, this just corresponds to a stretching or a shrinking (notice
-we set $a=b=1$ above. However, if we are thinking about a cylindrical
+we set $a=b=1$ above (see Figure \ref{fig:Figure-1}).
+
+\noindent \begin{center}
+%
+\begin{figure}
+\caption{\label{fig:Figure-1}The mapping between the global domain (right)
+and the parent element domain (left) using the shape functions. Figure
+taken from Hughes, Sec 3.2.}
+
+
+\noindent \begin{centering}
+\includegraphics[scale=0.45]{images/conman-fig1}
+\par\end{centering}
+\end{figure}
+
+\par\end{center}
+However, if we are thinking about a cylindrical
geometry, for example, we can use the Jacobian of the transformation
between the geometries. Lets look at two examples:
@@ -512,8 +534,8 @@
routine once and stored all the shape functions. As problem sizes
have grown, this took a lot of storage, so now we call it on the fly
for each element when needed. It is computationally more expensive
-but cuts storage. This strategy will also be necessary for a Lagrangian
-formulation or adaptive gridding.
+but cuts the storage requirement. This strategy will also be necessary
+for a Lagrangian formulation or adaptive gridding.
There are two domains to keep in mind when thinking about the finite
element method: the global domain and the parent element domain (Figure
@@ -526,21 +548,6 @@
between the input domain and the parent element domain, which is calculated
in routine \textbf{genshg}.
-\noindent \begin{center}
-%
-\begin{figure}
-\caption{\label{fig:Figure-1}The mapping between the global domain (right)
-and the parent element domain (left) using the shape functions. Figure
-taken from Hughes, Sec 3.2.}
-
-
-\noindent \begin{centering}
-\includegraphics[scale=0.45]{images/conman-fig1}
-\par\end{centering}
-\end{figure}
-
-\par\end{center}
-
For ConMan the choice was made to use bilinear quadrilaterals as the
parent elements (Figure \ref{fig:The-bilinear-shape}). Higher order
elements (i.e., biquadratic or bicubic-spline) require more computational
@@ -548,26 +555,15 @@
rather than using high-order elements, is the best strategy for an
efficient, accurate code for incompressible, advection-diffusion problems.
-\noindent \begin{center}
-%
-\begin{figure}
-\caption{\label{fig:The-bilinear-shape}The bilinear shape function for a single
-element (top) and the four elements whose shape functions combine
-to form the global shape function for node A (bottom). Figure taken
-from Hughes, Sec 3.2.}
-
-\noindent \begin{centering}
-\includegraphics[scale=0.4]{images/conman-fig2}
-\par\end{centering}
-\end{figure}
-
-\par\end{center}
-
Because of the changing between domains, it is necessary to define
several bookkeeping arrays to identify nodes and elements in each
of the domains.
+It is worth noting that the numbering of local elements always begins with
+1 in the lower left hand corner. There is no special reason, you just have
+to choose a convention.
+
\begin{description}
\item [{{id}}] transforms global nodes to equation numbers (Figure \ref{fig:Example-relationship-between}).
\item [{{ien}}] transforms element local node numbers to global node
@@ -636,12 +632,23 @@
examples:
Consider a 3 element by 3 element grid. I will number my elements
-and node starting in the lower left hand corner, and working up fastest.
+and node starting in the lower left hand corner, and the
+node and element numbers will increase linearly in the vertical direction.
-\begin{verbatim} For element 3: ien (3, 1) = 3 ien (3, 2) = 7 ien
-(3, 3) = 8 ien (3, 4) = 4 For element 5: ien (5, 1) = 6 ien (5, 2)
-= 10 ien (5, 3) = 11 ien (5, 4) = 7 \end{verbatim}
+\begin{verbatim}
+For element 3:
+ien (3, 1) = 3
+ien (3, 2) = 7
+ien (3, 3) = 8
+ien (3, 4) = 4
+For element 5:
+ien (5, 1) = 6
+ien (5, 2) = 10
+ien (5, 3) = 11
+ien (5, 4) = 7
+\end{verbatim}
+
There are three kinds of opperations we might think of wanting. The
first is taking values of some function (coordinates, velocities,
temperatures, stresses, etc.) defined on the global grid and getting
@@ -652,8 +659,8 @@
to the global value for that node, an {\em assembly} step.
All three operations, gather, scatter, and assemble are done by the
-routine local. Because it is called by the genshp routine above, it
-is a good example to look at.
+routine {\bf local}. Because it is called by the {\bf genshp} routine above, it
+is a good example to study.
\subsection{Equations}
@@ -799,9 +806,15 @@
for the temperature equation is stable, then this method is stable
and converges as $\Delta t\rightarrow0$.
-The element stiffness matrix (Figure \ref{fig:5-The-element-stiffness})
+The element stiffness matrix (equation 2.84)
is made up of the two terms from the left hand side of the integral
-equation. The integration is done using two by two gauss quadrature,
+equation.
+The full element stiffness matrix for the quadralilateral element
+is an 8 by 8 matrix made up of 16 of the 2 by 2 matrices as shown
+in Figure~\ref{fig:8-The-storage-for}. Because of symmetry
+we only need to form and store the upper triangular part of the
+matrix.
+The integration is done using two by two gauss quadrature,
which is exact when the elements are rectangular and bilinear shape
functions are used. The $\lambda$ term is under-integrated (one point
rule) to keep the large penalty value from effectively locking the
@@ -819,24 +832,24 @@
The stiffness matrix is formed in routine \textbf{f\_vstf} and the
right hand side is formed in routine \textbf{f\_tres}.
-\noindent \begin{center}
+%\noindent \begin{center}
%
-\begin{figure}
-\caption{\label{fig:5-The-element-stiffness}The element stiffness matrix for
-the 2D Cartesian stokes equation. The 8 by 8 matrix is make up of
-16 2 by 2 submatrices of the form shown below. The $\lambda$ and
-$\mu$ parts are shown separately for clarity.}
+%\begin{figure}
+%\caption{\label{fig:5-The-element-stiffness}The element stiffness matrix for
+%the 2D Cartesian stokes equation. The 8 by 8 matrix is make up of
+%16 2 by 2 submatrices of the form shown below. The $\lambda$ and
+%$\mu$ parts are shown separately for clarity.}
+%
+%
+%\noindent \begin{centering}
+%\includegraphics[scale=0.4]{images/conman-fig5}
+%\par\end{centering}
+%\end{figure}
+%
+%\par\end{center}
+%
+%I'm not sure if we need this too...
-
-\noindent \begin{centering}
-\includegraphics[scale=0.4]{images/conman-fig5}
-\par\end{centering}
-\end{figure}
-
-\par\end{center}
-
-I'm not sure if we need this too...
-
\noindent \begin{center}
%
\begin{figure}
@@ -867,7 +880,7 @@
\begin{equation}
T_{,j}n_{j}=q\,\,\,\,\,\, on\,\,\Gamma_{q}\label{eq:}\end{equation}
-
+\noindent
where $T$ is the temperature, $u_{i}$ is the velocity, $\kappa$
is the thermal diffusivity and $H$ is the internal heat source. The
weak form of the energy equation is given by
@@ -879,7 +892,7 @@
\begin{equation}
-\kappa\int_{\Omega}w_{,i}T_{,i}d\Omega+\int_{\Gamma_{q}}wT_{,j}n_{j}d\Gamma_{q}\label{eq:}\end{equation}
-
+\noindent
where $\dot{T}$ is the time derivative of temperature, $T_{,i}$
is the gradient of temperature, $w$ is the standard weighting function
and $\left(w+p\right)$ is the Petrov-Galerkin weighting function
@@ -901,7 +914,7 @@
\begin{equation}
\left(\xi u_{\xi}h_{\xi}+\eta u_{\eta}h_{\eta}\right)/2\label{eq:}\end{equation}
-
+\noindent
with
\begin{equation}
@@ -911,7 +924,7 @@
\begin{equation}
\eta=1-\frac{2\kappa}{u_{\eta}h_{\eta}}\label{eq:}\end{equation}
-
+\noindent
where $h_{\xi}$ and $h_{\eta}$ are the element lengths and $u_{\xi}$
and $u_{\eta}$ are the velocities in the local element coordinate
system ($\xi$ $\eta$ system) evaluated at the element center. This
@@ -971,7 +984,7 @@
\begin{equation}
\dot{T}_{n+1}^{\left(i+1\right)}=\dot{T}_{n+1}^{\left(i\right)}+\Delta\dot{T}_{n+1}^{\left(i\right)}\label{eq:}\end{equation}
-
+\noindent
where $i$ is the iteration number (for the corrector), $n$ is the
time step number, $T$ is the temperature, $\dot{T}$ is the derivative
of temperature with time, $\Delta\dot{T}$ is the correction to the
@@ -1478,19 +1491,26 @@
\item [{{dmhu(i)}}] . . . . . . internal heat source for ith material
group
\end{description}
+
\item [{{Activation~Energy~Card}}] \emph{tcon(1,i), i=1,numat}
+\begin{description}
+\item [{{tcon(1,i)}}] . . . . . activation energy for
+ith material group for temperature dependent viscosity (kJ/mole)
+\end{description}
+\item [{{Activation~Volume~Card}}] \emph{tcon(2,i), i=1,numat}
\begin{description}
-\item [{{tcon(1,i)}}] . . . . . nondimensional activation energy for
-ith material group for temperature dependent viscosity
+\item [{{tcon(2,i)}}] . . . . . activation volume for
+ith material group for temperature dependent viscosity (cm$^3$/mole)
\end{description}
-\item [{{Temperature~Dependent~Viscosity~Temperature~Offset~Card}}] \emph{tcon(2,i),
-i=1,numat}
+\item [{{Viscosity Cutoff~Card}}] \emph{tcon(3,i), i=1,numat}
\begin{description}
-\item [{{tcon(2,i)}}] . . . . . nondimen temperature offset for ith material
-group for temperature dependent viscosity
+\item [{{tcon(3,i)}}] . . . . . maximum value of the viscosity
+for the ith material group
\end{description}
+
+
\item [{{Surface~Force/Flux~Cards}}] - numsuf cards \emph{nel side
fnorm ftan flux}
@@ -1723,27 +1743,181 @@
a 1 by 1 square, constant viscosity with the Picard method. This is
Blankenbach 1a
-\begin{lyxcode}
-50~x~50~el.~plate~problem~from~Blankenbach~et~al.~1989~\#Nds~sdm~dof~X~Z~prc~ck~echo~rrst~wrst~nus~tdvf~tseq~nelg~sky~wr~2601~2~2~50~50~2~1~0~0~1~102~0~1~1~1~0~time~step~information~100~1~1.0~1.0~0.50000~output~information~100~100~100~100~velocity~boundary~condition~flags:~IFCMT,DELNXTLN~bnode~enode~incr~bcf1~bcf2~1~2551~51~0~1~2551~2601~1~1~0~1~51~1~1~0~51~2601~51~0~1~1~1~1~1~1~51~51~1~1~1~2551~2551~1~1~1~2601~2601~1~1~1~0~0~0~0~0~temperature~boundary~condition~flags~1~2551~51~1~51~2601~51~1~0~0~0~0~bndy~info~(top~-~bottom~rows)~1~2551~51~51~2601~51~0~0~0~bndy~info~(2nd~from~top~-~2nd~from~bottom~rows)~2~2552~51~50~2600~51~0~0~0~initial~condition~information~0.1~1.0~1.0~1.0~element~information~2~2500~4~4~1~2~0~5~0~0~viscosity~1.0e0~penalty~number~0.1E+08~diffusivity~(always~one)~1.0~Rayleigh~number~1.0e+04~internal~heating~parameter~0.0~0.0~0.0
-\end{lyxcode}
+\begin{verbatim}
+50 x 50 el. problem for Blankenbach Benchmark
+numnp nsd ndof nelx nelz prc ck iecho rrst wrst nus tdvf tseq nelg sky wr
+ 2601 2 2 50 50 2 1 0 0 1 102 0 1 1 1 0
+ time step information
+ 100 1 1.0 1.0 0.50000
+ output information
+ 100 100 100 100
+ velocity boundary condition flags
+ bnode enode incr bcf1 bcf2
+ 1 2551 51 0 1
+ 2551 2601 1 1 0
+ 1 51 1 1 0
+ 51 2601 51 0 1
+ 1 1 1 1 1
+ 51 51 1 1 1
+ 2551 2551 1 1 1
+ 2601 2601 1 1 1
+ 0 0 0 0 0
+ temperature boundary condition flags
+ 1 2551 51 1
+ 51 2601 51 1
+ 0 0 0 0
+ bndy info (top - bottom rows)
+ 1 2551 51
+ 51 2601 51
+ 0 0 0
+ bndy info (2nd from top - 2nd from bottom rows)
+ 2 2552 51
+ 50 2600 51
+ 0 0 0
+ initial condition information
+ 0.1 1.0 1.0 1.0
+ element information
+ 2 2500 4 4 1 2 0 5 0 0
+ viscosity
+ 1.0e0
+ penalty number
+ 0.1E+08
+ diffusivity (always one)
+ 1.0
+ Rayleigh number
+ 1.0e+04
+ internal heating parameter
+ 0.0
+ activation energy
+ 0.0
+ activation volume
+ 0.0
+ viscosity cutoff
+ 1.0e7
+\end{verbatim}
The lines below are a sample geometry file for the 50 by 50 element
problem.
-\begin{lyxcode}
-coordinates~1~4~0.0~0.0~2551~1~1.0~0.0~2601~1~1.0~1.0~51~1~0.0~1.0~50~51~50~1~0~0~0.0~0.0~velocity~boundary~conditions~(non-zero)~0~0~0.0~0.0~temperature~boundary~conditions~(non-zero)~1~2~1.0~2551~0~1.0~50~51~0~0~0.0~element~connectivity~and~material~groups~1~1~1~1~52~53~2~50~50~51~50~1~1~0~0~0~0~0~0~0
-\end{lyxcode}
+\begin{verbatim}
+coordinates
+ 1 4 0.0 0.0
+ 2551 1 1.0 0.0
+ 2601 1 1.0 1.0
+ 51 1 0.0 1.0
+ 50 51 50 1
+ 0 0 0.0 0.0
+ velocity boundary conditions (non-zero)
+ 0 0 0.0 0.0
+ temperature boundary conditions (non-zero)
+ 1 2 1.0
+ 2551 0 1.0
+ 50 51
+ 0 0 0.0
+ element connectivity and material groups
+ 1 1 1 1 52 53 2
+ 50 50 51 50 1 1
+ 0 0 0 0 0 0 0
+\end{verbatim}
\chapter{\label{cha:Output-Guide}Output Guide}
\section{The Output Files}
+The output files names are taken from the names in the runfile.
+The execution of conman proceeds by
+\noindent
+\% conman.pic < runfile
+\noindent
+where the runfile has a list of filenames that can be up to 80 characters long.
+The names in the runfile are attached to the following input or output files.
+
+\begin{verbatim}
+input
+geometry
+output
+restart input
+restart output
+time series output
+temperature (and velocity) field output
+stress (and viscosity) field output
+geoid output
+\end{verbatim}
+
+All files are ascii files. The {\em input} and {\em geometry} files are
+as described in the previous section.
+
+The {\em output} file is a formatted record of the input. If the variable
+{\em necho} is set to one, then values of the coordinates and boundary conditions
+are output and the file can become large. Near the end of the {\em output}
+file execution time is listed for various subparts of the code.
+
+The {\em restart input} file is an input file that is used if the variable
+{\em inrstr} is set to one.
+The file is ascii but formatted and the format statement can be found
+in the file {\bf input.F}.
+The first line contains the initial timestep and time,
+the second line is a header, and the third through {\em numnp} lines
+contains the node, temperature and time derivative of temperature.
+
+The {\em restart output} file is an output file that is used if the variable
+{\em iorstr} is set to one. We recommend this always be set to run.
+This file is overwritten every {\em nsdprt} steps and the output
+is written from the file {\bf timdrv.F}
+
+The {\em time series output} file is written every time step.
+The values are ascii and are heat flux at the bottom, heat flux
+at the top, kinetic energy, and time. All values are dimensionless.
+This file is written from the routine {\bf fluxke.F}
+
+The {\em Temperature and velocity output} file is an ascii file
+written from the routine {\bf print.F}. There are two header
+lines that contain, {\em nsd, nelx, nelz, numnp, nstep, time}.
+The second line labels the output and the third through
+numnp lines list the node, x, z, vx, vz, and temperature values at each node.
+This file is output every {\em nsvprt} steps and each successive
+set of values is appended to the end of the file.
+The unix command {\em split -{numnp+2} tempfile} can be used
+to split the file into files that are {\em numnp+2} lines long.
+
+The {\em Stress output} file is an ascii file
+written from the routine {\bf prtstr} which can be found in
+the file {\bf stress.F}. Like the temperature velocity file
+there are two header
+lines that contain, {\em nsd, nelx, nelz, numnp, nstep, time}.
+The second line labels the output and the third through
+numnp lines list the node, x, z, txx, tzz. txz. p and viscosity at each node.
+This file is output every {\em nsvprt} steps and each successive
+set of values is appended to the end of the file.
+The unix command {\em split -{numnp+2} tempfile} can be used
+to split the file into files that are {\em numnp+2} lines long.
+
+The {\em geoid output} file contains the horizontal coordinate,
+the dynamic topography, the geoid contribution from the temperature
+only and the total geoid (from surface and cmb topography, and
+internal (i.e, temperature) density contrasts).
+This is written from the file {\bf geoid.F}.
+The dimensional values are hard wired into the code and can
+be found at the top of geoid.F
+This version of the geoid code has not been tested with flow-through
+boundary conditions (i.e, {\em nwrap = nelz}).
+
+{\em Warning. As written, the code assumes a Rayleigh number of
+$10^7$ and does not adjust the parameters as the Rayleigh number
+changes. This would be a trivial fix but still may not yield the results
+intended. It is best to carefully check the {\bf geoid.F} source for
+the specific problem if interest.}
+
+
\chapter{\label{cha:The-Benchmark-Cases}The Benchmark Cases}
+We provide input and geometry files for the benchmark cases described
+below in the directories cookbook1 and cookbook2.
-\section{Constant Viscosity Benchmark for ConMan}
+\section{Cookbook 1: Constant Viscosity Benchmark for ConMan}
+
Here we reproduce the results from Blankenbach et al. (1989) for constant
viscosity in a unit-aspect ratio domain, with free-slip boundary conditions,
heated from below and cooled from above. The user should note that
@@ -1872,7 +2046,7 @@
-\section{Temperature-Dependent Viscosity Benchmark for ConMan}
+\section{Cookbook 2: Temperature-Dependent Viscosity Benchmark for ConMan}
Here we reproduce the results from Blankenbach et al. (1989) for temperature-dependent
viscosity in a unit-aspect ratio domain, with free-slip boundary conditions,
@@ -1913,7 +2087,35 @@
\end{table}
+\section{Cookbook 3: Constant Viscosity Driven Slab Problem for ConMan}
+This geometry is based on the subduction zone benchmark by van Keken et al. (2008).
+A thermal structure from this is presented on the cover. The purpose here is to illustrate
+how one can (with some struggle) implement a deformed and variable mesh with the
+fairly basic grid generator in ConMan. One can use a more advanced grid generator
+and input the coordinates of every node directly as well as the ien array for each
+element.
+
+This problem will not reproduce the exact result presented in van Keken et al (2008)
+because this version of ConMan does not have the discontinuous velocity field, which
+turns out to be critical to match the benchmark results.
+The resulting thermal structure should look like the Figure~\ref{fig:wedge}.
+
+\noindent \begin{center}
+%
+\begin{figure}
+\caption{\label{fig:wedge} Thermal structure for the wedge problem
+in Cookbook 3 example.}
+
+
+\begin{centering}
+\includegraphics[scale=0.6]{images/wedge}
+\par\end{centering}
+\end{figure}
+
+\par\end{center}
+
+
\appendix
%dummy comment inserted by tex2lyx to ensure that this paragraph is not empty
@@ -2308,10 +2510,13 @@
\bibitem[9]{Travis et al 1991}Travis, B.J., C. Anderson, J. Baumgardner,
C. Gable, B.H. Hager, P. Olson, R.J. O\textquoteright{}Connell,
-A. Raefsky, and G. Schubert, 1991. A benchmark comparison of numerical
+A. Raefsky, and G. Schubert, 1990. A benchmark comparison of numerical
methods for infinite Prandtl number convection in two-dimensional
-Cartesian geometry, \emph{Geophys. Astrophys. Fluid Dynamics,} in
-press.
+Cartesian geometry, \emph{Geophys. Astrophys. Fluid Dynamics,}
+55, 137-160.
+
+\bibitem[9]{van Keken et al 2008} van Keken, P.E., Currie, C., King, S.D., Behn, M.D., Canioncle, A., He, J., Katz, R.F., Lin, S.-C., Parmentier, E.M., Spiegelman, M., and Wang, K., 2008. A community benchmark for subduction zone modeling, Phys. Earth Planet. Int., in press.
+
\end{thebibliography}
\end{document}
More information about the cig-commits
mailing list