Simplify derivative test

Author: gassmoeller

source/material_model/visco_plastic.cc

@@ -183,7 +183,7 @@ namespace aspect
       // The first time this function is called (first iteration of first time step)
       // a specified "reference" strain rate is used as the returned value would
       // otherwise be zero.
-      const double edot_ii = ( (&(this->get_simulator()) != nullptr && this->get_timestep_number() == 0 && strain_rate.norm() <= std::numeric_limits<double>::min())
+      const double edot_ii = ( (this->get_timestep_number() == 0 && strain_rate.norm() <= std::numeric_limits<double>::min())
                                ?
                                ref_strain_rate
                                :
@@ -358,22 +358,11 @@ namespace aspect
       MaterialModel::MaterialModelDerivatives<dim> *derivatives;
       derivatives = out.template get_additional_output<MaterialModel::MaterialModelDerivatives<dim> >();
 
-      // When there is no simulator, we are in a test, so make sure we use the derivatives.
-      double derivative_scaling_factor = 1;
-      if (&(this->get_simulator()) != nullptr)
-        {
-          // If do not use use the Newton solver, get the derivative scaling factor,
-          // otherwise set it to zero.
-          if (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::Newton_Stokes)
-            {
-              derivative_scaling_factor = this->get_newton_handler().get_newton_derivative_scaling_factor();
-            }
-          else
-            {
-              derivative_scaling_factor = 0;
-            }
-        }
-
+      // If we use the Newton solver, get the derivative scaling factor,
+      // otherwise set it to zero.
+      double derivative_scaling_factor = 0.0;
+      if (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::Newton_Stokes)
+        derivative_scaling_factor = this->get_newton_handler().get_newton_derivative_scaling_factor();
 
       // Loop through points
       for (unsigned int i=0; i < in.temperature.size(); ++i)
@@ -458,7 +447,7 @@ namespace aspect
                     }
 
                   /**
-                   * Now compute the derivative of the viscoisty to the pressure
+                   * Now compute the derivative of the viscosity to the pressure
                    */
                   double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy);
 
