Rheology and conductivity modification for the entropy plugin

benchmarks/entropy_adiabat/entropy_adiabat.prm

@@ -163,8 +163,10 @@ subsection Material model
   set Model name = entropy model
 
   subsection Entropy model
-    set Data directory                   = $ASPECT_SOURCE_DIR/data/material-model/
-    set Material file name               = entropy-table/opxtable/material_table.txt
+    set Data directory                   = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/opxtable/
+    set Material file name               = material_table.txt
+    set Lateral viscosity file name      = constant_lateral_vis_prefactor.txt
+    set Thermal conductivity formulation = constant
     set Thermal conductivity             = 0.1
   end
 

benchmarks/entropy_adiabat/entropy_combined.prm

@@ -43,6 +43,7 @@ subsection Material model
   set Model name = entropy model
 
   subsection Entropy model
+    set Thermal conductivity formulation = constant
     set Thermal conductivity             = 50
   end
 end

benchmarks/entropy_adiabat/entropy_conduction.prm

@@ -128,10 +128,12 @@ subsection Material model
   set Model name = entropy model
 
   subsection Entropy model
-    set Data directory                   = $ASPECT_SOURCE_DIR/data/material-model/
-    set Material file name               = entropy-table/opxtable/material_table.txt
+    set Data directory                   = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/opxtable/
+    set Material file name               = material_table.txt
     set Maximum viscosity                = 1e30
     set Minimum viscosity                = 1e30
+    set Lateral viscosity file name      = constant_lateral_vis_prefactor.txt
+    set Thermal conductivity formulation = constant
     set Thermal conductivity             = 4.7
   end
 end

benchmarks/entropy_adiabat/entropy_half_space.prm

@@ -136,10 +136,12 @@ subsection Material model
   set Model name = entropy model
 
   subsection Entropy model
-    set Data directory                   = $ASPECT_SOURCE_DIR/data/material-model/
-    set Material file name               = entropy-table/opxtable/material_table.txt
+    set Data directory                   = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/opxtable/
+    set Material file name               = material_table.txt
     set Maximum viscosity                = 1e30
     set Minimum viscosity                = 1e30
+    set Lateral viscosity file name      = constant_lateral_vis_prefactor.txt
+    set Thermal conductivity formulation = constant
     set Thermal conductivity             = 4.7
   end
 end

benchmarks/entropy_adiabat/plugins/entropy_model.h

@@ -26,13 +26,14 @@
 #include <aspect/utilities.h>
 #include <aspect/simulator_access.h>
 #include <aspect/material_model/rheology/ascii_depth_profile.h>
