From 18cc77afa9bb6672805f20bcfe0d1f003c18a6d1 Mon Sep 17 00:00:00 2001 From: cedrict Date: Fri, 12 May 2017 17:58:54 -0400 Subject: [PATCH] tests for annulus & hollow_sphere benchmarks --- tests/annulus.cc | 562 ++++++++++++++++++++++++ tests/annulus.prm | 108 +++++ tests/annulus/screen-output | 44 ++ tests/hollow_sphere.cc | 696 ++++++++++++++++++++++++++++++ tests/hollow_sphere.prm | 99 +++++ tests/hollow_sphere/screen-output | 44 ++ 6 files changed, 1553 insertions(+) create mode 100644 tests/annulus.cc create mode 100644 tests/annulus.prm create mode 100644 tests/annulus/screen-output create mode 100644 tests/hollow_sphere.cc create mode 100644 tests/hollow_sphere.prm create mode 100644 tests/hollow_sphere/screen-output diff --git a/tests/annulus.cc b/tests/annulus.cc new file mode 100644 index 00000000000..5b9bcbbfeac --- /dev/null +++ b/tests/annulus.cc @@ -0,0 +1,562 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace aspect +{ + /** + * This is the "annulus" benchmark defined in the following paper: + * @code + * @Article{thph17, + * author = {C. Thieulot and E.G.P. Puckett and H. Lokavarapu}, + * title = {Stokes flow in an annulus: analytical solution and numerical benchmark}, + * journal = {xxx}, + * year = 2017, + * volume = xxx, + * number = {x}, + * publisher = {xxx}, + * pages = {xxx--xxx}} + * @endcode + * + */ + namespace AnnulusBenchmark + { + using namespace dealii; + + namespace AnalyticSolutions + { + Tensor<1,2> + Annulus_velocity (const Point<2> &pos, + const double k) + { + const double x = pos[0]; + const double y = pos[1]; + const double r = sqrt(x*x+y*y); + const double theta = atan2(y,x); + const double A=2.0, B=-3/log(2), C=-1; + const double f_r = A*r + B/r; + const double g_r = A*r/2 + B*log(r)/r + C/r; + const double v_r = g_r*k*sin(k*theta); + const double v_theta = f_r*cos(k*theta); + const double v_x = cos(theta)*v_r - sin(theta)*v_theta; + const double v_y = sin(theta)*v_r + cos(theta)*v_theta; + return Point<2> (v_x,v_y); + } + + double + Annulus_pressure (const Point<2> &pos, + const double k) + { + const double x = pos[0]; + const double y = pos[1]; + const double r=sqrt(x*x+y*y); + const double theta = atan2(y,x); + const double A=2.0, B=-3/log(2), C=-1; + const double f_r = 2*r + B/r; + const double g_r = A*r/2 + B*log(r)/r + C/r; + const double h_r=(2*g_r-f_r)/r; + const double l_r=0; + return k*h_r*sin(k*theta)+l_r; + } + + + /** + * The exact solution for the Annulus benchmark. + */ + template + class FunctionAnnulus : public Function + { + public: + FunctionAnnulus (const double beta) + : + Function(dim+2), + beta_(beta) + {} + + virtual void vector_value (const Point< dim > &pos, + Vector< double > &values) const + { + Assert (dim == 3, ExcNotImplemented()); + Assert (values.size() >= 4, ExcInternalError()); + + const Point<2> p (pos[0], pos[1]); + + const Tensor<1,2> v = AnalyticSolutions::Annulus_velocity (p, beta_); + values[0] = v[0]; + values[1] = v[1]; + values[2] = AnalyticSolutions::Annulus_pressure (p, beta_); + } + + private: + const double beta_; + }; + } + + + + template + class AnnulusBoundary : public BoundaryVelocity::Interface, public aspect::SimulatorAccess + { + public: + /** + * Constructor. + */ + AnnulusBoundary(); + + /** + * Return the boundary velocity as a function of position. + */ + virtual + Tensor<1,dim> + boundary_velocity (const types::boundary_id , + const Point &position) const; + + private: + const double beta; + }; + + + + + + + + /** + * @note This benchmark only talks about the flow field, not about a + * temperature field. All quantities related to the temperature are + * therefore set to zero in the implementation of this class. + _r* + * @ingroup MaterialModels + */ + template + class AnnulusMaterial : public MaterialModel::InterfaceCompatibility + { + public: + /** + * @name Physical parameters used in the basic equations + * @{ + */ + virtual double viscosity (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const SymmetricTensor<2,dim> &strain_rate, + const Point &position) const; + + virtual double density (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double compressibility (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double specific_heat (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double thermal_expansion_coefficient (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double thermal_conductivity (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + /** + * @} + */ + + /** + * @name Qualitative properties one can ask a material model + * @{ + */ + + /** + * Return whether the model is compressible or not. + * Incompressibility does not necessarily imply that the density is + * constant; rather, it may still depend on temperature or pressure. + * In the current context, compressibility means whether we should + * solve the contuity equation as $\nabla \cdot (\rho \mathbf u)=0$ + * (compressible Stokes) or as $\nabla \cdot \mathbf{u}=0$ + * (incompressible Stokes). + */ + virtual bool is_compressible () const; + /** + * @} + */ + /** + * Declare the parameters this class takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm); + + + + /** + * @name Reference quantities + * @{ + */ + virtual double reference_viscosity () const; + /** + * @} + */ + /** + * Returns the viscosity value in the inclusion + */ + double get_beta() const; + + /** + * viscosity value in the inclusion + */ + double beta; + }; + + template + double + AnnulusMaterial:: + viscosity (const double, + const double, + const std::vector &, /*composition*/ + const SymmetricTensor<2,dim> &, + const Point &p) const + { + return 1.; + } + + + template + double + AnnulusMaterial:: + reference_viscosity () const + { + return 1.; + } + + + template + double + AnnulusMaterial:: + specific_heat (const double, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0; + } + + + template + double + AnnulusMaterial:: + thermal_conductivity (const double, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0; + } + + + template + double + AnnulusMaterial:: + density (const double, + const double, + const std::vector &, /*composition*/ + const Point &pos) const + { + const double k=beta; + const double x = pos[0]; + const double y = pos[1]; + const double r=sqrt(x*x+y*y); + const double theta = atan2(y,x); + const double A=2.0, B=-3/log(2), C=-1; + const double f = A*r + B/r; + const double f_prime = 2 - B/std::pow(r,2.0); + const double g = A*r/2 + B*log(r)/r + C/r; + const double g_prime = A/2 - B*log(r)/std::pow(r,2.0) + B/std::pow(r,2.0) - C/std::pow(r,2.0); + const double g_prime_prime = -B/std::pow(r,3)*(3-2*log(r)) - 2./std::pow(r,3); + const double N = g_prime_prime - g_prime/r - (g*(std::pow(k,2) - 1))/std::pow(r,2.0) + f/std::pow(r,2.0) + f_prime/r; + const double density = N*k*sin(k*theta); + return density; + } + + + template + double + AnnulusMaterial:: + thermal_expansion_coefficient (const double temperature, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0; + } + + + template + double + AnnulusMaterial:: + compressibility (const double, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0.0; + } + + + template + bool + AnnulusMaterial:: + is_compressible () const + { + return false; + } + + template + void + AnnulusMaterial::declare_parameters (ParameterHandler &prm) + { + //create a global section in the parameter file for parameters + //that describe this benchmark. note that we declare them here + //in the material model, but other kinds of plugins (e.g., the gravity + //model below) may also read these parameters even though they do not + //declare them + prm.enter_subsection("Annulus benchmark"); + { + prm.declare_entry("Viscosity parameter", "0", + Patterns::Double (0), + "Viscosity parameter k in the Annulus benchmark."); + } + prm.leave_subsection(); + } + + + template + void + AnnulusMaterial::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Annulus benchmark"); + { + beta = prm.get_double ("Viscosity parameter"); + } + prm.leave_subsection(); + + // Declare dependencies on solution variables + this->model_dependence.viscosity = MaterialModel::NonlinearDependence::none; + this->model_dependence.density = MaterialModel::NonlinearDependence::none; + this->model_dependence.compressibility = MaterialModel::NonlinearDependence::none; + this->model_dependence.specific_heat = MaterialModel::NonlinearDependence::none; + this->model_dependence.thermal_conductivity = MaterialModel::NonlinearDependence::none; + } + + + template + double + AnnulusMaterial::get_beta() const + { + return beta; + } + + + + + template + AnnulusBoundary::AnnulusBoundary () + : + beta (0) + {} + + + + template <> + Tensor<1,2> + AnnulusBoundary<2>:: + boundary_velocity (const types::boundary_id , + const Point<2> &p) const + { + + const AnnulusMaterial<2> * + material_model + = dynamic_cast *>(&this->get_material_model()); + + return AnalyticSolutions::Annulus_velocity (p, material_model->get_beta()); + } + + + + template <> + Tensor<1,3> + AnnulusBoundary<3>:: + boundary_velocity (const types::boundary_id , + const Point<3> &p) const + { + Assert (false, ExcNotImplemented()); + return Tensor<1,3>(); + } + + + + + /** + *gravity model for the Annulus benchmark + */ + + template + class AnnulusGravity : public aspect::GravityModel::Interface + { + public: + virtual Tensor<1,dim> gravity_vector (const Point &pos) const; + + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm); + + double beta; + }; + + + /** + * A postprocessor that evaluates the accuracy of the solution. + * + * The implementation of error evaluators that correspond to the + * benchmarks defined in the paper Duretz et al. reference above. + */ + template + class AnnulusPostprocessor : public Postprocess::Interface, public ::aspect::SimulatorAccess + { + public: + /** + * Generate graphical output from the current solution. + */ + virtual + std::pair + execute (TableHandler &statistics); + }; + + template + std::pair + AnnulusPostprocessor::execute (TableHandler &statistics) + { + std_cxx1x::shared_ptr > ref_func; + { + const AnnulusMaterial * + material_model + = dynamic_cast *>(&this->get_material_model()); + + ref_func.reset (new AnalyticSolutions::FunctionAnnulus(material_model->get_beta())); + } + + const QGauss quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.velocities).degree+2); + + Vector cellwise_errors_u (this->get_triangulation().n_active_cells()); + Vector cellwise_errors_p (this->get_triangulation().n_active_cells()); + Vector cellwise_errors_ul2 (this->get_triangulation().n_active_cells()); + Vector cellwise_errors_pl2 (this->get_triangulation().n_active_cells()); + + ComponentSelectFunction comp_u(std::pair(0,dim), + dim+2); + ComponentSelectFunction comp_p(dim, dim+2); + + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_u, + quadrature_formula, + VectorTools::L1_norm, + &comp_u); + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_p, + quadrature_formula, + VectorTools::L1_norm, + &comp_p); + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_ul2, + quadrature_formula, + VectorTools::L2_norm, + &comp_u); + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_pl2, + quadrature_formula, + VectorTools::L2_norm, + &comp_p); + + const double u_l1 = Utilities::MPI::sum(cellwise_errors_u.l1_norm(),this->get_mpi_communicator()); + const double p_l1 = Utilities::MPI::sum(cellwise_errors_p.l1_norm(),this->get_mpi_communicator()); + const double u_l2 = std::sqrt(Utilities::MPI::sum(cellwise_errors_ul2.norm_sqr(),this->get_mpi_communicator())); + const double p_l2 = std::sqrt(Utilities::MPI::sum(cellwise_errors_pl2.norm_sqr(),this->get_mpi_communicator())); + + std::ostringstream os; + + os << std::scientific << u_l1 + << ", " << p_l1 + << ", " << u_l2 + << ", " << p_l2; + + return std::make_pair("Errors u_L1, p_L1, u_L2, p_L2:", os.str()); + } + + } +} + + + +// explicit instantiations +namespace aspect +{ + namespace AnnulusBenchmark + { + ASPECT_REGISTER_MATERIAL_MODEL(AnnulusMaterial, + "AnnulusMaterial", + "A material model that corresponds to the `Annulus' benchmark. " + "See the manual for more information.") + + ASPECT_REGISTER_BOUNDARY_VELOCITY_MODEL(AnnulusBoundary, + "AnnulusBoundary", + "Implementation of the velocity boundary conditions for the " + "`Annulus' benchmark. See the manual for more information about this " + "benchmark.") + + ASPECT_REGISTER_POSTPROCESSOR(AnnulusPostprocessor, + "AnnulusPostprocessor", + "A postprocessor that compares the solution of the `Annulus' benchmark " + "with the one computed by ASPECT " + "and reports the error. See the manual for more information.") + } +} + diff --git a/tests/annulus.prm b/tests/annulus.prm new file mode 100644 index 00000000000..1a1f37f12eb --- /dev/null +++ b/tests/annulus.prm @@ -0,0 +1,108 @@ +# The 2D Stokes in an annulus benchmark, for which an +# analytical solution is available. +# +# See the manual for more information. + +############### Global parameters +# We use a 2d setup. Since we are only interested +# in a steady state solution, we set the end time +# equal to the start time to force a single time +# step before the program terminates. + +set Additional shared libraries = ./libannulus.so +set Linear solver tolerance = 1e-12 +set Dimension = 2 +set Start time = 0 +set End time = 0 +set Use years in output instead of seconds = false +set Nonlinear solver scheme = Stokes only +set Output directory = output +set Pressure normalization = surface +set Number of cheap Stokes solver steps = 0 + +############### Parameters describing the model +# Because the temperature plays no role in this model we need not +# bother to describe temperature boundary conditions or +# the material parameters that pertain to the temperature. + +subsection Geometry model + set Model name = spherical shell + subsection Spherical shell + set Inner radius = 1 + set Outer radius = 2 + set Opening angle = 360 + end +end + +#Boundary conditions +subsection Model settings + set Tangential velocity boundary indicators = + set Prescribed velocity boundary indicators = 0 : AnnulusBoundary, \ + 1 : AnnulusBoundary +end + +#subsection Discretization +#set Stokes velocity polynomial degree = 1 +#set Use locally conservative discretization = true +#end + +subsection Material model + set Model name = AnnulusMaterial +end + +subsection Gravity model + set Model name = radial constant + subsection Radial constant + set Magnitude = 1 + end +end + +#Viscosity parameter is k, which controls the number +#of mantle convection cells + +subsection Annulus benchmark + set Viscosity parameter = 4 +end + +############### Parameters describing the temperature field +# As above, there is no need to set anything for the +# temperature boundary conditions. + +subsection Boundary temperature model + set Model name = box +end + +subsection Initial temperature model + set Model name = function + + subsection Function + set Function expression = 0 + end +end + + +############### Parameters describing the discretization +# The following parameters describe how often we want to refine +# the mesh globally and adaptively, what fraction of cells should +# be refined in each adaptive refinement step, and what refinement +# indicator to use when refining the mesh adaptively. + +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 2 + set Refinement fraction = 0.2 + set Strategy = velocity +end + +############### Parameters describing what to do with the solution +# The final section allows us to choose which postprocessors to +# run at the end of each time step. + +subsection Postprocess + set List of postprocessors = visualization, velocity statistics, AnnulusPostprocessor, depth average, memory statistics + + subsection Visualization + set List of output variables = density, viscosity, strain rate + end +end + diff --git a/tests/annulus/screen-output b/tests/annulus/screen-output new file mode 100644 index 00000000000..4ba79f7a9b7 --- /dev/null +++ b/tests/annulus/screen-output @@ -0,0 +1,44 @@ +----------------------------------------------------------------------------- +-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion. +-- . version 2.0.0-pre +-- . running in OPTIMIZED mode +-- . running with 1 MPI process +-- . using Trilinos +----------------------------------------------------------------------------- + +Loading shared library <./libannulus.so> + +Number of active cells: 192 (on 3 levels) +Number of degrees of freedom: 2,832 (1,728+240+864) + +*** Timestep 0: t=0 seconds + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+23 iterations. + Relative Stokes residual after nonlinear iteration 1: 1 + Solving Stokes system... 0+0 iterations. + Relative Stokes residual after nonlinear iteration 2: 2.09534e-13 + + Postprocessing: + Writing graphical output: output/solution/solution-00000 + RMS, max velocity: 1.08 m/s, 2.16 m/s + Errors u_L1, p_L1, u_L2, p_L2: 9.544368e-02, 6.064711e-01, 2.641477e-02, 3.399588e-01 + Writing depth average: output/depth_average.gnuplot + System matrix memory consumption: 1.00 MB + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 0.207s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assemble Stokes system | 2 | 0.0154s | 7.4% | +| Build Stokes preconditioner | 1 | 0.0127s | 6.2% | +| Solve Stokes system | 2 | 0.0339s | 16% | +| Initialization | 1 | 0.0341s | 16% | +| Postprocessing | 1 | 0.09s | 43% | +| Setup dof systems | 1 | 0.0101s | 4.9% | +| Setup initial conditions | 1 | 0.00497s | 2.4% | ++---------------------------------+-----------+------------+------------+ + diff --git a/tests/hollow_sphere.cc b/tests/hollow_sphere.cc new file mode 100644 index 00000000000..40bad3b608e --- /dev/null +++ b/tests/hollow_sphere.cc @@ -0,0 +1,696 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + + +namespace aspect +{ + /** + * This is the "HollowSphere" benchmark defined in the following paper: + * Analytical solution for viscous incompressible Stokes flow in a spherical shell + * C. Thieulot, in prep. + */ + namespace HollowSphereBenchmark + { + using namespace dealii; + + namespace AnalyticSolutions + { + Tensor<1,3> + hollow_sphere_velocity (const Point<3> &pos, + const double mmm) + { + const double gammma = 1.0; + + const std_cxx11::array spos = + aspect::Utilities::Coordinates::cartesian_to_spherical_coordinates(pos); + + const double r=spos[0]; + const double phi=spos[1]; + const double theta=spos[2]; + + const double R1 = 0.5; + const double R2 = 1.0; + + double alpha,beta,fr,gr; + + if (mmm == -1) + { + alpha=-gammma*(pow(R2,3)-pow(R1,3))/(pow(R2,3)*log(R1)-pow(R1,3)*log(R2)); + beta=-3*gammma*(log(R2)-log(R1))/(pow(R1,3)*log(R2)-pow(R2,3)*log(R1)) ; + fr=alpha/(r*r)+beta*r; + gr=-2/(r*r)*(alpha*log(r)+beta/3*pow(r,3)+gammma); + } + else + { + alpha=gammma*(mmm+1)*(pow(R1,-3)-pow(R2,-3))/(pow(R1,-mmm-4)-pow(R2,-mmm-4)); + beta=-3*gammma*(pow(R1,mmm+1)-pow(R2,mmm+1))/(pow(R1,mmm+4)-pow(R2,mmm+4)); + fr=alpha/pow(r,mmm+3)+beta*r; + gr=-2/(r*r)*(-alpha/(mmm+1)*pow(r,-mmm-1)+beta/3*pow(r,3)+gammma); + } + + const double v_r =gr*cos(theta); + const double v_theta=fr*sin(theta); + const double v_phi =fr*sin(theta); + const double v_x=sin(theta)*cos(phi)*v_r + cos(theta)*cos(phi)*v_theta-sin(phi)*v_phi; + const double v_y=sin(theta)*sin(phi)*v_r + cos(theta)*sin(phi)*v_theta+cos(phi)*v_phi; + const double v_z=cos(theta)*v_r - sin(theta)*v_theta; + + // create a Point<3> (because it has a constructor that takes + // three doubles) and return it (it automatically converts to + // the necessary Tensor<1,3>). + return Point<3> (v_x,v_y,v_z); + } + + double + hollow_sphere_pressure (const Point<3> &pos, + const double mmm) + { + const double gammma = 1.0; + const double mu0=1; + + const std_cxx11::array spos = + aspect::Utilities::Coordinates::cartesian_to_spherical_coordinates(pos); + + const double r=spos[0]; + const double phi=spos[1]; + const double theta=spos[2]; + + const double R1 = 0.5; + const double R2 = 1.0; + + double alpha,beta,fr,gr,hr,mur; + + + if (mmm == -1) + { + mur=mu0; + alpha=-gammma*(pow(R2,3)-pow(R1,3))/(pow(R2,3)*log(R1)-pow(R1,3)*log(R2)); + beta=-3*gammma*(log(R2)-log(R1))/(pow(R1,3)*log(R2)-pow(R2,3)*log(R1)) ; + fr=alpha/(r*r)+beta*r; + gr=-2/(r*r)*(alpha*log(r)+beta/3*pow(r,3)+gammma); + hr=2/r*gr*mur; + } + else + { + mur=mu0*pow(r,mmm+1); + alpha=gammma*(mmm+1)*(pow(R1,-3)-pow(R2,-3))/(pow(R1,-mmm-4)-pow(R2,-mmm-4)); + beta=-3*gammma*(pow(R1,mmm+1)-pow(R2,mmm+1))/(pow(R1,mmm+4)-pow(R2,mmm+4)); + fr=alpha/pow(r,mmm+3)+beta*r; + gr=-2/(r*r)*(-alpha/(mmm+1)*pow(r,-mmm-1)+beta/3*pow(r,3)+gammma); + hr=(mmm+3)/r*gr*mur; + } + + return hr*cos(theta); + } + + + /** + * The exact solution for the HollowSphere benchmark. + */ + template + class FunctionHollowSphere : public Function + { + public: + FunctionHollowSphere (const double mmm) + : + Function(dim+2), + mmm_(mmm) + {} + + virtual void vector_value (const Point< dim > &pos, + Vector< double > &values) const + { + Assert (dim == 3, ExcNotImplemented()); + Assert (values.size() >= 4, ExcInternalError()); + + const Point<3> p (pos[0], pos[1], pos[2]); + + const Tensor<1,3> v = AnalyticSolutions::hollow_sphere_velocity (p, mmm_); + values[0] = v[0]; + values[1] = v[1]; + values[2] = v[2]; + + values[3] = AnalyticSolutions::hollow_sphere_pressure (p, mmm_); + } + + private: + const double mmm_; + }; + } + + + + template + class HollowSphereBoundary : public BoundaryVelocity::Interface + { + public: + /** + * Constructor. + */ + HollowSphereBoundary(); + + /** + * Return the boundary velocity as a function of position. + */ + virtual + Tensor<1,dim> + boundary_velocity (const types::boundary_id , + const Point &position) const; + + + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm); + + + + private: + double mmm; + }; + + template + HollowSphereBoundary::HollowSphereBoundary () + : + mmm (0) + {} + + + + template <> + Tensor<1,2> + HollowSphereBoundary<2>:: + boundary_velocity (const types::boundary_id , + const Point<2> &p) const + { + Assert (false, ExcNotImplemented()); + return Tensor<1,2>(); + } + + + + template <> + Tensor<1,3> + HollowSphereBoundary<3>:: + boundary_velocity (const types::boundary_id , + const Point<3> &p) const + { + return AnalyticSolutions::hollow_sphere_velocity (p, mmm); + } + + + + + /** + * @note This benchmark only talks about the flow field, not about a + * temperature field. All quantities related to the temperature are + * therefore set to zero in the implementation of this class. + * + * @ingroup MaterialModels + */ + template + class HollowSphereMaterial : public MaterialModel::InterfaceCompatibility + { + public: + /** + * @name Physical parameters used in the basic equations + * @{ + */ + virtual double viscosity (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const SymmetricTensor<2,dim> &strain_rate, + const Point &position) const; + + virtual double density (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double compressibility (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double specific_heat (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double thermal_expansion_coefficient (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + + virtual double thermal_conductivity (const double temperature, + const double pressure, + const std::vector &compositional_fields, + const Point &position) const; + /** + * @} + */ + + /** + * @name Qualitative properties one can ask a material model + * @{ + */ + + /** + * Return whether the model is compressible or not. + * Incompressibility does not necessarily imply that the density is + * constant; rather, it may still depend on temperature or pressure. + * In the current context, compressibility means whether we should + * solve the contuity equation as $\nabla \cdot (\rho \mathbf u)=0$ + * (compressible Stokes) or as $\nabla \cdot \mathbf{u}=0$ + * (incompressible Stokes). + */ + virtual bool is_compressible () const; + /** + * @} + */ + /** + * Declare the parameters this class takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm); + + + + /** + * @name Reference quantities + * @{ + */ + virtual double reference_viscosity () const; + /** + * @} + */ + /** + * Returns the viscosity value in the inclusion + */ + double get_mmm() const; + + /** + * viscosity value in the inclusion + */ + double mmm; + }; + + template + double + HollowSphereMaterial:: + viscosity (const double, + const double, + const std::vector &, /*composition*/ + const SymmetricTensor<2,dim> &, + const Point &pos) const + { + const std_cxx11::array spos = + aspect::Utilities::Coordinates::cartesian_to_spherical_coordinates(pos); + + const double r=spos[0]; + const double mu = pow(r,mmm+1); + return mu; + } + + + template + double + HollowSphereMaterial:: + reference_viscosity () const + { + return 1.; + } + + + template + double + HollowSphereMaterial:: + specific_heat (const double, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0; + } + + + template + double + HollowSphereMaterial:: + thermal_conductivity (const double, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0; + } + + + template + double + HollowSphereMaterial:: + density (const double, + const double, + const std::vector &, /*composition*/ + const Point &p) const + { + return 1; + } + + + template + double + HollowSphereMaterial:: + thermal_expansion_coefficient (const double temperature, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0; + } + + + template + double + HollowSphereMaterial:: + compressibility (const double, + const double, + const std::vector &, /*composition*/ + const Point &) const + { + return 0.0; + } + + + template + bool + HollowSphereMaterial:: + is_compressible () const + { + return false; + } + + template + void + HollowSphereMaterial::declare_parameters (ParameterHandler &prm) + { + //create a global section in the parameter file for parameters + //that describe this benchmark. note that we declare them here + //in the material model, but other kinds of plugins (e.g., the gravity + //model below) may also read these parameters even though they do not + //declare them + prm.enter_subsection("HollowSphere benchmark"); + { + prm.declare_entry("Viscosity parameter", "-1", + Patterns::Double (), + "Viscosity in the HollowSphere benchmark."); + } + prm.leave_subsection(); + } + + + template + void + HollowSphereMaterial::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("HollowSphere benchmark"); + { + mmm = prm.get_double ("Viscosity parameter"); + } + prm.leave_subsection(); + + // Declare dependencies on solution variables + this->model_dependence.viscosity = MaterialModel::NonlinearDependence::none; + this->model_dependence.density = MaterialModel::NonlinearDependence::none; + this->model_dependence.compressibility = MaterialModel::NonlinearDependence::none; + this->model_dependence.specific_heat = MaterialModel::NonlinearDependence::none; + this->model_dependence.thermal_conductivity = MaterialModel::NonlinearDependence::none; + } + + template + void + HollowSphereBoundary::declare_parameters (ParameterHandler &prm) + { + //nothing to declare here. This plugin will however, read parameters + //declared by the material model in the "HollowSphere benchmark" section + } + + + template + void + HollowSphereBoundary::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("HollowSphere benchmark"); + { + mmm = prm.get_double ("Viscosity parameter"); + } + prm.leave_subsection(); + } + + + + template + double + HollowSphereMaterial::get_mmm() const + { + return mmm; + } + + /** + *gravity model for the HollowSphere benchmark + */ + + template + class HollowSphereGravity : public aspect::GravityModel::Interface + { + public: + virtual Tensor<1,dim> gravity_vector (const Point &pos) const; + + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm); + + double mmm; + }; + + template + Tensor<1,dim> + HollowSphereGravity:: + gravity_vector(const Point &pos) const + { + + const double x=pos[0]; + const double y=pos[1]; + const double z=pos[2]; + + const std_cxx11::array spos = + aspect::Utilities::Coordinates::cartesian_to_spherical_coordinates(pos); + + const double r=spos[0]; + const double phi=spos[1]; + const double theta=spos[2]; + + Tensor<1,dim> g; + + const double gammma = 1.0; + const double R1 = 0.5; + const double R2 = 1.0; + + double alpha,beta,rho; + + if (mmm == -1) + { + alpha = -gammma*(pow(R2,3)-pow(R1,3))/(pow(R2,3)*log(R1)-pow(R1,3)*log(R2)); + beta = -3*gammma*(log(R2)-log(R1))/(pow(R1,3)*log(R2)-pow(R2,3)*log(R1)) ; + rho=-(alpha/pow(r,4)*(8*log(r)-6) + 8./3.*beta/r+8*gammma/pow(r,4))*cos(theta); + } + else + { + alpha=gammma*(mmm+1)*(pow(R1,-3)-pow(R2,-3))/(pow(R1,-mmm-4)-pow(R2,-mmm-4)); + beta=-3*gammma*(pow(R1,mmm+1)-pow(R2,mmm+1))/(pow(R1,mmm+4)-pow(R2,mmm+4)); + rho=-(2*alpha*pow(r,-4)*(mmm+3)/(mmm+1)*(mmm-1)-2*beta/3*(mmm-1)*(mmm+3)*pow(r,mmm)-mmm*(mmm+5)*2*gammma*pow(r,mmm-3) )*cos(theta); + } + + const double g_x= -x/r; + const double g_y= -y/r; + const double g_z= -z/r; + + g[0]=rho*g_x; + g[1]=rho*g_y; + g[2]=rho*g_z; + + return g; + } + + template + void + HollowSphereGravity::declare_parameters (ParameterHandler &prm) + { + //nothing to declare here. This plugin will however, read parameters + //declared by the material model in the "HollowSphere benchmark" section + } + + template + void + HollowSphereGravity::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("HollowSphere benchmark"); + { + mmm = prm.get_double ("Viscosity parameter"); + } + prm.leave_subsection(); + } + + + /** + * A postprocessor that evaluates the accuracy of the solution. + * + * The implementation of error evaluators that correspond to the + * benchmarks defined in the paper Duretz et al. reference above. + */ + template + class HollowSpherePostprocessor : public Postprocess::Interface, public ::aspect::SimulatorAccess + { + public: + /** + * Generate graphical output from the current solution. + */ + virtual + std::pair + execute (TableHandler &statistics); + }; + + template + std::pair + HollowSpherePostprocessor::execute (TableHandler &statistics) + { + std_cxx1x::shared_ptr > ref_func; + { + const HollowSphereMaterial * + material_model + = dynamic_cast *>(&this->get_material_model()); + + ref_func.reset (new AnalyticSolutions::FunctionHollowSphere(material_model->get_mmm())); + } + + const QGauss quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.velocities).degree+2); + + Vector cellwise_errors_u (this->get_triangulation().n_active_cells()); + Vector cellwise_errors_p (this->get_triangulation().n_active_cells()); + Vector cellwise_errors_ul2 (this->get_triangulation().n_active_cells()); + Vector cellwise_errors_pl2 (this->get_triangulation().n_active_cells()); + + ComponentSelectFunction comp_u(std::pair(0,dim), + dim+2); + ComponentSelectFunction comp_p(dim, dim+2); + + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_u, + quadrature_formula, + VectorTools::L1_norm, + &comp_u); + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_p, + quadrature_formula, + VectorTools::L1_norm, + &comp_p); + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_ul2, + quadrature_formula, + VectorTools::L2_norm, + &comp_u); + VectorTools::integrate_difference (this->get_mapping(),this->get_dof_handler(), + this->get_solution(), + *ref_func, + cellwise_errors_pl2, + quadrature_formula, + VectorTools::L2_norm, + &comp_p); + + const double u_l1 = Utilities::MPI::sum(cellwise_errors_u.l1_norm(),this->get_mpi_communicator()); + const double p_l1 = Utilities::MPI::sum(cellwise_errors_p.l1_norm(),this->get_mpi_communicator()); + const double u_l2 = std::sqrt(Utilities::MPI::sum(cellwise_errors_ul2.norm_sqr(),this->get_mpi_communicator())); + const double p_l2 = std::sqrt(Utilities::MPI::sum(cellwise_errors_pl2.norm_sqr(),this->get_mpi_communicator())); + + std::ostringstream os; + + os << std::scientific << u_l1 + << ", " << p_l1 + << ", " << u_l2 + << ", " << p_l2; + + return std::make_pair("Errors u_L1, p_L1, u_L2, p_L2:", os.str()); + } + + } +} + + + +// explicit instantiations +namespace aspect +{ + namespace HollowSphereBenchmark + { + ASPECT_REGISTER_MATERIAL_MODEL(HollowSphereMaterial, + "HollowSphereMaterial", + "A material model that corresponds to the `HollowSphere' benchmark. " + "See the manual for more information.") + + ASPECT_REGISTER_BOUNDARY_VELOCITY_MODEL(HollowSphereBoundary, + "HollowSphereBoundary", + "Implementation of the velocity boundary conditions for the " + "`HollowSphere' benchmark. See the manual for more information about this " + "benchmark.") + + ASPECT_REGISTER_POSTPROCESSOR(HollowSpherePostprocessor, + "HollowSpherePostprocessor", + "A postprocessor that compares the solution of the `HollowSphere' benchmark " + "with the one computed by ASPECT " + "and reports the error. See the manual for more information.") + ASPECT_REGISTER_GRAVITY_MODEL(HollowSphereGravity, + "HollowSphereGravity", + "A gravity model in corresponding to the `HollowSphere' benchmark. " + "See the manual for more information.") + } +} + diff --git a/tests/hollow_sphere.prm b/tests/hollow_sphere.prm new file mode 100644 index 00000000000..a8d0efdb4fd --- /dev/null +++ b/tests/hollow_sphere.prm @@ -0,0 +1,99 @@ +# The 3D HollowSphere benchmark, for which an +# analytical solution is available (Thieulot, in prep). + + +############### Global parameters + +set Additional shared libraries = ./libhollow_sphere.so +set Linear solver tolerance = 1e-12 +set Dimension = 3 +set Start time = 0 +set End time = 0 +set Use years in output instead of seconds = false +set Nonlinear solver scheme = Stokes only +set Output directory = output +set Pressure normalization = surface +set Number of cheap Stokes solver steps = 0 + +############### Parameters describing the model +# The setup is a 3D shell +# Because the temperature plays no role in this model we need not +# bother to describe temperature boundary conditions or +# the material parameters that pertain to the temperature. + +subsection Geometry model + set Model name = spherical shell + + subsection Spherical shell + set Inner radius = 0.5 + set Outer radius = 1.0 + set Cells along circumference = 12 + end + +end + +#Boundary conditions +subsection Model settings + + set Prescribed velocity boundary indicators = inner : HollowSphereBoundary, \ + outer : HollowSphereBoundary +end + + +subsection Material model + set Model name = HollowSphereMaterial +end + + +subsection Gravity model + set Model name = HollowSphereGravity +end + +subsection HollowSphere benchmark + #Viscosity parameter is m + set Viscosity parameter = -1 +end + +############### Parameters describing the temperature field +# As above, there is no need to set anything for the +# temperature boundary conditions. + +subsection Boundary temperature model + set Model name = box +end + +subsection Initial temperature model + set Model name = function + + subsection Function + set Function expression = 0 + end +end + + +############### Parameters describing the discretization +# The following parameters describe how often we want to refine +# the mesh globally and adaptively, what fraction of cells should +# be refined in each adaptive refinement step, and what refinement +# indicator to use when refining the mesh adaptively. + +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 2 + set Refinement fraction = 0.2 + set Strategy = velocity +end + + +############### Parameters describing what to do with the solution +# The final section allows us to choose which postprocessors to +# run at the end of each time step. + +subsection Postprocess + set List of postprocessors = visualization, velocity statistics, HollowSpherePostprocessor, depth average, memory statistics + + subsection Visualization + set List of output variables = density, viscosity, strain rate + end +end + diff --git a/tests/hollow_sphere/screen-output b/tests/hollow_sphere/screen-output new file mode 100644 index 00000000000..c00fb4287e1 --- /dev/null +++ b/tests/hollow_sphere/screen-output @@ -0,0 +1,44 @@ +----------------------------------------------------------------------------- +-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion. +-- . version 2.0.0-pre +-- . running in OPTIMIZED mode +-- . running with 1 MPI process +-- . using Trilinos +----------------------------------------------------------------------------- + +Loading shared library <./libhollow_sphere.so> + +Number of active cells: 768 (on 3 levels) +Number of degrees of freedom: 28,690 (20,790+970+6,930) + +*** Timestep 0: t=0 seconds + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+46 iterations. + Relative Stokes residual after nonlinear iteration 1: 1 + Solving Stokes system... 0+0 iterations. + Relative Stokes residual after nonlinear iteration 2: 2.566e-13 + + Postprocessing: + Writing graphical output: output/solution/solution-00000 + RMS, max velocity: 1.55 m/s, 4.84 m/s + Errors u_L1, p_L1, u_L2, p_L2: 5.111117e-01, 1.866412e-01, 1.767652e-01, 1.552251e-01 + Writing depth average: output/depth_average.gnuplot + System matrix memory consumption: 53.68 MB + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 15.4s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assemble Stokes system | 2 | 2.26s | 15% | +| Build Stokes preconditioner | 1 | 0.531s | 3.5% | +| Solve Stokes system | 2 | 3.58s | 23% | +| Initialization | 1 | 0.031s | 0.2% | +| Postprocessing | 1 | 8.54s | 56% | +| Setup dof systems | 1 | 0.29s | 1.9% | +| Setup initial conditions | 1 | 0.0522s | 0.34% | ++---------------------------------+-----------+------------+------------+ +