[cig-commits] commit 2022 by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Nov 20 14:47:09 PST 2013
Revision 2022
Add another test.
A trunk/aspect/tests/compressibility/
A trunk/aspect/tests/compressibility/screen-output
A trunk/aspect/tests/compressibility.cc
A trunk/aspect/tests/compressibility.prm
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2022&peg=2022
Diff:
Added: trunk/aspect/tests/compressibility/screen-output
===================================================================
--- trunk/aspect/tests/compressibility/screen-output (rev 0)
+++ trunk/aspect/tests/compressibility/screen-output 2013-11-20 22:46:40 UTC (rev 2022)
@@ -0,0 +1,22 @@
+Loading shared library <./libcompressibility.so>
+
+Running with 1 MPI task.
+Number of active cells: 1,024 (on 6 levels)
+Number of degrees of freedom: 13,764 (8,450+1,089+4,225)
+
+*** Timestep 0: t=0 seconds
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system... 26 iterations.
+ Nonlinear Stokes residual: 21.2064
+
+ Postprocessing:
+ Top/bottom flux: 1/2.718
+ Writing graphical output: output-compressibility/solution-00000
+
+Termination requested by criterion: end time
+
+
++---------------------------------------------+------------+------------+
++---------------------------------+-----------+------------+------------+
++---------------------------------+-----------+------------+------------+
+
Added: trunk/aspect/tests/compressibility.cc
===================================================================
--- trunk/aspect/tests/compressibility.cc (rev 0)
+++ trunk/aspect/tests/compressibility.cc 2013-11-20 22:46:40 UTC (rev 2022)
@@ -0,0 +1,212 @@
+#include <aspect/material_model/simple.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ using namespace dealii;
+
+ template <int dim>
+ class Compressibility : public MaterialModel::Simple<dim>
+ {
+ public:
+ /**
+ * @name Physical parameters used in the basic equations
+ * @{
+ */
+ virtual double density (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position) const;
+
+ virtual double compressibility (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const Point<dim> &position) const;
+
+ /**
+ * Return true if the compressibility() function returns something that
+ * is not zero.
+ */
+ virtual bool
+ is_compressible () const;
+ };
+
+ }
+}
+
+namespace aspect
+{
+ namespace MaterialModel
+ {
+
+ template <int dim>
+ double
+ Compressibility<dim>::
+ density (const double temperature,
+ const double pressure,
+ const std::vector<double> &composition,
+ const Point<dim> &position) const
+ {
+ return 1.0 + compressibility(temperature,pressure,composition,position) * pressure;
+ }
+
+ template <int dim>
+ double
+ Compressibility<dim>::
+ compressibility (const double ,
+ const double,
+ const std::vector<double> &,
+ const Point<dim> &) const
+ {
+ return 1.0;
+ }
+
+ template <int dim>
+ bool
+ Compressibility<dim>::
+ is_compressible () const
+ {
+ return true;
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ ASPECT_REGISTER_MATERIAL_MODEL(Compressibility,
+ "compressibility",
+ "A simple material model that is like the "
+ "'Simple' model, but has a non-zero compressibility.")
+ }
+}
+
+
+
+#include <aspect/postprocess/interface.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ template <int dim>
+ class Compressibility : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+ {
+ public:
+ virtual
+ std::pair<std::string,std::string>
+ execute (TableHandler &statistics);
+ };
+ }
+}
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_values.h>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ template <int dim>
+ std::pair<std::string,std::string>
+ Compressibility<dim>::execute (TableHandler &statistics)
+ {
+ // create a quadrature formula based on the temperature element alone.
+ // be defensive about determining that what we think is the temperature
+ // element, is it in fact
+ Assert (this->get_fe().n_base_elements() == 3+(this->n_compositional_fields()>0 ? 1 : 0),
+ ExcNotImplemented());
+ const QGauss<dim-1> quadrature_formula (this->get_fe().base_element(2).degree+1);
+
+ FEFaceValues<dim> fe_face_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature_formula,
+ update_gradients | update_values |
+ update_q_points | update_JxW_values);
+
+ std::vector<std::vector<double> > composition_values (this->n_compositional_fields(),std::vector<double> (quadrature_formula.size()));
+ std::vector<Tensor<1,dim> > velocity_values(quadrature_formula.size());
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->get_dof_handler().begin_active(),
+ endc = this->get_dof_handler().end();
+
+ typename MaterialModel::Interface<dim>::MaterialModelInputs in(fe_face_values.n_quadrature_points, this->n_compositional_fields());
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs out(fe_face_values.n_quadrature_points);
+
+ // compute the integral of the viscosity. since we're on a unit box,
+ // this also is the average value
+ double bottom_flux_integral = 0;
+ double top_flux_integral = 0;
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ for (unsigned int f=0;f<2*dim;++f)
+ if (cell->at_boundary(f)
+ &&
+ ((cell->face(f)->boundary_indicator() == 2)
+ ||
+ (cell->face(f)->boundary_indicator() == 3)))
+ {
+ fe_face_values.reinit (cell,f);
+ fe_face_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(),
+ in.temperature);
+ fe_face_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(),
+ in.pressure);
+ fe_face_values[this->introspection().extractors.velocities].get_function_symmetric_gradients (this->get_solution(),
+ in.strain_rate);
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ fe_face_values[this->introspection().extractors.compositional_fields[c]].get_function_values(this->get_solution(),
+ composition_values[c]);
+ fe_face_values[this->introspection().extractors.velocities].get_function_values (this->get_solution(),
+ velocity_values);
+
+ in.position = fe_face_values.get_quadrature_points();
+ for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
+ {
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ in.composition[i][c] = composition_values[c][i];
+ }
+
+ this->get_material_model().evaluate(in, out);
+
+ if(cell->face(f)->boundary_indicator() == 2)
+ for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
+ bottom_flux_integral += out.densities[q] * velocity_values[q][1] * fe_face_values.JxW(q);
+ if(cell->face(f)->boundary_indicator() == 3)
+ for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
+ top_flux_integral += out.densities[q] * velocity_values[q][1] * fe_face_values.JxW(q);
+ }
+
+ std::ostringstream screen_text1;
+ std::ostringstream screen_text2;
+ screen_text1.precision(4);
+ screen_text2.precision(4);
+ screen_text1 << Utilities::MPI::sum(bottom_flux_integral, this->get_mpi_communicator());
+ screen_text2 << Utilities::MPI::sum(top_flux_integral, this->get_mpi_communicator());
+
+ return std::pair<std::string, std::string> ("Top/bottom flux:",
+ screen_text2.str() + "/" + screen_text1.str());
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Postprocess
+ {
+ ASPECT_REGISTER_POSTPROCESSOR(Compressibility,
+ "compressibility",
+ "A postprocessor that computes some statistics about "
+ "the mass fluxes.")
+ }
+}
Added: trunk/aspect/tests/compressibility.prm
===================================================================
--- trunk/aspect/tests/compressibility.prm (rev 0)
+++ trunk/aspect/tests/compressibility.prm 2013-11-20 22:46:40 UTC (rev 2022)
@@ -0,0 +1,134 @@
+# A testcase that demonstrates compressibility. It solves the
+# equations on a box with Dirichlet boundary conditions equal to
+# (0,1) at the bottom boundary, a tangential velocity at the sides
+# and a free boundary at the top.
+#
+# We then have a constant compressibilty that should lead to an
+# exponential density profile. The mass flux (product of density
+# and velocity) into the box and out of the box should be the same.
+#
+# We can only tell that the correct mass flux is computed in a
+# postprocessor, which we implement in the .cc file.
+#
+#
+# The correct solution has an exponentially increasing pressure and
+# density (top to bottom) where p(z)=e^(1-z)-1 (depth=1-z) and
+# rho(z)=e^(1-z). To conserve mass, this then requires a velocity
+# in z-direction equal to
+# v(z) = e^(-(1-z))
+# However, this is the product of the equation
+# div(rho u) = 0
+# which we instead solve in each iteration as
+# div u = -1/rho drho/dp u^* . g
+# where u^* is a linearization point. In the first iteration, we start
+# with u^*=0, so we solve
+# div u = 0
+# and the velocity field is of course completely wrong (it does not
+# accelerate towards the top as expected but is in fact constant).
+#
+# This testcase therefore verifies that we compute this wrong solution
+# in the first iteration. We use the Stokes only/max_nonlinear_iterations=1
+# option for the solver scheme.
+
+set Dimension = 2
+set CFL number = 1.0
+set End time = 0
+set Output directory = output-compressibility
+set Start time = 0
+set Adiabatic surface temperature = 0
+set Surface pressure = 0
+set Use years in output instead of seconds = false # default: true
+
+set Nonlinear solver scheme = Stokes only
+set Max nonlinear iterations = 1
+
+
+set Additional shared libraries = ./libcompressibility.so
+
+
+subsection Boundary temperature model
+ set Model name = box
+end
+
+
+
+subsection Gravity model
+ set Model name = vertical
+end
+
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 1
+ set Y extent = 1
+ set Z extent = 1
+ end
+end
+
+
+# temperature field doesn't matter. set it to zero
+subsection Initial conditions
+ set Model name = function
+ subsection Function
+ set Function expression = 0
+ end
+end
+
+
+# no gravity. the pressure will equal just the dynamic component
+subsection Gravity model
+ set Model name = vertical
+ subsection Vertical
+ set Magnitude = 1.0
+ end
+end
+
+
+subsection Material model
+ set Model name = compressibility
+
+ subsection Simple model
+ set Reference density = 1 # default: 3300
+ set Reference specific heat = 1250
+ set Reference temperature = 0 # default: 293
+ set Thermal conductivity = 1e-6 # default: 4.7
+ set Thermal expansion coefficient = 0
+ set Viscosity = 1 # default: 5e24
+ end
+end
+
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 5
+end
+
+
+subsection Model settings
+ set Fixed temperature boundary indicators =
+ set Tangential velocity boundary indicators = 0, 1
+ set Zero velocity boundary indicators =
+ set Prescribed velocity boundary indicators = 2: function
+ set Include shear heating = false
+end
+
+subsection Boundary velocity model
+ subsection Function
+ set Variable names = x,y
+ set Function expression = 0;1
+ end
+end
+
+subsection Postprocess
+ set List of postprocessors = compressibility, visualization
+
+ subsection Visualization
+
+ set List of output variables = density
+ set Number of grouped files = 0
+ set Output format = vtu
+ set Time between graphical output = 0 # default: 1e8
+ end
+end
More information about the CIG-COMMITS
mailing list