[cig-commits] [commit] master: Implemented new pressure compatibility correction. (8d254c2)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Aug 4 08:06:24 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/78abb106b73747fb364d8481c9df27e32aed4597...ea45ef11705f505f2d2d63c77aa04f2d7be9a4c1
>---------------------------------------------------------------
commit 8d254c21b5348977a6d9a005c518806ba4486f46
Author: Rene Gassmoeller <R.Gassmoeller at mailbox.org>
Date: Sat Jun 7 13:01:05 2014 +0200
Implemented new pressure compatibility correction.
>---------------------------------------------------------------
8d254c21b5348977a6d9a005c518806ba4486f46
source/simulator/core.cc | 9 +-
source/simulator/helper_functions.cc | 59 ++++++++++---
source/simulator/solver.cc | 4 +-
tests/pressure-compatibility-2.prm | 99 ++++++++++++++++++++++
tests/pressure-compatibility-2/screen-output | 50 +++++++++++
tests/pressure-compatibility-2/statistics | 21 +++++
.../pressure-compatibility.prm | 44 ++++------
tests/pressure-compatibility/screen-output | 60 +++++++++++++
tests/pressure-compatibility/statistics | 23 +++++
9 files changed, 323 insertions(+), 46 deletions(-)
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 97af832..ca3c130 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -323,14 +323,11 @@ namespace aspect
p != parameters.tangential_velocity_boundary_indicators.end();
++p)
open_velocity_boundary_indicators.erase (*p);
-//TODO: The "correct" condition would be to not do the correction for compressible
- // models if there are either open boundaries, or if there are boundaries
- // where the prescribed velocity is so that in- and outflux do not exactly
- // match
+
+ // we need to do the rhs compatibility modification, if the model is
+ // compressible, and there is no open boundary to balance the pressure
do_pressure_rhs_compatibility_modification = (material_model->is_compressible()
&&
- (parameters.prescribed_velocity_boundary_indicators.size() == 0)
- &&
(open_velocity_boundary_indicators.size() == 0));
// make sure that we don't have to fill every column of the statistics
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index 1630bbb..c88d483 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -809,18 +809,57 @@ namespace aspect
if (parameters.use_locally_conservative_discretization)
AssertThrow(false, ExcNotImplemented());
- if (do_pressure_rhs_compatibility_modification)
- {
- // TODO: currently does not work if velocity and
- // pressure are in the same block.
- Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
- ExcNotImplemented());
+ // TODO: currently does not work if velocity and
+ // pressure are in the same block.
+ Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
+ ExcNotImplemented());
- const double mean = vector.block(introspection.block_indices.pressure).mean_value();
- const double correction = -mean*vector.block(introspection.block_indices.pressure).size()/global_volume;
+ const QGauss<dim-1> quadrature_formula (parameters.stokes_velocity_degree+1);
- vector.block(introspection.block_indices.pressure).add(correction, pressure_shape_function_integrals.block(introspection.block_indices.pressure));
- }
+ FEFaceValues<dim> fe_face_values (mapping,
+ finite_element,
+ quadrature_formula,
+ update_normal_vectors |
+ update_q_points |
+ update_JxW_values);
+
+ double local_normal_velocity_integral = 0;
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ // for every surface face that is part of a geometry boundary
+ // and that is owned by this processor,
+ // integrate the normal velocity magnitude.
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->face(f)->at_boundary() &&
+ parameters.prescribed_velocity_boundary_indicators.find(cell->face(f)->boundary_indicator())!=
+ parameters.prescribed_velocity_boundary_indicators.end())
+ {
+ fe_face_values.reinit (cell, f);
+
+ // for each of the quadrature points, evaluate the
+ // normal velocity and add it to the integral
+ for (unsigned int q=0; q<quadrature_formula.size(); ++q)
+ {
+ const Tensor<1,dim> velocity =
+ velocity_boundary_conditions.find(cell->face(f)->boundary_indicator())->second->boundary_velocity(fe_face_values.quadrature_point(q));
+
+ const double normal_velocity = velocity * fe_face_values.normal_vector(q);
+ local_normal_velocity_integral += normal_velocity * fe_face_values.JxW(q);
+ }
+ }
+
+ const double global_normal_velocity_integral =
+ Utilities::MPI::sum (local_normal_velocity_integral,mpi_communicator);
+
+ const double mean = vector.block(introspection.block_indices.pressure).mean_value();
+ const double correction = (local_normal_velocity_integral - mean * vector.block(introspection.block_indices.pressure).size()) / global_volume;
+
+ vector.block(introspection.block_indices.pressure).add(correction, pressure_shape_function_integrals.block(introspection.block_indices.pressure));
}
diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc
index 11a8197..23c2fcc 100644
--- a/source/simulator/solver.cc
+++ b/source/simulator/solver.cc
@@ -485,7 +485,7 @@ namespace aspect
}
distributed_stokes_solution.compress(VectorOperation::insert);
- if (material_model->is_compressible ())
+ if (do_pressure_rhs_compatibility_modification)
make_pressure_rhs_compatible(system_rhs);
// we need a temporary vector for the residual (even if we don't care about it)
@@ -582,7 +582,7 @@ namespace aspect
// if the model is compressible then we need to adjust the right hand
// side of the equation to make it compatible with the matrix on the
// left
- if (material_model->is_compressible ())
+ if (do_pressure_rhs_compatibility_modification)
make_pressure_rhs_compatible(system_rhs);
// (ab)use the distributed solution vector to temporarily put a residual in
diff --git a/tests/pressure-compatibility-2.prm b/tests/pressure-compatibility-2.prm
new file mode 100644
index 0000000..e181c92
--- /dev/null
+++ b/tests/pressure-compatibility-2.prm
@@ -0,0 +1,99 @@
+# This is a compressible test case that shall
+# test the behaviour of the pressure compatibility
+# mechanism in case of prescribed normal velocities
+# integrated over the boundaries. It is converging
+# in this setup, but if the solver tolerance is
+# decreased or the prescribed normal velocity is
+# increased it fails to converge, which is exactly
+# the behaviour we intend.
+
+set CFL number = 1
+set End time = 400
+set Adiabatic surface temperature = 0.5
+set Use years in output instead of seconds = false
+set Linear solver tolerance = 1e-7
+set Pressure normalization = no
+
+subsection Boundary temperature model
+ set Model name = box
+ subsection Box
+ set Top temperature = 0
+ set Bottom temperature = 2
+ end
+end
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 1
+ set Y extent = 1
+ end
+end
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 1
+ end
+end
+
+
+subsection Initial conditions
+ set Model name = harmonic perturbation
+ subsection Harmonic perturbation
+ set Magnitude = 0.01
+ end
+end
+
+subsection Material model
+ set Model name = simple compressible
+
+ subsection Simple compressible model
+ set Reference density = 1
+ set Reference specific heat = 1
+ set Thermal expansion coefficient = 1
+ set Viscosity = 1
+ set Thermal conductivity = 1e-7
+ set Reference compressibility = 1e-6
+ end
+
+end
+
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+
+ set Initial global refinement = 3
+
+ set Refinement fraction = 0.0
+ set Coarsening fraction = 0.0
+
+ set Strategy = velocity
+
+ set Time steps between mesh refinement = 0
+end
+
+
+subsection Model settings
+ set Fixed temperature boundary indicators = 2,3
+
+ set Include adiabatic heating = true
+ set Include shear heating = true
+
+ set Tangential velocity boundary indicators = 0,1,3
+ set Prescribed velocity boundary indicators = 2:function
+end
+
+subsection Boundary velocity model
+ subsection Function
+ set Variable names = x,y
+ set Function expression = 0;1e-9
+ end
+end
+
+
+subsection Postprocess
+ set List of postprocessors = velocity statistics, pressure statistics, temperature statistics
+end
diff --git a/tests/pressure-compatibility-2/screen-output b/tests/pressure-compatibility-2/screen-output
new file mode 100644
index 0000000..a673ade
--- /dev/null
+++ b/tests/pressure-compatibility-2/screen-output
@@ -0,0 +1,50 @@
+-----------------------------------------------------------------------------
+-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
+-- . version 1.2.pre
+-- . running in DEBUG mode
+-- . running with 1 MPI process
+-- . using Trilinos
+-----------------------------------------------------------------------------
+
+Number of active cells: 64 (on 4 levels)
+Number of degrees of freedom: 948 (578+81+289)
+
+*** Timestep 0: t=0 seconds
+ Solving temperature system... 0 iterations.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 17 iterations.
+
+ Postprocessing:
+ RMS, max velocity: 6.6e-05 m/s, 0.000142 m/s
+ Pressure min/avg/max: -0.002229 Pa, 0.5 Pa, 0.998 Pa
+ Temperature min/avg/max: 0 K, 0.8632 K, 2 K
+
+*** Timestep 1: t=400 seconds
+ Solving temperature system... 11 iterations.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 17 iterations.
+
+ Postprocessing:
+ RMS, max velocity: 6.66e-05 m/s, 0.000141 m/s
+ Pressure min/avg/max: -0.007536 Pa, 0.5 Pa, 0.9937 Pa
+ Temperature min/avg/max: 0 K, 0.864 K, 2 K
+
+Termination requested by criterion: end time
+
+
++---------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start | 0.74s | |
+| | | |
+| Section | no. calls | wall time | % of total |
++---------------------------------+-----------+------------+------------+
+| Assemble Stokes system | 2 | 0.187s | 25% |
+| Assemble temperature system | 2 | 0.149s | 20% |
+| Build Stokes preconditioner | 2 | 0.13s | 18% |
+| Build temperature preconditioner| 2 | 0.00236s | 0.32% |
+| Solve Stokes system | 2 | 0.0214s | 2.9% |
+| Solve temperature system | 2 | 0.00236s | 0.32% |
+| Initialization | 2 | 0.077s | 10% |
+| Postprocessing | 2 | 0.0268s | 3.6% |
+| Setup dof systems | 1 | 0.0577s | 7.8% |
++---------------------------------+-----------+------------+------------+
+
diff --git a/tests/pressure-compatibility-2/statistics b/tests/pressure-compatibility-2/statistics
new file mode 100644
index 0000000..4bd5258
--- /dev/null
+++ b/tests/pressure-compatibility-2/statistics
@@ -0,0 +1,21 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Number of mesh cells
+# 4: Number of Stokes degrees of freedom
+# 5: Number of temperature degrees of freedom
+# 6: Iterations for temperature solver
+# 7: Iterations for Stokes solver
+# 8: Velocity iterations in Stokes preconditioner
+# 9: Schur complement iterations in Stokes preconditioner
+# 10: Time step size (seconds)
+# 11: RMS velocity (m/s)
+# 12: Max. velocity (m/s)
+# 13: Minimal pressure (Pa)
+# 14: Average pressure (Pa)
+# 15: Maximal pressure (Pa)
+# 16: Minimal temperature (K)
+# 17: Average temperature (K)
+# 18: Maximal temperature (K)
+# 19: Average nondimensional temperature (K)
+0 0.0000e+00 64 659 289 0 17 18 18 4.0000e+02 6.60174323e-05 1.42344742e-04 -2.22938226e-03 5.00000168e-01 9.97983627e-01 0.00000000e+00 8.63191646e-01 2.00000000e+00 4.31595823e-01
+1 4.0000e+02 64 659 289 11 17 18 18 4.4360e+02 6.66435963e-05 1.40632995e-04 -7.53585560e-03 5.00000168e-01 9.93683210e-01 0.00000000e+00 8.64018267e-01 2.00000000e+00 4.32009133e-01
diff --git a/doc/manual/cookbooks/platelike-boundary/platelike.prm b/tests/pressure-compatibility.prm
similarity index 51%
copy from doc/manual/cookbooks/platelike-boundary/platelike.prm
copy to tests/pressure-compatibility.prm
index a275c2d..943324e 100644
--- a/doc/manual/cookbooks/platelike-boundary/platelike.prm
+++ b/tests/pressure-compatibility.prm
@@ -1,18 +1,17 @@
-############### Global parameters
+# A test, whether we apply the right pressure compatibility
+# correction in the case of prescribed boundary conditions.
+# This test is essentially similar to the platelike-boundary
+# cookbook but uses a strongly compressible material model.
+# Previously we did not correct the pressure RHS in this setup,
+# which prevented it from converging.
+
set Dimension = 2
set Start time = 0
-set End time = 20
+set End time = 0.1
set Use years in output instead of seconds = false
-set Output directory = output
-
+set Adiabatic surface temperature = 0
-############### Parameters describing the model
-# Let us here choose again a box domain of size 2x1
-# where we fix the temperature at the bottom and top,
-# allow free slip along the bottom, left and right,
-# and prescribe the velocity along the top using the
-# `function' description.
subsection Geometry model
set Model name = box
@@ -32,8 +31,6 @@ subsection Model settings
end
-# We then set the temperature to one at the bottom and zero
-# at the top:
subsection Boundary temperature model
set Model name = box
@@ -44,8 +41,6 @@ subsection Boundary temperature model
end
-# The velocity along the top boundary models a spreading
-# center that is moving left and right:
subsection Boundary velocity model
subsection Function
set Variable names = x,z,t
@@ -55,10 +50,6 @@ subsection Boundary velocity model
end
-# We then choose a vertical gravity model and describe the
-# initial temperature with a vertical gradient. The default
-# strength for gravity is one. The material model is the
-# same as before.
subsection Gravity model
set Model name = vertical
end
@@ -75,29 +66,26 @@ end
subsection Material model
- set Model name = simple
+ set Model name = simple compressible
- subsection Simple model
+ subsection Simple compressible model
set Thermal conductivity = 1e-6
set Thermal expansion coefficient = 1e-4
+ set Reference compressibility = 7e-5
set Viscosity = 1
end
end
-# The final part of this input file describes how many times the
-# mesh is refined and what to do with the solution once computed
subsection Mesh refinement
set Initial adaptive refinement = 0
- set Initial global refinement = 5
+ set Initial global refinement = 3
set Time steps between mesh refinement = 0
end
subsection Postprocess
- set List of postprocessors = visualization, temperature statistics, heat flux statistics
-
- subsection Visualization
- set Time between graphical output = 0.1
- end
+ set List of postprocessors = temperature statistics, heat flux statistics, velocity statistics
end
+
+
diff --git a/tests/pressure-compatibility/screen-output b/tests/pressure-compatibility/screen-output
new file mode 100644
index 0000000..dd8a341
--- /dev/null
+++ b/tests/pressure-compatibility/screen-output
@@ -0,0 +1,60 @@
+-----------------------------------------------------------------------------
+-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
+-- . version 1.2.pre
+-- . running in DEBUG mode
+-- . running with 1 MPI process
+-- . using Trilinos
+-----------------------------------------------------------------------------
+
+Number of active cells: 64 (on 4 levels)
+Number of degrees of freedom: 948 (578+81+289)
+
+*** Timestep 0: t=0 seconds
+ Solving temperature system... 0 iterations.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 19 iterations.
+
+ Postprocessing:
+ Temperature min/avg/max: 0 K, 0.5 K, 1 K
+ Heat fluxes through boundary parts: 3.58e-21 W, 8.411e-21 W, -2e-06 W, 2e-06 W
+ RMS, max velocity: 0.353 m/s, 1.02 m/s
+
+*** Timestep 1: t=0.0625 seconds
+ Solving temperature system... 10 iterations.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 15 iterations.
+
+ Postprocessing:
+ Temperature min/avg/max: 0 K, 0.4995 K, 1 K
+ Heat fluxes through boundary parts: 1.745e-08 W, 1.801e-08 W, -1.999e-06 W, 2.17e-06 W
+ RMS, max velocity: 0.352 m/s, 1.02 m/s
+
+*** Timestep 2: t=0.1 seconds
+ Solving temperature system... 10 iterations.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 19 iterations.
+
+ Postprocessing:
+ Temperature min/avg/max: 0 K, 0.4994 K, 1 K
+ Heat fluxes through boundary parts: 2.867e-08 W, 2.956e-08 W, -2.002e-06 W, 2.265e-06 W
+ RMS, max velocity: 0.352 m/s, 1.02 m/s
+
+Termination requested by criterion: end time
+
+
++---------------------------------------------+------------+------------+
+| Total wallclock time elapsed since start | 1.05s | |
+| | | |
+| Section | no. calls | wall time | % of total |
++---------------------------------+-----------+------------+------------+
+| Assemble Stokes system | 3 | 0.278s | 26% |
+| Assemble temperature system | 3 | 0.232s | 22% |
+| Build Stokes preconditioner | 3 | 0.195s | 19% |
+| Build temperature preconditioner| 3 | 0.00337s | 0.32% |
+| Solve Stokes system | 3 | 0.0465s | 4.4% |
+| Solve temperature system | 3 | 0.00354s | 0.34% |
+| Initialization | 2 | 0.0741s | 7% |
+| Postprocessing | 3 | 0.046s | 4.4% |
+| Setup dof systems | 1 | 0.0586s | 5.6% |
++---------------------------------+-----------+------------+------------+
+
diff --git a/tests/pressure-compatibility/statistics b/tests/pressure-compatibility/statistics
new file mode 100644
index 0000000..fdee379
--- /dev/null
+++ b/tests/pressure-compatibility/statistics
@@ -0,0 +1,23 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Number of mesh cells
+# 4: Number of Stokes degrees of freedom
+# 5: Number of temperature degrees of freedom
+# 6: Iterations for temperature solver
+# 7: Iterations for Stokes solver
+# 8: Velocity iterations in Stokes preconditioner
+# 9: Schur complement iterations in Stokes preconditioner
+# 10: Time step size (seconds)
+# 11: Minimal temperature (K)
+# 12: Average temperature (K)
+# 13: Maximal temperature (K)
+# 14: Average nondimensional temperature (K)
+# 15: Outward heat flux through boundary with indicator 0 (W)
+# 16: Outward heat flux through boundary with indicator 1 (W)
+# 17: Outward heat flux through boundary with indicator 2 (W)
+# 18: Outward heat flux through boundary with indicator 3 (W)
+# 19: RMS velocity (m/s)
+# 20: Max. velocity (m/s)
+0 0.0000e+00 64 659 289 0 19 20 20 6.2500e-02 0.00000000e+00 5.00000000e-01 1.00000000e+00 5.00000000e-01 3.58008376e-21 8.41071040e-21 -2.00000000e-06 2.00000000e-06 3.52685352e-01 1.02187233e+00
+1 6.2500e-02 64 659 289 10 15 16 16 3.7500e-02 0.00000000e+00 4.99525428e-01 1.00000000e+00 4.99525428e-01 1.74548475e-08 1.80063133e-08 -1.99882536e-06 2.16999665e-06 3.52353376e-01 1.02484047e+00
+2 1.0000e-01 64 659 289 10 19 20 20 6.2500e-02 0.00000000e+00 4.99387316e-01 1.00000000e+00 4.99387316e-01 2.86707882e-08 2.95642484e-08 -2.00163236e-06 2.26455734e-06 3.52495631e-01 1.02339473e+00
More information about the CIG-COMMITS
mailing list