Merge pull request #6366 from danieldouglas92/correctly_compute_mass_fractions_tian

Correctly Compute Reactions in the Reactive Fluid Transport Model

Author: jdannberg

benchmarks/solubility/plugin/solubility.cc

@@ -72,9 +72,14 @@ namespace aspect
          * @param in Object that contains the current conditions.
          * @param melt_fractions Vector of doubles that is filled with the
          * allowable free fluid fraction for each given input conditions.
+         * @param out Optional pointer to the material properties provided by the
+         * material model. By default, this variable is a nullptr. If the melt
+         * fractions depend on material model properties, then this parameter
+         * must be set to a valid pointer to a MaterialModelOutputs object.
          */
         void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                             std::vector<double> &melt_fractions) const override;
+                             std::vector<double> &melt_fractions,
+                             const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const override;
 
         double reference_darcy_coefficient () const override;
 
@@ -357,7 +362,8 @@ namespace aspect
     void
     Volatiles<dim>::
     melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                    std::vector<double> &melt_fractions) const
+                    std::vector<double> &melt_fractions,
+                    const MaterialModel::MaterialModelOutputs<dim> *) const
     {
       for (unsigned int q=0; q<in.n_evaluation_points(); ++q)
         {

cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm

@@ -61,6 +61,10 @@ subsection Solver parameters
     set Linear solver tolerance             = 1e-10
     set Number of cheap Stokes solver steps = 1000
   end
+  subsection Operator splitting parameters
+    set Reaction solver type = fixed step
+    set Reaction time step   = 10
+  end
 end
 
 # Initialize the compositional fields. When using the tian approximation

doc/modules/changes/20250614_danieldouglas92

@@ -0,0 +1,5 @@
+Fixed: The conservation of water mass when being partitioned
+within the solid phase and the free fluid phase in the
+reactive fluid transport model.
+<br>
+(Daniel Douglas, 2025/06/14)

include/aspect/material_model/latent_heat_melt.h

@@ -73,7 +73,8 @@ namespace aspect
          */
 
         void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                             std::vector<double> &melt_fractions) const override;
+                             std::vector<double> &melt_fractions,
+                             const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const override;
 
 
         /**

include/aspect/material_model/melt_boukare.h

@@ -90,9 +90,14 @@ namespace aspect
          * @param in Object that contains the current conditions.
          * @param melt_fractions Vector of doubles that is filled with the
          * equilibrium melt fraction for each given input conditions.
+         * @param out Optional pointer to the material properties provided by the
+         * material model. By default, this variable is a nullptr. If the melt
+         * fractions depend on material model properties, then this parameter
+         * must be set to a valid pointer to a MaterialModelOutputs object.
          */
         void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                             std::vector<double> &melt_fractions) const override;
+                             std::vector<double> &melt_fractions,
+                             const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const override;
 
         /**
          * @name Reference quantities

include/aspect/material_model/melt_global.h

@@ -67,9 +67,14 @@ namespace aspect
          * @param in Object that contains the current conditions.
          * @param melt_fractions Vector of doubles that is filled with the
          * equilibrium melt fraction for each given input conditions.
+         * @param out Optional pointer to the material properties provided by the
+         * material model. By default, this variable is a nullptr. If the melt
+         * fractions depend on material model properties, then this parameter
+         * must be set to a valid pointer to a MaterialModelOutputs object.
          */
         void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                             std::vector<double> &melt_fractions) const override;
+                             std::vector<double> &melt_fractions,
+                             const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const override;
 
         /**
          * @name Reference quantities

include/aspect/material_model/melt_simple.h

@@ -78,7 +78,8 @@ namespace aspect
                       typename Interface<dim>::MaterialModelOutputs &out) const override;
 
         void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                             std::vector<double> &melt_fractions) const override;
+                             std::vector<double> &melt_fractions,
+                             const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const override;
 
         /**
          * @name Reference quantities

include/aspect/material_model/reaction_model/tian2019_solubility.h

@@ -61,12 +61,12 @@ namespace aspect
            * @p in and @p melt_fraction need to have the same size.
            *
            * @param in Object that contains the current conditions.
-           * @param porosity_idx the index of the "porosity" composition
+           * @param mass_frac_porosity the mass fraction of the "porosity" composition
            * @param q the quadrature point index
            */
           double
           melt_fraction (const MaterialModel::MaterialModelInputs<dim> &in,
-                         const unsigned int porosity_idx,
+                         const double mass_frac_porosity,
                          unsigned int q) const;
 
           /**

include/aspect/material_model/reactive_fluid_transport.h

@@ -68,9 +68,14 @@ namespace aspect
          * @param in Object that contains the current conditions.
          * @param melt_fractions Vector of doubles that is filled with the
          * allowable free fluid fraction for each given input conditions.
+         * @param out Optional pointer to the material properties provided by the
+         * material model. By default, this variable is a nullptr. If the melt
+         * fractions depend on material model properties, then this parameter
+         * must be set to a valid pointer to a MaterialModelOutputs object.
          */
         void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                             std::vector<double> &melt_fractions) const override;
+                             std::vector<double> &melt_fractions,
+                             const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const override;
 
         /**
          * Initialize the base model at the beginning of the model run

include/aspect/melt.h

@@ -153,9 +153,14 @@ namespace aspect
          * @param in Object that contains the current conditions.
          * @param melt_fractions Vector of doubles that is filled with the
          * equilibrium melt fraction for each given input conditions.
+         * @param out Optional pointer to the material properties provided by the
+         * material model. By default, this variable is a nullptr. If the melt
+         * fractions depend on material model properties, then this parameter
+         * must be set to a valid pointer to a MaterialModelOutputs object.
          */
         virtual void melt_fractions (const MaterialModel::MaterialModelInputs<dim> &in,
-                                     std::vector<double> &melt_fractions) const = 0;
+                                     std::vector<double> &melt_fractions,
+                                     const MaterialModel::MaterialModelOutputs<dim> *out = nullptr) const = 0;
 
         /**
          * Return whether an object provided as argument is of a class that is