Merge pull request #6131 from lhy11009/discrete_phase_funcion_2

Discrete phase function

data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt

@@ -0,0 +1,8 @@
+# This is a data output from HeFESTo
+# Independent variables are T(K) and P(bar)
+# POINTS: 2 2
+T(K)                P(bar)              s,J/K/kg            rho,kg/m3           alpha,1/K           cp,J/K/kg           vp,km/s             vs,km/s             h,J/kg  dominant_phase_index,1
+250.0               0.0                 509.74740059944844  3332.8967443100287  1.9781287327429287e-05 748.5636509347647   8.116376883289359   4.7495273882440685  -13445454.390915 0.0
+4000.0              0.0                 3804.799798825367   2879.31713365359    0.00012376961201237974 1865.6003296108431  5.89690417370827    3.068046128481904   -8401781.805825338 0.0 
+250.0               400781.25           326.1989261514726   4738.970221094591   9.15389351114028e-06 565.4531513648064   11.974450477467947  6.779618381631337   -3722250.1404551915 1.0                
+4000.0              400781.25           3394.028743328079   4373.768154238075   2.6303932404721963e-05 1362.8625515701538  10.986106548987486  5.8231315598976465  873079.6281181354   1.0                

doc/modules/changes/20241107_lhy11009

@@ -0,0 +1,6 @@
+Added: there is now a new class of phase function that handles discrete phase transitions 
+by looking up the most dominant phases in a lookup table. This function can be used
+to make the rheology of the visco-plastic material model dependent on the dominant
+mineral phase.
+<br>
+(Haoyuan Li, 2024/11/07)

include/aspect/material_model/utilities.h

@@ -549,6 +549,191 @@ namespace aspect
         unsigned int phase_index;
       };
 
+      /**
+       * A class that bundles functionality to look up the dominant phase in
+       * tables for each respective composition and export the values
+       * of phase functions. The class can handle arbitrary numbers of
+       * dominant phases for each composition, but the calling side
+       * has to determine how to use the return values of this object
+       * (e.g. in terms of density or viscosity).
+       */
+      template <int dim>
+      class PhaseFunctionDiscrete: public ::aspect::SimulatorAccess<dim>
+      {
+        public:
+
+          /**
+           * The initialization process loads the contents of the material files
+           * for the respective compositions.
+           */
+          void initialize();
+
+          /**
+           * Percentage of material that has already undergone the phase
+           * transition to the higher-pressure material. For this class
+           * this function only returns 1.0 or 0.0, depending on whether
+           * the selected phase transition has been crossed or not.
+           */
+          double compute_value (const PhaseFunctionInputs<dim> &in) const;
+
+          /**
+           * No valid implementation exists for this function, as the derivative of a
+           * discrete function is undefined at locations of phase jumps. This function raises an
+           * error to ensure that a phase derivative request is not made for this phase function.
+           */
+          double compute_derivative () const;
+
+          /**
+           * Return the total number of phase transitions.
+           */
+          unsigned int n_phase_transitions () const;
+
+          /**
+           * Return the total number of phases.
+           */
+          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 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;
+
+          /**
+           * Declare the parameters this class takes through input files.
+           * Note that this class does not declare its own subsection,
+           * i.e. the parameters will be declared in the subsection that
+           * was active before calling this function.
+           */
+          static
+          void
+          declare_parameters (ParameterHandler &prm);
+
+          /**
+           * Read the parameters this class declares from the parameter file.
+           * Note that this class does not declare its own subsection,
+           * i.e. the parameters will be parsed from the subsection that
+           * was active before calling this function.
+           */
+          void
+          parse_parameters (ParameterHandler &prm);
+
+
+        private:
+          /**
+          * Directory path where data files are stored.
+          */
+          std::string data_directory;
+
+          /**
+           * List of file names containing material data for each composition.
+          */
+          std::vector<std::string> material_file_names;
+
+          /**
+          * Minimum temperature values for each composition in the P-T table.
+          */
+          std::vector<double> minimum_temperature;
+
+          /**
+          * Maximum temperature values for each composition in the P-T table.
+          */
+          std::vector<double> maximum_temperature;
+
+          /**
+          * Temperature intervals used for each composition in the P-T table.
+          */
+          std::vector<double> interval_temperature;
+
+          /**
+          * Minimum pressure values for each composition in the P-T table.
+          */
+          std::vector<double> minimum_pressure;
+
+          /**
+          * Maximum pressure values for each composition in the P-T table.
+          */
+          std::vector<double> maximum_pressure;
+
+          /**
+          * Pressure intervals used for each composition in the P-T table.
+          */
+          std::vector<double> interval_pressure;
+
+          /**
+           * List of pointers to objects that read and process data we get from
+           * material data files. There is one pointer/object per lookup file.
+           */
+          std::vector<std::unique_ptr<Utilities::StructuredDataLookup<2>>> material_lookup;
+
+          /**
+           * List of phase indicators of the most dominant phases in the material data files
+           * to construct the different phase transitions in this class. For a description of
+           * the use of the phase indicators, please see the documentation of the input parameter
+           * 'Phase transition indicators' in the function declare_parameters().
+           */
+          std::vector<unsigned int> transition_indicators;
+
+          /**
+           * A vector that stores how many phase transitions there are for each compositional field.
+           */
+          std::unique_ptr<std::vector<unsigned int>> n_phase_transitions_per_composition;
+
+          /**
+           * A vector that stores how many phases there are for each compositional field.
+           */
+          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;
+      };
+
       /**
        * A class that bundles functionality to compute the values and
        * derivatives of phase functions. The class can handle arbitrary

include/aspect/material_model/visco_plastic.h

@@ -181,6 +181,12 @@ namespace aspect
     class ViscoPlastic : public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
     {
       public:
+        /**
+         * Initialization function. Loads the material data and sets up
+         * pointers if it is required.
+         */
+        void
+        initialize () override;
 
         void evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
                       MaterialModel::MaterialModelOutputs<dim> &out) const override;
@@ -261,6 +267,16 @@ namespace aspect
          */
         MaterialUtilities::PhaseFunction<dim> phase_function;
 
+        /**
+         * Determines whether to look up the dominant phases for each composition in its respective lookup table.
+         */
+        bool use_dominant_phase_for_viscosity;
+
+        /**
+         * Object that handles discrete phase transitions for the rheology if requested by the variable use_dominant_phase_for_viscosity.
+         */
+        std::unique_ptr<MaterialUtilities::PhaseFunctionDiscrete<dim>> phase_function_discrete;
+
     };
 
   }

include/aspect/structured_data.h

@@ -248,6 +248,15 @@ namespace aspect
          */
         double get_maximum_component_value(const unsigned int component) const;
 
+        /**
+         * Retrieve the number of table points for a given dimension.
+         * Equivalent to calling get_interpolation_point_coordinates().size().
+         *
+         * @param dimension The index of the dimension for which to get the number of table points.
+         * @return The number of points along the specified dimension.
+         */
+        unsigned int get_number_of_coordinates(const unsigned int dimension) const;
+
       private:
         /**
          * The number of data components read in (=columns in the data file).