[cig-commits] commit 1986 by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Sat Oct 19 20:30:18 PDT 2013
Revision 1986
Add
A trunk/aspect/tests/shear-thinning/
A trunk/aspect/tests/shear-thinning/screen-output
A trunk/aspect/tests/shear-thinning.cc
A trunk/aspect/tests/shear-thinning.prm
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1986&peg=1986
Diff:
Added: trunk/aspect/tests/shear-thinning/screen-output
===================================================================
--- trunk/aspect/tests/shear-thinning/screen-output (rev 0)
+++ trunk/aspect/tests/shear-thinning/screen-output 2013-10-20 03:29:54 UTC (rev 1986)
@@ -0,0 +1,21 @@
+Loading shared library <./libshear-thinning.so>
+
+Running with 1 MPI task.
+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... 26 iterations.
+
+ Postprocessing:
+ Average viscosity: 0.5858
+
+Termination requested by criterion: end time
+
+
++---------------------------------------------+------------+------------+
++---------------------------------+-----------+------------+------------+
++---------------------------------+-----------+------------+------------+
+
Added: trunk/aspect/tests/shear-thinning.cc
===================================================================
--- trunk/aspect/tests/shear-thinning.cc (rev 0)
+++ trunk/aspect/tests/shear-thinning.cc 2013-10-20 03:29:54 UTC (rev 1986)
@@ -0,0 +1,182 @@
+#include <aspect/material_model/simple.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ using namespace dealii;
+
+ template <int dim>
+ class ShearThinning : public MaterialModel::Simple<dim>
+ {
+ public:
+ /**
+ * @name Physical parameters used in the basic equations
+ * @{
+ */
+ virtual double viscosity (const double temperature,
+ const double pressure,
+ const std::vector<double> &compositional_fields,
+ const SymmetricTensor<2,dim> &strain_rate,
+ const Point<dim> &position) const;
+
+ /**
+ * Return true if the viscosity() function returns something that
+ * may depend on the variable identifies by the argument.
+ */
+ virtual bool
+ viscosity_depends_on (const NonlinearDependence::Dependence dependence) const;
+ };
+
+ }
+}
+
+namespace aspect
+{
+ namespace MaterialModel
+ {
+
+ template <int dim>
+ double
+ ShearThinning<dim>::
+ viscosity (const double ,
+ const double,
+ const std::vector<double> &,
+ const SymmetricTensor<2,dim> &strain_rate,
+ const Point<dim> &) const
+ {
+ return 1./(1+strain_rate.norm());
+ }
+
+ template <int dim>
+ bool
+ ShearThinning<dim>::
+ viscosity_depends_on (const NonlinearDependence::Dependence dependence) const
+ {
+ return ((dependence & NonlinearDependence::strain_rate) != NonlinearDependence::none);
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ ASPECT_REGISTER_MATERIAL_MODEL(ShearThinning,
+ "shear thinning",
+ "A simple material model that is like the "
+ "'Simple' model, but has a viscosity equal to 1/|strain rate|.")
+ }
+}
+
+
+
+#include <aspect/postprocess/interface.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ template <int dim>
+ class ShearThinning : 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>
+ ShearThinning<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> quadrature_formula (this->get_fe().base_element(2).degree+1);
+
+ FEValues<dim> fe_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()));
+
+ 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_values.n_quadrature_points, this->n_compositional_fields());
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs out(fe_values.n_quadrature_points);
+
+ // compute the integral of the viscosity. since we're on a unit box,
+ // this also is the average value
+ double viscosity_integral = 0;
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit (cell);
+ fe_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(),
+ in.temperature);
+ fe_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(),
+ in.pressure);
+ fe_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_values[this->introspection().extractors.compositional_fields[c]].get_function_values(this->get_solution(),
+ composition_values[c]);
+
+ in.position = fe_values.get_quadrature_points();
+ for (unsigned int i=0; i<fe_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);
+
+
+ for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
+ viscosity_integral += out.viscosities[q] * fe_values.JxW(q);
+ }
+
+ std::ostringstream screen_text;
+ screen_text.precision(4);
+ screen_text << Utilities::MPI::sum(viscosity_integral, this->get_mpi_communicator());
+
+ return std::pair<std::string, std::string> ("Average viscosity:",
+ screen_text.str());
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Postprocess
+ {
+ ASPECT_REGISTER_POSTPROCESSOR(ShearThinning,
+ "shear thinning",
+ "A postprocessor that computes some statistics about "
+ "the viscosity.")
+ }
+}
Added: trunk/aspect/tests/shear-thinning.prm
===================================================================
--- trunk/aspect/tests/shear-thinning.prm (rev 0)
+++ trunk/aspect/tests/shear-thinning.prm 2013-10-20 03:29:54 UTC (rev 1986)
@@ -0,0 +1,107 @@
+# A testcase that demonstrates shear thinning. It solves the
+# equations on a box with Dirichlet boundary conditions equal to
+# (z,0), which then is also the velocity everywhere. This yields
+# a constant strain rate [[0,1/2],[1/2,0]] with norm |dot eps|=1/sqrt(2).
+#
+# We then have a viscosity that depends on
+# the strain rate as eta=1/(1+|dot eps|). Because the strain rate
+# is constant, so is the viscosity.
+#
+# In this version of the testcase, we only run a single IMPES step,
+# so the viscosity we use for the first velocity step is eta=1.
+# We can only tell that the correct viscosity is computed in a
+# postprocessor, which we implement in the .cc file. The correct
+# value that needs to be computed is viscosity=1/(1+|dot eps|),
+# i.e., viscosity=0.585786
+
+set Dimension = 2
+set CFL number = 1.0
+set End time = 0
+set Output directory = output-shear-thinning
+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 = IMPES
+
+set Additional shared libraries = ./libshear-thinning.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 = x
+ end
+end
+
+
+# no gravity. the pressure will equal just the dynamic component
+subsection Gravity model
+ set Model name = vertical
+ subsection Vertical
+ set Magnitude = 0
+ end
+end
+
+
+subsection Material model
+ set Model name = shear thinning
+
+ subsection Simple model
+ set Reference density = 1 # default: 3300
+ set Reference specific heat = 1250
+ set Reference temperature = 1 # default: 293
+ set Thermal conductivity = 1e-6 # default: 4.7
+ set Thermal expansion coefficient = 2e-5
+ set Viscosity = 1 # default: 5e24
+ end
+end
+
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 3
+end
+
+
+subsection Model settings
+ set Fixed temperature boundary indicators = 0, 1, 2, 3
+ set Tangential velocity boundary indicators =
+ set Zero velocity boundary indicators =
+ set Prescribed velocity boundary indicators = 0: function, 1: function, 2: function, 3: function
+end
+
+subsection Boundary velocity model
+ subsection Function
+ set Variable names = x,z
+ set Function expression = z;0
+ end
+end
+
+subsection Postprocess
+ set List of postprocessors = shear thinning
+end
+
More information about the CIG-COMMITS
mailing list