feat: Change molarity concentration to molality concentration

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr4108-17561-eb62184
+  baseline: integratedTests/baseline_integratedTests-pr3996-17639-718fc69
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -4,6 +4,10 @@ Notes
 This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
+PR #3996 (2026-09-20) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3996-17639-718fc69.tar.gz>
+=====================
+Switch the single-phase reactive transport solver to a molality-based primary
+species concentration by introducing solventMassFraction.
 
 PR #3994 (2026-09-03) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3994-17525-4ae3593.tar.gz>
 =====================

inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml

@@ -7,7 +7,7 @@
       name="SinglePhaseReactiveFlow"
       discretization="singlePhaseTPFA"
       targetRegions="{ Medium, Source }"
-      initialDt="1e-5"
+      initialDt="1e-10"
       writeLinearSystem="0"
       logLevel="1">
       <NonlinearSolverParameters
@@ -18,7 +18,7 @@
         lineSearchAction="None"
         timeStepIncreaseFactor="2"
         logLevel="1"
-        maxAllowedResidualNorm="1e100"/>
+        maxAllowedResidualNorm="1e10"/>
       <LinearSolverParameters
         directParallel="1"/>
     </SinglePhaseReactiveTransport>
@@ -37,14 +37,16 @@
   </ElementRegions>
 
   <Constitutive>
+    <!-- synthetic benchmark: solventMassFraction * defaultDensity = 1, so molality equals concentration per m^3 -->
     <ReactiveCompressibleSinglePhaseFluid
       name="water"
       defaultDensity="1000"
       defaultViscosity="0.001"
       referencePressure="0.0"
       compressibility="0.0"
       viscosibility="0.0"
-      chemicalSystemType="chainSerialAllKinetic"/>
+      chemicalSystemType="chainSerialAllKinetic"
+      solventMassFraction="0.001"/>
 
     <CompressibleSolidConstantPermeability
       name="rock"

inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml

@@ -5,6 +5,7 @@
     <File name="./ChainSerial_base.xml" />
     <File name="./ChainSerial_1DAllKinetic_base.xml" />
     <File name="./ChainSerial_initialAggregate.xml" />
+    <File name="./ChainSerial_initialLogConc.xml" />
     <File name="./ChainSerial_logConc.xml" />
   </Included>
 
@@ -25,17 +26,17 @@
     maxTime="400">
     <PeriodicEvent
       name="outputs"
-      cycleFrequency="5"
+      timeFrequency="20"
       target="/Outputs/vtkOutput"/>
 
     <PeriodicEvent
       name="logPrimarySpeciesHistoryCollection"
-      cycleFrequency="5"
+      timeFrequency="1"
       target="/Tasks/logPrimarySpeciesCollection"/>
 
     <PeriodicEvent
       name="logPrimarySpeciesTimeHistoryOutput"
-      cycleFrequency="1"
+      timeFrequency="1"
       target="/Outputs/logPrimarySpeciesHistoryOutput"/>
 
     <PeriodicEvent
@@ -45,7 +46,7 @@
 
     <PeriodicEvent
       name="restarts"
-      cycleFrequency="500"
+      timeFrequency="100"
       target="/Outputs/restartOutput"/>
   </Events>
 

inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml

@@ -0,0 +1,43 @@
+<?xml version="1.0" ?>
+
+<Problem>
+  <Included>
+    <File name="./ChainSerial_base.xml" />
+    <File name="./ChainSerial_1DAllKinetic_base.xml" />
+    <File name="./ChainSerial_initialAggregate.xml" />
+    <File name="./ChainSerial_initialLogConc.xml" />
+    <File name="./ChainSerial_logConc.xml" />
+  </Included>
+
+  <Mesh>
+    <InternalMesh
+      name="mesh1"
+      elementTypes="{ C3D8 }"
+      xCoords="{ -1e-4, 0, 40.0 }"
+      yCoords="{ 0, 1 }"
+      zCoords="{ 0, 1 }"
+      nx="{ 1, 40 }"
+      ny="{ 1 }"
+      nz="{ 1 }"
+      cellBlockNames="{ cb-0_0_0, cb-1_0_0 }"/>
+  </Mesh>
+
+  <Events
+    maxTime="20">
+    <PeriodicEvent
+      name="outputs"
+      timeFrequency="5"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="solverApplications"
+      maxEventDt="1"
+      target="/Solvers/SinglePhaseReactiveFlow"/>
+
+    <PeriodicEvent
+      name="restarts"
+      cycleFrequency="20"
+      target="/Outputs/restartOutput"/>
+  </Events>
+
+</Problem>

inputFiles/singlePhaseReactiveTransport/ChainSerial_initialLogConc.xml