@@ -634,14 +623,6 @@ namespace aspect
     void
     ViscoPlastic<dim>::declare_parameters (ParameterHandler &prm)
     {
-      prm.enter_subsection("Compositional fields");
-      {
-        prm.declare_entry ("Number of fields", "0",
-                           Patterns::Integer (0),
-                           "The number of fields that will be advected along with the flow field, excluding "
-                           "velocity, pressure and temperature.");
-      }
-      prm.leave_subsection();
       prm.enter_subsection("Material model");
       {
         prm.enter_subsection ("Visco Plastic");
@@ -827,13 +808,7 @@ namespace aspect
     ViscoPlastic<dim>::parse_parameters (ParameterHandler &prm)
     {
       // increment by one for background:
-      unsigned int n_fields = 0;
-      prm.enter_subsection("Compositional fields");
-      {
-        n_fields = prm.get_integer("Number of fields")+1;//this->n_compositional_fields() + 1;
-      }
-      prm.leave_subsection();
-
+      const unsigned int n_fields = this->n_compositional_fields() + 1;
 
       // number of required compositional fields for full finite strain tensor
       const unsigned int s = Tensor<2,dim>::n_independent_components;
@@ -842,7 +817,6 @@ namespace aspect
       {
         prm.enter_subsection ("Visco Plastic");
         {
-
           // Reference and minimum/maximum values
           reference_T = prm.get_double("Reference temperature");
           min_strain_rate = prm.get_double("Minimum strain rate");
@@ -873,6 +847,7 @@ namespace aspect
           if (use_strain_weakening)
             AssertThrow(this->n_compositional_fields() >= 1,
                         ExcMessage("There must be at least one compositional field. "));
+
           use_finite_strain_tensor  = prm.get_bool ("Use finite strain tensor");
           if (use_finite_strain_tensor)
             AssertThrow(this->n_compositional_fields() >= s,

tests/visco_plastic_derivatives.cc

@@ -1,18 +1,26 @@
 #include <aspect/simulator.h>
-#include <deal.II/grid/tria.h>
+#include <aspect/simulator_access.h>
+
 #include <aspect/material_model/interface.h>
 #include <aspect/material_model/visco_plastic.h>
 #include <aspect/simulator_access.h>
 #include <aspect/newton.h>
 
+#include <deal.II/grid/tria.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/std_cxx11/shared_ptr.h>
+
 #include <iostream>
 
-int f(double parameter)
+template <int dim>
+void f(const aspect::SimulatorAccess<dim> &simulator_access,
+       aspect::Assemblers::Manager<dim> &)
 {
+  // Prepare NewtonHandler for test
+  const_cast<aspect::NewtonHandler<dim> &> (simulator_access.get_newton_handler()).set_newton_derivative_scaling_factor(1.0);
 
-  std::cout << std::endl << "Test for p = " << parameter << std::endl;
+  std::cout << std::endl << "Testing ViscoPlastic derivatives against analytical derivatives" << std::endl;
 
-  const int dim=2;
   using namespace aspect::MaterialModel;
   MaterialModelInputs<dim> in_base(5,3);
   in_base.composition[0][0] = 0;
@@ -38,9 +46,9 @@ int f(double parameter)
   in_base.pressure[4] = 5e8;
 
   /**
-   * We can't take to small strain-rates, because then the difference in the
-   * visocisty will be too small for the double accuracy which stores
-   * the visocity solutions and the finite diference solution.
+   * We can't take too small strain-rates, because then the difference in the
+   * viscosity will be too small for the double accuracy which stores
+   * the viscosity solutions and the finite difference solution.
    */
   in_base.strain_rate[0] = SymmetricTensor<2,dim>();
   in_base.strain_rate[0][0][0] = 1e-12;
@@ -125,48 +133,18 @@ int f(double parameter)
   MaterialModelOutputs<dim> out_dviscositydstrainrate_oneone(5,3);
   MaterialModelOutputs<dim> out_dviscositydtemperature(5,3);
 
-  if (out_base.get_additional_output<MaterialModelDerivatives<dim> >() != NULL)
-    throw "error";
-
-  out_base.additional_outputs.push_back(std::make_shared<MaterialModelDerivatives<dim> > (5));
-
-  ViscoPlastic<dim> mat;
-  ParameterHandler prm;
-  mat.declare_parameters(prm);
-
-  prm.enter_subsection("Compositional fields");
-  {
-    prm.set ("Number of fields", "3");
-  }
-  prm.leave_subsection();
-
-  prm.enter_subsection("Material model");
-  {
-    prm.enter_subsection ("Visco Plastic");
-    {
-//      prm.enter_subsection ("Viscosity");
-//      {
-      //prm.set ("Reference strain rate", "1e-20");
-      prm.set ("Angles of internal friction", "30");
-      //    }
-//     prm.leave_subsection();
-    }
-    prm.leave_subsection();
-  }
-  prm.leave_subsection();
-
-  mat.parse_parameters(prm);
+  out_base.additional_outputs.push_back(std_cxx11::make_shared<MaterialModelDerivatives<dim> > (5));
 
-  mat.evaluate(in_base, out_base);
-  mat.evaluate(in_dviscositydpressure, out_dviscositydpressure);
-  mat.evaluate(in_dviscositydstrainrate_zerozero, out_dviscositydstrainrate_zerozero);
-  mat.evaluate(in_dviscositydstrainrate_onezero, out_dviscositydstrainrate_onezero);
-  mat.evaluate(in_dviscositydstrainrate_oneone, out_dviscositydstrainrate_oneone);
-  mat.evaluate(in_dviscositydtemperature, out_dviscositydtemperature);
+  simulator_access.get_material_model().evaluate(in_base, out_base);
+  simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure);
+  simulator_access.get_material_model().evaluate(in_dviscositydstrainrate_zerozero, out_dviscositydstrainrate_zerozero);
+  simulator_access.get_material_model().evaluate(in_dviscositydstrainrate_onezero, out_dviscositydstrainrate_onezero);
+  simulator_access.get_material_model().evaluate(in_dviscositydstrainrate_oneone, out_dviscositydstrainrate_oneone);
+  simulator_access.get_material_model().evaluate(in_dviscositydtemperature, out_dviscositydtemperature);
 
   // set up additional output for the derivatives
   MaterialModelDerivatives<dim> *derivatives;
-  derivatives = out_base.get_additional_output<MaterialModelDerivatives<dim> >();
+  derivatives = out_base.template get_additional_output<MaterialModelDerivatives<dim> >();
 
   double temp;
   for (unsigned int i = 0; i < 5; i++)
@@ -223,7 +201,7 @@ int f(double parameter)
       std::cout << "onezero at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][0]   << std::endl;
       if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][0]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][0])) )
         {
-          std::cout << "   Error: The derivative of the viscosity to the strain rate is too different from the analitical value." << std::endl;
+          std::cout << "   Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl;
           Error = true;
         }
     }