+#include <aspect/material_model/rheology/drucker_prager.h>
+#include <aspect/material_model/steinberger.h>
 
 namespace aspect
 {
   namespace MaterialModel
   {
     using namespace dealii;
-
     /**
      * A material model that is designed to use pressure and entropy (rather
      * than pressure and temperature) as independent variables. It will look up
@@ -44,6 +45,7 @@ namespace aspect
      * the Stokes equation).
      * @ingroup MaterialModels
      */
+
     template <int dim>
     class EntropyModel: public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
     {
@@ -110,30 +112,67 @@ namespace aspect
 
 
       private:
+
+        double angle_of_internal_friction;
+        double cohesion;
+
         /**
          * Minimum/Maximum viscosity and lateral viscosity variations.
          */
-        double lateral_viscosity_prefactor;
         double min_eta;
         double max_eta;
+        double max_lateral_eta_variation;
 
         /**
-         * The value for thermal conductivity. This model only
-         * implements a constant thermal conductivity for the whole domain.
+         * The value for thermal conductivity. It can be a constant
+         * for the whole domain, or P-T dependent.
          */
         double thermal_conductivity_value;
+        double thermal_conductivity (const double temperature,
+                                     const double pressure,
+                                     const Point<dim> &position) const;
+
+        enum ConductivityFormulation
+        {
+          constant,
+          p_T_dependent
+        } conductivity_formulation;
+
+        /**
+         * Parameters for the temperature- and pressure dependence of the
+         * thermal conductivity.
+         */
+        std::vector<double> conductivity_transition_depths;
+        std::vector<double> reference_thermal_conductivities;
+        std::vector<double> conductivity_pressure_dependencies;
+        std::vector<double> conductivity_reference_temperatures;
+        std::vector<double> conductivity_exponents;
+        std::vector<double> saturation_scaling;
+        double maximum_conductivity;
 
         /**
          * Information about the location of data files.
          */
         std::string data_directory;
         std::string material_file_name;
+        std::string lateral_viscosity_file_name;
 
         /**
          * Pointer to the StructuredDataLookup object that holds the material data.
          */
         std::unique_ptr<Utilities::StructuredDataLookup<2>> material_lookup;
 
+        /**
+         * Pointer to an object that reads and processes data for the lateral
+         * temperature dependency of viscosity.
+         */
+        std::unique_ptr<internal::LateralViscosityLookup> lateral_viscosity_prefactor_lookup;
+
+        /**
+         * Objects for computing plastic stresses, viscosities, and additional outputs
+         */
+        Rheology::DruckerPrager<dim> drucker_prager_plasticity;
+
         /**
          * Pointer to the rheology model used for depth-dependence from an
          * ascii file

data/material-model/entropy-table/opxtable/constant_lateral_vis_prefactor.txt

@@ -0,0 +1,3 @@
+# H/nR [K], depth [km]
+0 0.0000000e+00
+0 3.0000000e+06

data/material-model/entropy-table/opxtable/constant_lateral_vis_prefactor_50000.txt

@@ -0,0 +1,3 @@
+# H/nR [K], depth [km]
+50000 0.0000000e+00
+50000 3.0000000e+06

doc/modules/changes/20230712_ranpengli

@@ -0,0 +1,4 @@
+Changed: The entropy plugin now includes a plasticity rheology, which can be set by cohesion and frictional angle. It can also use a viscosity prefactor profile, set the lateral viscosity variation limit, and choose p-T dependent conductivity
+
+<br>
+(Ranpeng Li & Juliane Dannberg, 2023/07/13)

tests/entropy_plasticity.cc

@@ -0,0 +1,24 @@
+/*
+  Copyright (C) 2022 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 "../benchmarks/entropy_adiabat/plugins/compute_entropy_profile.cc"
+#include "../benchmarks/entropy_adiabat/plugins/entropy_advection.cc"
+#include "../benchmarks/entropy_adiabat/plugins/entropy_model.cc"
+#include "../benchmarks/entropy_adiabat/plugins/prescribed_temperature.cc"

tests/entropy_plasticity.prm

@@ -0,0 +1,159 @@
+# A box model for convection under early-Earth conditions, with plasticity,
+# reference viscosity profile, and p-T dependent conductivity
+
+set Dimension                              = 2
+set Use years in output instead of seconds = true
+set End time                               = 4e5
+set Output directory                       = output-entropy-plasticity
+set Nonlinear solver scheme                = iterated Advection and Stokes
+
+subsection Formulation
+  set Mass conservation = projected density field
+end
+
+subsection Discretization
+  subsection Stabilization parameters
+    set List of compositional fields with disabled boundary entropy viscosity = entropy
+  end
+end
+
+subsection Material model
+  set Material averaging = harmonic average only viscosity
+end
+
+subsection Geometry model
+  set Model name = spherical shell
+
+  subsection Spherical shell
+    set Inner radius = 3481000
+    set Outer radius = 6371000
+  end
+end
+
+subsection Heating model
+  set List of model names = shear heating
+end
+
+# Set the reference profiles, using a constant entropy
+subsection Adiabatic conditions model
+  set Model name = compute entropy profile
+
+  subsection Compute entropy profile
+    # This is specific entropy, equivalent to T=1600K according to the used table
+    set Surface entropy = 2535
+  end
+end
+
+subsection Initial temperature model
+  set Model name = prescribed temperature
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = top, bottom
+  set List of model names =  initial temperature
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = top, bottom
+end
+
+subsection Nullspace removal
+  set Remove nullspace = net rotation
+end
+
+subsection Gravity model
+  set Model name = radial constant
+
+  subsection Radial constant
+    set Magnitude = 9.81
+  end
+end
+
+# The temperature field is computed from the entropy
+subsection Temperature field
+  set Temperature method = prescribed field with diffusion
+end
+
+# Density is computed from entropy
+subsection Compositional fields
+  set Number of fields = 2
+  set Names of fields = entropy, density_field
+  set Compositional field methods = field, prescribed field
+end
+
+# Set the initial temperature profile
+subsection Initial composition model
+  set List of model names = function
+
+  subsection Function
+    set Coordinate system = spherical
+    set Variable names = r,phi
+    set Function expression = 2535.08 + 10 * sin(2*phi)*sin(pi*(r-6371000)/2890000) - 1878.618 * erfc((6371000-r)/(2*sqrt(1.152e-6*5e7*31557600))) + \
+    486.368 * max(erfc((r-3481000)/(2*sqrt(1.152e-6*5e8*31557600))), exp(-((phi - 3/4*pi)^2)/(2*(1/35*pi)^2) - ((r-3481000)^2)/(2*400000^2))  );0
+
+    #set Function expression = 2687.748 + 10 * sin(2*phi)*sin(pi*(r-6371000)/2890000) - 2031.286 * erfc((6371000-r)/(2*sqrt(1.152e-6*5e7*31557600))) + 333.7 * erfc((r-3481000)/(2*sqrt(1.152e-6*5e8*31557600)));0
+  end
+end
+
+# Set the boundary entropy to values corresponding to
+# realistic temperature.
+subsection Boundary composition model
+  set Fixed composition boundary indicators = top, bottom
+  set List of model names = spherical constant
+
+  subsection Spherical constant
+    # Corresponding to T=300 K for p=0 GPa according to used table
+    set Outer composition = 600
+    # Corresponding to T=3400 K for p=140 GPa according to used table
+    set Inner composition = 2950
+  end
+end
+
+# This material model reads in the Hefesto table
+subsection Material model
+  set Model name = entropy model
+
+  subsection Entropy model
+    set Data directory                   = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/opxtable/
+    set Material file name               = material_table.txt
+    set Lateral viscosity file name      = constant_lateral_vis_prefactor_50000.txt
+
+    set Thermal conductivity formulation = p-T-dependent
+
+    set Maximum viscosity                = 1e30
+    set Minimum viscosity                = 1e12
+    set Maximum lateral viscosity variation          = 2000
+
+    set Cohesion                         = 8e7
+    set Angle of internal friction       = 0.01
+  end
+
+  subsection Depth dependent viscosity
+    set Data directory                   = $ASPECT_SOURCE_DIR/benchmarks/entropy_adiabat/data/
+    set Data file name                   = constant_viscosity.txt
+  end
+end
+
+subsection Mesh refinement
+  set Initial global refinement                = 3
+  set Initial adaptive refinement              = 0
+  set Minimum refinement level                 = 2
+  set Time steps between mesh refinement       = 0
+  set Strategy                                 = composition threshold
+
+  subsection Composition threshold
+    set Compositional field thresholds = 2700, 1e6
+  end
+end
+
+subsection Checkpointing
+  set Time between checkpoint = 600
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization, velocity statistics, temperature statistics, mass flux statistics, composition statistics
+end
+
+subsection Boundary composition model
+   set Allow fixed composition on outflow boundaries = true
+end