@@ -0,0 +1,33 @@
+<?xml version="1.0" ?>
+
+<Problem>
+  <FieldSpecifications>
+  <!-- INITIAL LOG CONCENTRATIONS AT MEDIUM, FOR SYSTEMS WITH NO EQUILIBRIUM REACTIONS -->
+    <FieldSpecification
+      name="initialLog_C1_medium"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/Medium"
+      fieldName="logPrimarySpeciesConcentration"
+      component="0"
+      scale="-18.420680744"/>
+
+    <FieldSpecification
+      name="initialLog_C2_medium"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/Medium"
+      fieldName="logPrimarySpeciesConcentration"
+      component="1"
+      scale="-18.420680744"/>
+
+    <FieldSpecification
+      name="initialLog_C3_medium"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/Medium"
+      fieldName="logPrimarySpeciesConcentration"
+      component="2"
+      scale="-18.420680744"/>
+  </FieldSpecifications>
+</Problem>

inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DAdvective_base.xml

@@ -43,14 +43,16 @@
   </ElementRegions>
 
   <Constitutive>
+    <!-- synthetic benchmark: solventMassFraction * defaultDensity = 1, so molality equals concentration per m^3 -->
     <ReactiveCompressibleSinglePhaseFluid
       name="water"
       defaultDensity="1000"
       defaultViscosity="0.001"
       referencePressure="0.0"
       compressibility="0.0"
       viscosibility="0.0"
-      chemicalSystemType="momasEasy"/>
+      chemicalSystemType="momasEasy"
+      solventMassFraction="0.001"/>
 
     <CompressibleSolidConstantPermeability
       name="rock_a"

inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml

@@ -43,14 +43,16 @@
   </ElementRegions>
 
   <Constitutive>
+    <!-- synthetic benchmark: solventMassFraction * defaultDensity = 1, so molality equals concentration per m^3 -->
     <ReactiveCompressibleSinglePhaseFluid
       name="water"
       defaultDensity="1000"
       defaultViscosity="0.001"
       referencePressure="0.0"
       compressibility="0.0"
       viscosibility="0.0"
-      chemicalSystemType="momasEasy"/>
+      chemicalSystemType="momasEasy"
+      solventMassFraction="0.001"/>
 
     <CompressibleSolidConstantPermeability
       name="rock_a"

inputFiles/singlePhaseReactiveTransport/MoMaSEasy_2DAdvective_base.xml

@@ -55,14 +55,16 @@
   </ElementRegions>
 
   <Constitutive>
+    <!-- synthetic benchmark: solventMassFraction * defaultDensity = 1, so molality equals concentration per m^3 -->
     <ReactiveCompressibleSinglePhaseFluid
       name="water"
       defaultDensity="1000"
       defaultViscosity="0.001"
       referencePressure="0.0"
       compressibility="0.0"
       viscosibility="0.0"
-      chemicalSystemType="momasEasy"/>
+      chemicalSystemType="momasEasy"
+      solventMassFraction="0.001"/>
 
     <CompressibleSolidConstantPermeability
       name="rock_a"

inputFiles/singlePhaseReactiveTransport/MoMaSMedium_1DAdvective_base.xml

@@ -43,14 +43,16 @@
   </ElementRegions>
 
   <Constitutive>
+    <!-- synthetic benchmark: solventMassFraction * defaultDensity = 1, so molality equals concentration per m^3 -->
     <ReactiveCompressibleSinglePhaseFluid
       name="water"
       defaultDensity="1000"
       defaultViscosity="0.001"
       referencePressure="0.0"
       compressibility="0.0"
       viscosibility="0.0"
-      chemicalSystemType="momasMedium"/>
+      chemicalSystemType="momasMedium"
+      solventMassFraction="0.001"/>
 
     <CompressibleSolidConstantPermeability
       name="rock_a"

inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats

@@ -23,6 +23,15 @@ decks = [
         check_step=20,
         restartcheck_params=RestartcheckParameters(**restartcheck_params)),
 
+    TestDeck(
+        name="ChainSerial_1DAllKinetic_smoke",
+        description=
+        '1D serial chain reaction, all kinetic with no equilibrium reactions (1D, kineticReaction)',
+        partitions=((1, 1, 1), (2, 1, 1)),
+        restart_step=20,
+        check_step=20,
+        restartcheck_params=RestartcheckParameters(**restartcheck_params)),
+
     TestDeck(
         name="mixedReactionCarbonateSystem_1D_smoke",
         description=

src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp

@@ -55,7 +55,7 @@ DECLARE_FIELD( primarySpeciesAggregateConcentration_n,
                "primarySpeciesAggregateConcentration_n",
                array3dLayoutSpecies,
                1e-16,
-               LEVEL_0,
+               NOPLOT,
                WRITE_AND_READ,
                "primarySpeciesAggregateConcentration at the previous timestep" );
 
@@ -71,17 +71,17 @@ DECLARE_FIELD( dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentra
                "dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations",
                array4dLayoutSpecies_dC,
                0,
-               LEVEL_0,
+               NOPLOT,
                WRITE_AND_READ,
-               "Deivatives of primarySpeciesAggregateConcentration w.r.t log primary species concentration" );
+               "Derivatives of primarySpeciesAggregateConcentration w.r.t log primary species concentration" );
 
 DECLARE_FIELD( dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations,
                "dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations",
                array4dLayoutSpecies_dC,
                0,
-               LEVEL_0,
+               NOPLOT,
                WRITE_AND_READ,
-               "Deivatives of primarySpeciesMobileAggregateConcentration w.r.t log primary species concentration" );
+               "Derivatives of primarySpeciesMobileAggregateConcentration w.r.t log primary species concentration" );
 
 DECLARE_FIELD( secondarySpeciesConcentration,
                "secondarySpeciesConcentration",
@@ -111,9 +111,9 @@ DECLARE_FIELD( dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations,
                "dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations",
                array4dLayoutSpecies_dC,
                0,
-               LEVEL_0,
+               NOPLOT,
                WRITE_AND_READ,
-               "Deivatives of aggregate concentration rates w.r.t log primary species concentration" );
+               "Derivatives of aggregate concentration rates w.r.t log primary species concentration" );
 }
 
 }

