Merge branch 'MolarVolume' into 'master'

[CL] Set reactant molality alternatively

See merge request !3372

ChemistryLib/PhreeqcIO.cpp

@@ -70,17 +70,48 @@ void setAqueousSolution(std::vector<double> const& concentrations,
 template <typename Reactant>
 void setReactantMolality(Reactant& reactant,
                          GlobalIndexType const& chemical_system_id,
-                         MaterialPropertyLib::Phase const& solid_phase,
+                         MaterialPropertyLib::Medium const* medium,
                          ParameterLib::SpatialPosition const& pos,
                          double const t)
 {
+    auto const& solid_phase = medium->phase("Solid");
     auto const& solid_constituent = solid_phase.component(reactant.name);
 
-    auto const molality =
-        solid_constituent.property(MaterialPropertyLib::PropertyType::molality)
-            .template initialValue<double>(pos, t);
+    auto const& liquid_phase = medium->phase("AqueousLiquid");
+
+    if (solid_constituent.hasProperty(
+            MaterialPropertyLib::PropertyType::molality))
+    {
+        auto const molality =
+            solid_constituent
+                .property(MaterialPropertyLib::PropertyType::molality)
+                .template initialValue<double>(pos, t);
 
-    (*reactant.molality)[chemical_system_id] = molality;
+        (*reactant.molality)[chemical_system_id] = molality;
+    }
+    else
+    {
+        auto const volume_fraction =
+            solid_constituent
+                .property(MaterialPropertyLib::PropertyType::volume_fraction)
+                .template initialValue<double>(pos, t);
+
+        auto const fluid_density =
+            liquid_phase.property(MaterialPropertyLib::PropertyType::density)
+                .template initialValue<double>(pos, t);
+
+        auto const porosity =
+            medium->property(MaterialPropertyLib::PropertyType::porosity)
+                .template initialValue<double>(pos, t);
+
+        auto const molar_volume =
+            solid_constituent
+                .property(MaterialPropertyLib::PropertyType::molar_volume)
+                .template initialValue<double>(pos, t);
+
+        (*reactant.molality)[chemical_system_id] =
+            volume_fraction / fluid_density / porosity / molar_volume;
+    }
 }
 
 template <typename Reactant>
@@ -171,17 +202,16 @@ void PhreeqcIO::initializeChemicalSystemConcrete(
     setAqueousSolution(concentrations, chemical_system_id,
                        *_chemical_system->aqueous_solution);
 
-    auto const& solid_phase = medium->phase("Solid");
     for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
     {
-        setReactantMolality(kinetic_reactant, chemical_system_id, solid_phase,
-                            pos, t);
+        setReactantMolality(kinetic_reactant, chemical_system_id, medium, pos,
+                            t);
     }
 
     for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
     {
-        setReactantMolality(equilibrium_reactant, chemical_system_id,
-                            solid_phase, pos, t);
+        setReactantMolality(equilibrium_reactant, chemical_system_id, medium,
+                            pos, t);
     }
 }
 

MaterialLib/MPL/PropertyType.h

@@ -59,6 +59,7 @@ enum PropertyType : int
     longitudinal_dispersivity,
     molality,
     molar_mass,
+    molar_volume,
     mole_fraction,
     /// used to compute the hydrodynamic dispersion tensor.
     molecular_diffusion,
@@ -88,6 +89,7 @@ enum PropertyType : int
     transversal_dispersivity,
     vapor_pressure,
     viscosity,
+    volume_fraction,
     number_of_properties
 };
 
@@ -182,7 +184,6 @@ inline PropertyType convertStringToProperty(std::string const& inString)
     {
         return PropertyType::longitudinal_dispersivity;
     }
-    // TODO (renchao): add property "volume fraction"
     if (boost::iequals(inString, "molality"))
     {
         return PropertyType::molality;
@@ -191,6 +192,10 @@ inline PropertyType convertStringToProperty(std::string const& inString)
     {
         return PropertyType::molar_mass;
     }
+    if (boost::iequals(inString, "molar_volume"))
+    {
+        return PropertyType::molar_volume;
+    }
     if (boost::iequals(inString, "mole_fraction"))
     {
         return PropertyType::mole_fraction;
@@ -295,6 +300,10 @@ inline PropertyType convertStringToProperty(std::string const& inString)
     {
         return PropertyType::viscosity;
     }
+    if (boost::iequals(inString, "volume_fraction"))
+    {
+        return PropertyType::volume_fraction;
+    }
 
     OGS_FATAL(
         "The property name '{:s}' does not correspond to any known property",
@@ -328,6 +337,7 @@ static const std::array<std::string, PropertyType::number_of_properties>
                              "longitudinal_dispersivity",
                              "molality",
                              "molar_mass",
+                             "molar_volume",
                              "mole_fraction",
                              "molecular_diffusion",
                              "name",
@@ -353,7 +363,8 @@ static const std::array<std::string, PropertyType::number_of_properties>
                              "transport_porosity",
                              "transversal_dispersivity",
                              "vapor_pressure",
-                             "viscosity"}};
+                             "viscosity",
+                             "volume_fraction"}};
 
 /// This data type is based on a std::array. It can hold pointers to objects of
 /// class Property or its inheritors. The size of this array is determined by

Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/KineticReactant/1d_isofrac_flag_formula.prj

@@ -111,9 +111,14 @@
                             <name>Productd</name>
                             <properties>
                                 <property>
-                                    <name>molality</name>
+                                    <name>volume_fraction</name>
                                     <type>Constant</type>
-                                    <value>1e-6</value>
+                                    <value>0.1</value>
+                                </property>
+                                <property>
+                                    <name>molar_volume</name>
+                                    <type>Constant</type>
+                                    <value>1e2</value>
                                 </property>
                             </properties>
                         </component>