@@ -241,34 +219,40 @@ int f(double parameter)
       std::cout << "oneone at point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][1]  << std::endl;
       if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][1]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][1])) )
         {
-          std::cout << "   Error: The derivative of the viscosity to the strain rate is too different from the analitical value." << std::endl;
+          std::cout << "   Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl;
           Error = true;
         }
 
     }
 
   if (Error)
     {
-      std::cout << "Some parts of the test where not succesful." << std::endl;
+      std::cout << "Some parts of the test where not successful." << std::endl;
     }
   else
     {
       std::cout << "OK" << std::endl;
     }
 
-  return 42;
+  // Restore NewtonHandler state after test
+  const_cast<aspect::NewtonHandler<dim> &> (simulator_access.get_newton_handler()).set_newton_derivative_scaling_factor(0.0);
+
+}
+
+template <>
+void f(const aspect::SimulatorAccess<3> &simulator_access,
+       aspect::Assemblers::Manager<3> &)
+{
+  AssertThrow(false,dealii::ExcInternalError());
 }
 
-int exit_function()
+template <int dim>
+void signal_connector (aspect::SimulatorSignals<dim> &signals)
 {
-  exit(0);
-  return 42;
+  std::cout << "* Connecting signals" << std::endl;
+  signals.set_assemblers.connect (&f<dim>);
 }
-// run this function by initializing a global variable by it
-int ik = f(-1); // Testing harmonic mean
-int ji = f(0); // Testing geometric mean
-int jj = f(1); // Testing arithmetic mean
-int kj = f(1000); // Testing max function
-int kl = exit_function();
 
+ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>,
+                                  signal_connector<3>)
 

tests/visco_plastic_derivatives.prm

@@ -2,3 +2,106 @@
 # material model are the same as the finite difference derivatives 
 # created by the drucker prager material model.
 set Additional shared libraries = tests/libvisco_plastic_derivatives.so
+
+set Dimension                              = 2
+set End time                               = 0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme                = Newton Stokes
+set Max nonlinear iterations               = 1
+set Number of cheap Stokes solver steps    = 0
+
+# Model geometry (100x100 km, 10 km spacing)
+subsection Geometry model
+  set Model name = box
+  subsection Box
+    set X repetitions = 10
+    set Y repetitions = 10
+    set X extent      = 100e3
+    set Y extent      = 100e3
+  end
+end
+
+# Mesh refinement specifications 
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 0
+  set Time steps between mesh refinement = 0
+end
+
+
+# Boundary classifications (fixed T boundaries, prescribed velocity) 
+subsection Model settings
+  set Fixed temperature boundary indicators   = bottom, top, left, right
+  set Prescribed velocity boundary indicators = bottom y: function, top y: function, left x: function, right x: function
+end
+
+# Velocity on boundaries characterized by functions
+subsection Boundary velocity model
+  subsection Function
+    set Variable names      = x,y
+    set Function constants  = m=0.0005, year=1
+    set Function expression = if (x<50e3 , -1*m/year, 1*m/year); if (y<50e3 , 1*m/year, -1*m/year);
+  end
+end
+
+# Temperature boundary and initial conditions
+subsection Boundary temperature model
+  set List of model names = box
+  subsection Box
+    set Bottom temperature = 273
+    set Left temperature   = 273
+    set Right temperature  = 273
+    set Top temperature    = 273
+  end
+end
+subsection Initial temperature model
+  set Model name = function
+  subsection Function
+    set Function expression = 273
+  end
+end
+
+# Compositional fields used to track finite strain invariant
+subsection Compositional fields
+  set Number of fields = 3
+end
+
+# We prescribe some initial strain at the center of the domain
+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Variable names      = x,y
+    set Function expression = if(x>=45e3&x<=55e3&y>=45.3e3&y<=55.e3,0.2,0);0;0
+  end
+end
+
+# Boundary composition specification
+subsection Boundary composition model
+  set Model name = initial composition
+end
+
+# Material model (values for background material)
+subsection Material model
+  set Model name = visco plastic
+  subsection Visco Plastic
+    set Angles of internal friction = 30.
+  end
+end
+
+# Gravity model
+subsection Gravity model
+  set Model name = vertical
+  subsection Vertical
+    set Magnitude = 10.0
+  end
+end
+
+# Post processing
+# named additional outputs includes the weakened cohesions and friction angles
+subsection Postprocess
+  set List of postprocessors = velocity statistics, mass flux statistics, visualization
+  subsection Visualization
+    set List of output variables = viscosity, strain rate, named additional outputs
+    set Output format            = gnuplot
+  end
+end