src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp

@@ -30,8 +30,6 @@ namespace constitutive
 namespace reactivefluid
 {
 
-using namespace hpcReact::bulkGeneric;
-
 template< typename BASE >
 ReactiveSinglePhaseFluid< BASE >::
 ReactiveSinglePhaseFluid( string const & name, Group * const parent ):
@@ -43,6 +41,13 @@ ReactiveSinglePhaseFluid( string const & name, Group * const parent ):
     setDescription( "Chemical System type. Available options are: "
                     "``" + EnumStrings< ChemicalSystemType >::concat( "|" ) + "``" );
 
+  this->registerWrapper( viewKeyStruct::solventMassFractionString(), &m_solventMassFraction ).
+    setApplyDefaultValue( 1.0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "Mass fraction of solvent in the solution [-]. Species molality [mol/kg solvent] "
+                    "times this fraction is the amount per kg of solution, and times the fluid density "
+                    "the amount per m^3 of solution. The default of 1 is a dilute aqueous solution." );
+
   this->template registerField< fields::reactivefluid::initialPrimarySpeciesConcentration >( &m_initialPrimarySpeciesConcentration );
   this->template registerField< fields::reactivefluid::secondarySpeciesConcentration >( &m_secondarySpeciesConcentration );
   this->template registerField< fields::reactivefluid::primarySpeciesAggregateConcentration >( &m_primarySpeciesAggregateConcentration );
@@ -69,6 +74,7 @@ deliverClone( string const & name, Group * const parent ) const
   newConstitutiveRelation.m_numPrimarySpecies = m_numPrimarySpecies;
   newConstitutiveRelation.m_numSecondarySpecies = m_numSecondarySpecies;
   newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions;
+  newConstitutiveRelation.m_solventMassFraction = m_solventMassFraction;
 
   return clone;
 }
@@ -116,6 +122,11 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization()
       m_numKineticReactions = 0;
       break;
   }
+
+  GEOS_THROW_IF( m_solventMassFraction <= 0.0 || m_solventMassFraction > 1.0,
+                 GEOS_FMT( "{}: invalid value of attribute '{}', expected a value in (0, 1]",
+                           this->getDataContext(), viewKeyStruct::solventMassFractionString() ),
+                 InputError );
 }
 
 template< typename BASE >

src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp

@@ -113,6 +113,15 @@ class ReactiveSinglePhaseFluid : public BASE
 
   integer numKineticReactions() const { return m_numKineticReactions; }
 
+  /**
+   * @brief Mass fraction of solvent in the solution [-].
+   *
+   * HPCReact is a molality based library: concentrations, equilibrium constants and mass-action
+   * quotients are all on the molal scale [mol/kg solvent]. With w this fraction and rho the solution
+   * density, moles per kg of solution are m w and moles per m^3 of solution are m w rho.
+   */
+  real64 solventMassFraction() const { return m_solventMassFraction; }
+
   /**
    * @brief Kernel wrapper class for ReactiveSinglePhaseFluid.
    */
@@ -235,7 +244,6 @@ class ReactiveSinglePhaseFluid : public BASE
   {
     using namespace hpcReact::geochemistry;
     using namespace hpcReact::MoMasBenchmark;
-    using namespace hpcReact::bulkGeneric;
     using namespace hpcReact::ChainGeneric;
     switch( m_chemicalSystemType )
     {
@@ -330,6 +338,7 @@ class ReactiveSinglePhaseFluid : public BASE
   struct viewKeyStruct : ConstitutiveBase::viewKeyStruct
   {
     static constexpr char const * chemicalSystemNameString() { return "chemicalSystemType"; }
+    static constexpr char const * solventMassFractionString() { return "solventMassFraction"; }
   };
 
 protected:
@@ -365,6 +374,10 @@ class ReactiveSinglePhaseFluid : public BASE
   array4d< real64, constitutive::reactivefluid::LAYOUT_SPECIES_DC >  m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations;
 
   ChemicalSystemType m_chemicalSystemType;
+
+  /// Mass fraction of solvent in the solution [-], constant; for the carbonate brine EQ3/6 gives 0.898
+  /// TODO: depends on concentration, w = 1 / (1 + sum_i m_i M_i)
+  real64 m_solventMassFraction;
 };
 
 // these aliases are useful in constitutive dispatch