Add compressible EOS option to visco plastic

Author: naliboff

include/aspect/material_model/visco_plastic.h

@@ -24,6 +24,7 @@
 #include <aspect/simulator_access.h>
 #include <aspect/material_model/interface.h>
 #include <aspect/material_model/equation_of_state/multicomponent_incompressible.h>
+#include <aspect/material_model/equation_of_state/multicomponent_compressible.h>
 #include <aspect/material_model/rheology/visco_plastic.h>
 
 #include<deal.II/fe/component_mask.h>
@@ -258,7 +259,17 @@ namespace aspect
         /**
          * Object for computing the equation of state.
          */
-        EquationOfState::MulticomponentIncompressible<dim> equation_of_state;
+        EquationOfState::MulticomponentIncompressible<dim> equation_of_state_incompressible;
+
+        /**
+         * Object for computing the equation of state.
+         */
+        EquationOfState::MulticomponentCompressible<dim> equation_of_state_compressible;
+
+        /**
+         * Whether to use a compressible or incompressible equation of state.
+         */
+        bool use_compressible_equation_of_state;
 
         /**
          * Object that handles phase transitions.

source/material_model/visco_plastic.cc

@@ -68,7 +68,7 @@ namespace aspect
           else
             {
               EquationOfStateOutputs<dim> eos_outputs_all_phases (n_phases);
-              equation_of_state.evaluate(in, 0, eos_outputs_all_phases);
+              equation_of_state_incompressible.evaluate(in, 0, eos_outputs_all_phases);
               reference_density = eos_outputs_all_phases.densities[0];
             }
 
@@ -121,7 +121,7 @@ namespace aspect
       for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
         {
           // First compute the equation of state variables and thermodynamic properties
-          equation_of_state.evaluate(in, i, eos_outputs_all_phases);
+          equation_of_state_incompressible.evaluate(in, i, eos_outputs_all_phases);
 
           const double gravity_norm = this->get_gravity_model().gravity_vector(in.position[i]).norm();
           const double reference_density = (this->get_adiabatic_conditions().is_initialized())
@@ -308,7 +308,7 @@ namespace aspect
     ViscoPlastic<dim>::
     is_compressible () const
     {
-      return equation_of_state.is_compressible();
+      return equation_of_state_incompressible.is_compressible();
     }
 
 
@@ -342,6 +342,15 @@ namespace aspect
 
           EquationOfState::MulticomponentIncompressible<dim>::declare_parameters (prm);
 
+          EquationOfState::MulticomponentCompressible<dim>::declare_parameters (prm);
+
+          prm.declare_entry ("Use compressible equation of state","false",
+                             Patterns::Bool (),
+                             "Whether to use an incompressible or compressible equation of state. "
+                             "If set to true, the material model will switch from using the "
+                             "multicomponent incompressible multicomponent compressible equation "
+                             "of state model. ");
+
           Rheology::ViscoPlastic<dim>::declare_parameters(prm);
 
           // Equation of state parameters
@@ -389,9 +398,26 @@ namespace aspect
           n_phases = phase_function.n_phases_over_all_chemical_compositions();
 
           // Equation of state parameters
-          equation_of_state.initialize_simulator (this->get_simulator());
-          equation_of_state.parse_parameters (prm,
-                                              std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+          use_compressible_equation_of_state = prm.get_bool ("Use compressible equation of state");
+          if (use_compressible_equation_of_state == false)
+            {
+              equation_of_state_incompressible.initialize_simulator (this->get_simulator());
+              equation_of_state_incompressible.parse_parameters (prm,
+                                                                 std::make_unique<std::vector<unsigned int>>(n_phases_for_each_chemical_composition));
+            }
+          else
+            {
+              prm.enter_subsection ("Equation of State");
+              {
+                prm.enter_subsection ("Multicomponent Compressible");
+                {
+                  equation_of_state_compressible.initialize_simulator (this->get_simulator());
+                  equation_of_state_compressible.parse_parameters (prm);
+                }
+                prm.leave_subsection();
+              }
+              prm.leave_subsection();
+            }
 
           // Make options file for parsing maps to double arrays
           std::vector<std::string> chemical_field_names = this->introspection().chemical_composition_field_names();