Merge pull request #5998 from bobmyhill/composite_phases

added phases test for CompositeViscoPlastic

tests/composite_viscous_outputs_phases.cc

@@ -0,0 +1,261 @@
+/*
+  Copyright (C) 2024 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include <aspect/simulator.h>
+#include <aspect/material_model/utilities.h>
+#include <aspect/material_model/rheology/composite_visco_plastic.h>
+#include <aspect/simulator_signals.h>
+
+template <int dim>
+void f(const aspect::SimulatorAccess<dim> &simulator_access,
+       aspect::Assemblers::Manager<dim> &)
+{
+  // This function tests whether the composite creep rheology is producing
+  // the correct composite viscosity and partial strain rates corresponding to
+  // the different creep mechanisms incorporated into the rheology.
+  // It is assumed that each individual creep mechanism has already been tested.
+
+  using namespace aspect::MaterialModel;
+
+  // First, we set up a few objects which are used by the rheology model.
+  aspect::ParameterHandler prm;
+  const std::vector<std::string> list_of_composition_names = simulator_access.introspection().get_composition_names();
+  MaterialUtilities::PhaseFunction<dim> phase_function;
+  phase_function.initialize_simulator (simulator_access.get_simulator());
+  phase_function.declare_parameters (prm);
+  prm.set("Define transition by depth instead of pressure", "false");
+  prm.set("Phase transition pressures", "3e9");
+  prm.set("Phase transition pressure widths", "1e9");
+  prm.set("Phase transition temperatures", "273");
+  prm.set("Phase transition Clapeyron slopes", "0");
+  phase_function.parse_parameters (prm);
+
+  std::vector<unsigned int> n_phases_for_each_composition = phase_function.n_phases_for_each_composition();
+  // Currently, phase_function.n_phases_for_each_composition() returns a list of length
+  // equal to the total number of compositions, whether or not they are chemical compositions.
+  // The equation_of_state (multicomponent incompressible) requires a list only for
+  // chemical compositions.
+  std::vector<unsigned int> n_phases_for_each_chemical_composition = {n_phases_for_each_composition[0]};
+  std::vector<unsigned int> n_phase_transitions_for_each_chemical_composition = {n_phases_for_each_composition[0] - 1};
+  unsigned int n_phases = n_phases_for_each_composition[0];
+  for (auto i : simulator_access.introspection().chemical_composition_field_indices())
+    {
+      n_phases_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1]);
+      n_phase_transitions_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1] - 1);
+      n_phases += n_phases_for_each_composition[i+1];
+    }
+
+  const unsigned int composition = 0;
+  const std::vector<double> volume_fractions = {1.};
+  std::vector<double> phase_function_values = {0.};
+  const std::vector<unsigned int> n_phase_transitions_per_composition = n_phase_transitions_for_each_chemical_composition;
+
+  // Next, we initialise instances of the composite rheology and
+  // individual creep mechanisms.
+  std::unique_ptr<Rheology::CompositeViscoPlastic<dim>> composite_creep;
+  composite_creep = std::make_unique<Rheology::CompositeViscoPlastic<dim>>();
+  composite_creep->initialize_simulator (simulator_access.get_simulator());
+  composite_creep->declare_parameters(prm);
+  MaterialUtilities::PhaseFunction<dim>::declare_parameters(prm);
+  prm.set("Viscosity averaging scheme", "isostrain");
+  prm.set("Include diffusion creep in composite rheology", "true");
+  prm.set("Include dislocation creep in composite rheology", "true");
+  prm.set("Include Peierls creep in composite rheology", "true");
+  prm.set("Include Drucker Prager plasticity in composite rheology", "true");
+  prm.set("Peierls creep flow law", "viscosity approximation");
+  prm.set("Cohesions", "background:1e9|5e8");
+  prm.set("Maximum yield stress", "5e10");
+  composite_creep->parse_parameters(prm, std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+
+  std::unique_ptr<Rheology::DiffusionCreep<dim>> diffusion_creep;
+  diffusion_creep = std::make_unique<Rheology::DiffusionCreep<dim>>();
+  diffusion_creep->initialize_simulator (simulator_access.get_simulator());
+  diffusion_creep->declare_parameters(prm);
+  diffusion_creep->parse_parameters(prm, std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+
+  std::unique_ptr<Rheology::DislocationCreep<dim>> dislocation_creep;
+  dislocation_creep = std::make_unique<Rheology::DislocationCreep<dim>>();
+  dislocation_creep->initialize_simulator (simulator_access.get_simulator());
+  dislocation_creep->declare_parameters(prm);
+  dislocation_creep->parse_parameters(prm, std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+
+  std::unique_ptr<Rheology::PeierlsCreep<dim>> peierls_creep;
+  peierls_creep = std::make_unique<Rheology::PeierlsCreep<dim>>();
+  peierls_creep->initialize_simulator (simulator_access.get_simulator());
+  peierls_creep->declare_parameters(prm);
+  peierls_creep->parse_parameters(prm, std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+
+  std::unique_ptr<Rheology::DruckerPragerPower<dim>> drucker_prager_power;
+  drucker_prager_power = std::make_unique<Rheology::DruckerPragerPower<dim>>();
+  drucker_prager_power->initialize_simulator (simulator_access.get_simulator());
+  drucker_prager_power->declare_parameters(prm);
+  prm.set("Cohesions", "background:1e9|5e8");
+  prm.set("Maximum yield stress", "5e10");
+  drucker_prager_power->parse_parameters(prm, std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+
+  // The creep components are arranged in series with each other.
+  // This package of components is then arranged in parallel with
+  // a strain rate limiter with a constant viscosity lim_visc.
+  // The whole system is then arranged in series with a viscosity limiter with
+  // viscosity max_visc.
+  // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc)
+  double min_visc = prm.get_double("Minimum viscosity");
+  double max_visc = prm.get_double("Maximum viscosity");
+  double lim_visc = (min_visc*max_visc)/(max_visc - min_visc);
+
+  // Assign values to the variables which will be passed to compute_viscosity
+  // The test involves pure shear calculations at variable pressure and temperature
+  double temperature;
+  double pressure;
+  const double grain_size = 1.e-3;
+  SymmetricTensor<2,dim> strain_rate;
+  strain_rate[0][0] = -1e-11;
+  strain_rate[0][1] = 0.;
+  strain_rate[1][1] = 1e-11;
+  strain_rate[2][0] = 0.;
+  strain_rate[2][1] = 0.;
+  strain_rate[2][2] = 0.;
+
+  std::cout << "temperature (K)   phase transition progress   eta (Pas)   creep stress (Pa)   edot_ii (/s)   edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl;
+
+  // Loop through strain rates, tracking whether there is a discrepancy in
+  // the decomposed strain rates.
+  bool error = false;
+  double viscosity;
+  double total_strain_rate;
+  double creep_strain_rate;
+  double creep_stress;
+  double diff_stress;
+  double disl_stress;
+  double prls_stress;
+  double drpr_stress;
+  std::vector<double> partial_strain_rates(5, 0.);
+
+  for (unsigned int i=0; i <= 2; i++)
+    {
+      pressure = 1.e9 + i*2.e9;
+      std::cout << "pressure: " << pressure / 1.e9 << " GPa" << std::endl;
+      for (unsigned int j = 0; j <= 10; j++)
+        {
+          temperature = 1000. + j*100.;
+
+          // Compute the phase function values
+          // The depth and gravity are set to zero because they are unused
+          // when phase functions are calculated by pressure.
+          // The phase index is set to invalid_unsigned_int, because it is only used internally
+          // in phase_average_equation_of_state_outputs to loop over all existing phases
+          MaterialUtilities::PhaseFunctionInputs<dim> phase_inputs(temperature,
+                                                                   pressure,
+                                                                   0., 0.,
+                                                                   numbers::invalid_unsigned_int);
+
+          // Compute value of phase functions
+          for (unsigned int j=0; j < phase_function.n_phase_transitions(); ++j)
+            {
+              phase_inputs.phase_index = j;
+              phase_function_values[j] = phase_function.compute_value(phase_inputs);
+            }
+
+          Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition);
+
+          // Compute the viscosity
+          viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates, phase_function_values, n_phase_transitions_per_composition);
+          total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
+
+          // The creep strain rate is calculated by subtracting the strain rate
+          // of the max viscosity dashpot from the total strain rate
+          // The creep stress is then calculated by subtracting the stress running
+          // through the strain rate limiter from the total stress
+          creep_strain_rate = total_strain_rate - partial_strain_rates[4];
+          creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate);
+
+          // Print the output
+          std::cout << temperature << ' ' << phase_function_values[0] << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate;
+          for (unsigned int i=0; i < partial_strain_rates.size(); ++i)
+            {
+              std::cout << ' ' << partial_strain_rates[i]/total_strain_rate;
+            }
+          std::cout << std::endl;
+
+          // The following lines test that each individual creep mechanism
+          // experiences the same creep stress
+
+          // Each creep mechanism should experience the same stress
+          diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition);
+          disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition);
+          prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition);
+          if (partial_strain_rates[3] > 0.)
+            {
+              drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion,
+                            p.angle_internal_friction,
+                            pressure,
+                            partial_strain_rates[3],
+                            p.max_yield_stress);
+            }
+          else
+            {
+              drpr_stress = creep_stress;
+            }
+
+          if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6)
+              || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6)
+              || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6)
+              || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6))
+            {
+              error = true;
+              std::cout << "   creep stress: " << creep_stress;
+              std::cout << " diffusion stress: " << diff_stress;
+              std::cout << " dislocation stress: " << disl_stress;
+              std::cout << " peierls stress: " << prls_stress;
+              std::cout << " drucker prager stress: " << drpr_stress << std::endl;
+            }
+        }
+    }
+
+  if (error)
+    {
+      std::cout << "   Error: The individual creep stresses differ by more than the required tolerance." << std::endl;
+      std::cout << "Some parts of the test were not successful." << std::endl;
+    }
+  else
+    {
+      std::cout << "OK" << std::endl;
+    }
+}
+
+template <>
+void f(const aspect::SimulatorAccess<2> &,
+       aspect::Assemblers::Manager<2> &)
+{
+  AssertThrow(false,dealii::ExcInternalError());
+}
+
+template <int dim>
+void signal_connector (aspect::SimulatorSignals<dim> &signals)
+{
+  using namespace dealii;
+  std::cout << "* Connecting signals" << std::endl;
+  signals.set_assemblers.connect (std::bind(&f<dim>,
+                                            std::placeholders::_1,
+                                            std::placeholders::_2));
+}
+
+ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>,
+                                  signal_connector<3>)

tests/composite_viscous_outputs_phases.prm

@@ -0,0 +1,80 @@
+# This test checks whether the composite viscous rheology
+# plugin produces the correct composite viscosity and
+# decomposed strain rates when a phase transition is included.
+
+set Additional shared libraries = tests/libcomposite_viscous_outputs_phases.so
+set Dimension                              = 3
+set End time                               = 0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme                = no Advection, no Stokes
+
+# 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)
+
+# Temperature boundary and initial conditions
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators   = bottom, top, left, right
+  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
+
+# Velocity on boundaries characterized by functions
+subsection Boundary velocity model
+  set Prescribed velocity boundary indicators = bottom y: function, top y: function, left x: function, right x: function
+
+  subsection Function
+    set Variable names      = x,y,z
+    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);0
+  end
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 273
+  end
+end
+
+# Material model
+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

tests/composite_viscous_outputs_phases/screen-output

@@ -0,0 +1,53 @@
+
+Loading shared library <./libcomposite_viscous_outputs_phases.debug.so>
+
+* Connecting signals
+temperature (K)   phase transition progress   eta (Pas)   creep stress (Pa)   edot_ii (/s)   edot_ii fractions (diff, disl, prls, drpr, max)
+pressure: 1 GPa
+1000 0.0179862 4.49858e+19 8.97717e+08 1e-11 1.6913e-06 0.00916831 0.99083 5.3662e-13 4.49858e-09
+1100 0.0179862 2.98905e+19 5.9581e+08 1e-11 7.23388e-05 0.836312 0.163616 6.73098e-22 2.98905e-09
+1200 0.0179862 7.70568e+18 1.52114e+08 1e-11 0.000594407 0.999399 6.44893e-06 1.51744e-51 7.70568e-10
+1300 0.0179862 2.39282e+18 4.58564e+07 1e-11 0.00338081 0.996619 3.19443e-09 1.38894e-77 2.39282e-10
+1400 0.0179862 9.18168e+17 1.63634e+07 1e-11 0.0149609 0.985039 1.26589e-11 5.83892e-100 9.18169e-11
+1500 0.0179862 4.32083e+17 6.64167e+06 1e-11 0.0538305 0.94617 2.09885e-13 1.5366e-119 4.32084e-11
+1600 0.0179862 2.47246e+17 2.94493e+06 1e-11 0.161077 0.838923 8.74429e-15 3.3607e-137 2.47246e-11
+1700 0.0179862 1.67379e+17 1.34758e+06 1e-11 0.397345 0.602655 5.42942e-16 3.551e-154 1.67378e-11
+1800 0.0179862 1.28447e+17 568948 1e-11 0.749975 0.250025 2.2682e-17 6.70367e-173 1.28447e-11
+1900 0.0179862 1.09563e+17 191256 1e-11 0.962698 0.0373021 2.59166e-19 1.42377e-196 1.09563e-11
+2000 0.0179862 1.02964e+17 59280.1 1e-11 0.99654 0.00345962 2.26225e-21 5.22568e-222 1.02965e-11
+pressure: 3 GPa
+1000 0.5 5.3691e+19 1.07182e+09 1e-11 4.76869e-07 0.000587744 0.99515 0.00426168 5.3691e-09
+1100 0.5 4.29097e+19 8.56193e+08 1e-11 2.79907e-05 0.139287 0.860685 5.65005e-08 4.29097e-09
+1200 0.5 1.70541e+19 3.39081e+08 1e-11 0.000398001 0.998688 0.000914421 4.35181e-28 1.70541e-09
+1300 0.5 4.9077e+18 9.6154e+07 1e-11 0.00233578 0.997664 1.07007e-07 1.86999e-55 4.9077e-10
+1400 0.5 1.72879e+18 3.25758e+07 1e-11 0.0106234 0.989377 1.4219e-10 5.86483e-79 1.72879e-10
+1500 0.5 7.33471e+17 1.26694e+07 1e-11 0.0392318 0.960768 1.03038e-12 1.82512e-99 7.3347e-11
+1600 0.5 3.72285e+17 5.4457e+06 1e-11 0.120854 0.879146 2.33681e-14 8.43651e-118 3.72285e-11
+1700 0.5 2.2329e+17 2.46581e+06 1e-11 0.311075 0.688925 9.90838e-16 5.26533e-135 2.2329e-11
+1800 0.5 1.53962e+17 1.07924e+06 1e-11 0.638071 0.361929 3.94829e-17 6.01736e-153 1.53962e-11
+1900 0.5 1.19579e+17 391585 1e-11 0.922158 0.0778423 5.25907e-19 5.8183e-175 1.19579e-11
+2000 0.5 1.06072e+17 121442 1e-11 0.992096 0.0079044 3.18167e-21 2.19751e-200 1.06071e-11
+pressure: 5 GPa
+1000 0.982014 4.06599e+19 8.11198e+08 1e-11 8.52316e-08 7.64121e-06 0.00625559 0.993737 4.06599e-09
+1100 0.982014 4.0594e+19 8.0988e+08 1e-11 7.12921e-06 0.00536781 0.0784945 0.916131 4.0594e-09
+1200 0.982014 3.57677e+19 7.13354e+08 1e-11 0.000251506 0.814985 0.183155 0.0016077 3.57677e-09
+1300 0.982014 1.01801e+19 2.01602e+08 1e-11 0.00161362 0.998378 7.89782e-06 5.82951e-31 1.01801e-09
+1400 0.982014 3.34138e+18 6.48275e+07 1e-11 0.00754067 0.992459 3.15133e-09 1.3456e-55 3.34138e-10
+1500 0.982014 1.30692e+18 2.41383e+07 1e-11 0.0285575 0.971443 9.10293e-12 4.74544e-77 1.30691e-10
+1600 0.982014 6.01678e+17 1.00336e+07 1e-11 0.0903463 0.909654 1.03399e-13 4.10849e-96 6.01679e-11
+1700 0.982014 3.23232e+17 4.46465e+06 1e-11 0.240983 0.759017 2.73137e-15 1.07226e-113 3.23232e-11
+1800 0.982014 1.99404e+17 1.98807e+06 1e-11 0.527179 0.472821 9.07293e-17 2.90092e-131 1.99404e-11
+1900 0.982014 1.38804e+17 776080 1e-11 0.855044 0.144956 1.43049e-18 1.08706e-151 1.38805e-11
+2000 0.982014 1.12372e+17 247431 1e-11 0.982283 0.0177172 7.09191e-21 1.63545e-176 1.12372e-11
+OK
+Number of active cells: 100 (on 1 levels)
+Number of degrees of freedom: 5,534 (3,969+242+1,323)
+
+*** Timestep 0:  t=0 years, dt=0 years
+
+   Postprocessing:
+
+Termination requested by criterion: end time
+
+
+