Merge pull request #5996 from gassmoeller/simplify_use_of_phase_function

Simplify use of phase function for computing volume fractions

include/aspect/material_model/utilities.h

@@ -592,6 +592,11 @@ namespace aspect
            */
           unsigned int n_phases () const;
 
+          /**
+           * Return the total number of phases over all chemical compositions.
+           */
+          unsigned int n_phases_over_all_chemical_compositions () const;
+
           /**
            * Return the Clapeyron slope (dp/dT of the transition) for
            * phase transition number @p phase_index.
@@ -603,14 +608,32 @@ namespace aspect
            */
           double get_transition_depth (const unsigned int phase_index) const;
 
+          /**
+           * Return how many phase transitions there are for each chemical composition.
+           */
+          const std::vector<unsigned int> &
+          n_phase_transitions_for_each_chemical_composition () const;
+
+          /**
+           * Return how many phases there are for each chemical composition.
+           */
+          const std::vector<unsigned int> &
+          n_phases_for_each_chemical_composition () const;
+
           /**
            * Return how many phase transitions there are for each composition.
+           * Note, that most likely you only need the number of phase transitions
+           * for each chemical composition, so use the function above instead.
+           * This function is only kept for backward compatibility.
            */
           const std::vector<unsigned int> &
           n_phase_transitions_for_each_composition () const;
 
           /**
            * Return how many phases there are for each composition.
+           * Note, that most likely you only need the number of phase transitions
+           * for each chemical composition, so use the function above instead.
+           * This function is only kept for backward compatibility.
            */
           const std::vector<unsigned int> &
           n_phases_for_each_composition () const;
@@ -667,10 +690,25 @@ namespace aspect
            */
           std::vector<unsigned int> n_phases_per_composition;
 
+          /**
+           * A vector that stores how many phase transitions there are for each chemical compositional field.
+           */
+          std::vector<unsigned int> n_phase_transitions_per_chemical_composition;
+
+          /**
+           * A vector that stores how many phases there are for each chemical compositional field.
+           */
+          std::vector<unsigned int> n_phases_per_chemical_composition;
+
           /**
            * Total number of phases over all compositional fields
            */
           unsigned int n_phases_total;
+
+          /**
+           * Total number of phases over all compositional fields
+           */
+          unsigned int n_phases_total_chemical_compositions;
       };
     }
   }

source/material_model/utilities.cc

@@ -1296,6 +1296,8 @@ namespace aspect
           return transition_pressures.size();
       }
 
+
+
       template <int dim>
       unsigned int
       PhaseFunction<dim>::
@@ -1304,13 +1306,27 @@ namespace aspect
         return n_phases_total;
       }
 
+
+
+      template <int dim>
+      unsigned int
+      PhaseFunction<dim>::
+      n_phases_over_all_chemical_compositions () const
+      {
+        return n_phases_total_chemical_compositions;
+      }
+
+
+
       template <int dim>
       const std::vector<unsigned int> &
       PhaseFunction<dim>::n_phase_transitions_for_each_composition () const
       {
         return *n_phase_transitions_per_composition;
       }
 
+
+
       template <int dim>
       const std::vector<unsigned int> &
       PhaseFunction<dim>::n_phases_for_each_composition () const
@@ -1320,6 +1336,24 @@ namespace aspect
 
 
 
+      template <int dim>
+      const std::vector<unsigned int> &
+      PhaseFunction<dim>::n_phase_transitions_for_each_chemical_composition () const
+      {
+        return n_phase_transitions_per_chemical_composition;
+      }
+
+
+
+      template <int dim>
+      const std::vector<unsigned int> &
+      PhaseFunction<dim>::n_phases_for_each_chemical_composition () const
+      {
+        return n_phases_per_chemical_composition;
+      }
+
+
+
       template <int dim>
       double
       PhaseFunction<dim>::
@@ -1516,6 +1550,18 @@ namespace aspect
             n_phases_per_composition.push_back(n+1);
             n_phases_total += n+1;
           }
+
+        Assert(has_background_field == true, ExcInternalError());
+        // The background field is always the first composition
+        n_phases_per_chemical_composition = {n_phases_per_composition[0]};
+        n_phase_transitions_per_chemical_composition = {n_phases_per_composition[0] - 1};
+        n_phases_total_chemical_compositions = n_phases_per_composition[0];
+        for (auto i : this->introspection().chemical_composition_field_indices())
+          {
+            n_phases_per_chemical_composition.push_back(n_phases_per_composition[i+1]);
+            n_phase_transitions_per_chemical_composition.push_back(n_phases_per_composition[i+1] - 1);
+            n_phases_total_chemical_compositions += n_phases_per_composition[i+1];
+          }
       }
     }
   }

source/material_model/visco_plastic.cc

@@ -384,21 +384,9 @@ namespace aspect
           phase_function.initialize_simulator (this->get_simulator());
           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]};
-          n_phase_transitions_for_each_chemical_composition = {n_phases_for_each_composition[0] - 1};
-          n_phases = n_phases_for_each_composition[0];
-          for (auto i : this->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 std::vector<unsigned int> n_phases_for_each_chemical_composition = phase_function.n_phases_for_each_chemical_composition();
+          n_phase_transitions_for_each_chemical_composition = phase_function.n_phase_transitions_for_each_chemical_composition();
+          n_phases = phase_function.n_phases_over_all_chemical_compositions();
 
           // Equation of state parameters
           equation_of_state.initialize_simulator (this->get_simulator());