From d4ff8af67d1d47441a0a3b25735a990b3bc9de3a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 17 Aug 2009 13:51:21 +0000 Subject: [PATCH] Discuss end times. git-svn-id: https://svn.dealii.org/trunk@19284 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/Makefile | 6 ++-- deal.II/examples/step-32/doc/intro.dox | 42 ++++++++++++++------------ deal.II/examples/step-32/step-32.cc | 4 +-- 3 files changed, 28 insertions(+), 24 deletions(-) diff --git a/deal.II/examples/step-32/Makefile b/deal.II/examples/step-32/Makefile index 94a9f0780f..770b78a94b 100644 --- a/deal.II/examples/step-32/Makefile +++ b/deal.II/examples/step-32/Makefile @@ -14,7 +14,7 @@ target = $(basename $(shell echo step-*.cc)) # run-time checking of parameters and internal states is performed, so # you should set this value to `on' while you develop your program, # and to `off' when running production computations. -debug-mode = on +debug-mode = off # As third field, we need to give the path to the top-level deal.II @@ -45,13 +45,13 @@ clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk # settings include $D/common/Make.global_options -################################################################ +################################################################ # This example program will only work if Trilinos is installed. If this # is not the case, then simply redefine the main targets to do nothing ifneq ($(USE_CONTRIB_TRILINOS),yes) default run clean: @echo - @echo "===========================================================" + @echo "===========================================================" @echo "= This program cannot be compiled without Trilinos. Make=" @echo "= sure you have Trilinos installed and detected during =" @echo "= configuration of deal.II =" diff --git a/deal.II/examples/step-32/doc/intro.dox b/deal.II/examples/step-32/doc/intro.dox index 240a1ad131..30c5baf428 100644 --- a/deal.II/examples/step-32/doc/intro.dox +++ b/deal.II/examples/step-32/doc/intro.dox @@ -333,7 +333,7 @@ SI system) even if we list them here in other units. As a reminder, the equations we want to solve are these: @f{eqnarray*} - -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& + -\nabla \cdot (2 \eta \varepsilon ({\mathbf u})) + \nabla p &=& -\rho \; \beta \; T \mathbf{g}, \\ \nabla \cdot {\mathbf u} &=& 0, @@ -355,7 +355,7 @@ the following quantities: continental and oceanic - plates) to a depth of 2890km (where the + plates) to a depth of 2890km (where the outer earth core starts). The radii are therefore $R_0=(6371-2890)\text{km}, R_1=(6371-35)\text{km}$. This domain is conveniently generated using the @@ -408,8 +408,8 @@ the following quantities: proportional to the velocity difference, leading to a boundary condition of the form @f{align*} - \mathbf{n}\cdot [2\eta \varepsilon(\mathbf v)] - &= + \mathbf{n}\cdot [2\eta \varepsilon(\mathbf v)] + &= s \mathbf{n} \times [\mathbf v - \mathbf v_0], \\ \mathbf{n} \cdot \mathbf v &= 0, @@ -429,30 +429,30 @@ the following quantities: \hat{\mathbf v}$. That's certainly not a good situation that we would like to avoid. The traditional way to work around this is to pick an arbitrary point on the boundary and call this your fixed point by choosing the - velocity to be zero in all components there. + velocity to be zero in all components there.
  • To first order, the gravity vector always points downward. The question for a body as big as the earth is just: where is "up". The answer of course is "radially inward, towards the center of the earth". So at the surface of the earth, we have @f[ - \mathbf g - = + \mathbf g + = -9.81 \frac{\text{m}}{\text{s}^2} \frac{\mathbf x}{\|\mathbf x\|}, @f] where $9.81 \frac{\text{m}}{\text{s}^2}$ happens to be the average gravity acceleration at the earth surface. But in the earth interior, the question becomes a bit more complicated: at the center of the earth, for example, you have matter pulling equally hard in all directions, and so $\mathbf g=0$. In - between, the net force is described as follows: let us define the + between, the net force is described as follows: let us define the gravity potential by @f[ \varphi(\mathbf x) = - \int_{\text{earth}} - -G \frac{\rho(\mathbf y)}{\|\mathbf x-\mathbf y\|} + \int_{\text{earth}} + -G \frac{\rho(\mathbf y)}{\|\mathbf x-\mathbf y\|} \ \text{d}y, @f] then $\mathbf g(\mathbf x) = -\nabla \varphi(\mathbf x)$. If we assume that @@ -464,12 +464,12 @@ the following quantities: partial differential equation in all of ${\mathbb R}^3$ exploiting the radial symmetry): @f[ - \mathbf g(\mathbf x) = + \mathbf g(\mathbf x) = \left\{ \begin{array}{ll} -\frac{4}{3}\pi G \rho \|\mathbf x\| \frac{\mathbf x}{\|\mathbf x\|} & \text{for} \ \|\mathbf x\|The second to last parameter we need to specify is the viscosity @@ -533,10 +533,10 @@ the following quantities: viscosity of the material that flows into the area vacated under the rebounding continental plates. - Using this technique, values around $\eta=10^{21} \text{Pa\; s} + Using this technique, values around $\eta=10^{21} \text{Pa\; s} = 10^{21} \frac{\text{N\; s}}{\text{m}^2} = 10^{21} \frac{\text{kg}}{\text{m\; s}}$ have been found as the most - likely, though the error bar on this is at least one order of magnitude. + likely, though the error bar on this is at least one order of magnitude. While we will use this value, we again have to caution that there are many physical reasons to assume that this is not the correct value. First, it @@ -563,7 +563,7 @@ the following quantities: conductivity, $\rho$ the density, and $c_p$ the specific heat. For this, the literature indicates that it increases from around $0.7$ in the upper mantle to around $1.7 \frac{\text{mm}^2}{\text{s}}$ in the lower - mantle, though the exact value + mantle, though the exact value is not really all that important: heat transport through convection is several orders of magnitude more important than through thermal conduction. It may be of interest to know that perovskite, the most abundant @@ -573,10 +573,14 @@ the following quantities: transport through radiative transfer is more efficient than through thermal conduction. - In view of these considerations, let us choose + In view of these considerations, let us choose $\kappa=1 \frac{\text{mm}^2}{\text{s}} =10^{-6} \frac{\text{m}^2}{\text{s}}$ for the purpose of this program. All of these pieces of equation data are defined in the program in the -EquationData namespace. +EquationData namespace. When run, the program produces +maximal velocities around 10 centimeters per year, approximately the +physically correct order of magnitude. We will set the end time to 10 +million years. + diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 2bb9e1d82d..8a94f3537d 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -104,7 +104,7 @@ namespace EquationData const double T1 = 700+273; const double year_in_seconds = 60*60*24*365.2425; - + const double end_time = 1e7 * year_in_seconds; template Tensor<1,dim> gravity_vector (const Point &p) @@ -2941,7 +2941,7 @@ void BoussinesqFlowProblem::run () old_old_temperature_solution = old_temperature_solution; old_temperature_solution = temperature_solution; } - while (time <= 1e8 * EquationData::year_in_seconds); + while (time <= EquationData::end_time); } -- 2.39.5