Skip to content

feat: Change molarity concentration to molality concentration #3996

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
frankfeifan wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

feat: Change molarity concentration to molality concentration #3996

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0a5f798
switch from molarity concentration to molality
frankfeifan Mar 14, 2026
befd8a9
make the timestep in unit test smaller as the jacobian is very sensit…
frankfeifan Mar 16, 2026
243c99c
Merge remote-tracking branch 'origin/develop' into feat/frankfei/mola…
frankfeifan Mar 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/coreComponents/constitutive/HPCReact
Open in desktop
Submodule HPCReact updated 5 files
+1 −0 src/common/constants.hpp
+6 −2 src/reactions/exampleSystems/BulkGeneric.hpp
+2 −1 src/reactions/exampleSystems/ChainGeneric.hpp
+14 −2 src/reactions/exampleSystems/MoMasBenchmark.hpp
+7 −2 src/reactions/reactionsSystems/Parameters.hpp
10 changes: 10 additions & 0 deletions src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,9 @@ namespace reactivefluid
{

using namespace hpcReact::bulkGeneric;
using namespace hpcReact::geochemistry;
using namespace hpcReact::ChainGeneric;
using namespace hpcReact::MoMasBenchmark;

template< typename BASE >
ReactiveSinglePhaseFluid< BASE >::
Expand Down Expand Up @@ -69,6 +72,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_solventDensity = m_solventDensity;

return clone;
}
Expand All @@ -84,36 +88,42 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization()
m_numPrimarySpecies = 9;
m_numSecondarySpecies = 16;
m_numKineticReactions = 5;
m_solventDensity = ultramaficSystem.getSolventDensity();
break;

case ChemicalSystemType::carbonate:
m_numPrimarySpecies = 7;
m_numSecondarySpecies = 10;
m_numKineticReactions = 1;
m_solventDensity = carbonateSystem.getSolventDensity();
break;

case ChemicalSystemType::carbonateAllEquilibrium:
m_numPrimarySpecies = 7;
m_numSecondarySpecies = 11;
m_numKineticReactions = 0;
m_solventDensity = carbonateSystemAllEquilibrium.getSolventDensity();
break;

case ChemicalSystemType::chainSerialAllKinetic:
m_numPrimarySpecies = 3;
m_numSecondarySpecies = 0;
m_numKineticReactions = 3;
m_solventDensity = serialAllKineticParams.getSolventDensity();
break;

case ChemicalSystemType::momasMedium:
m_numPrimarySpecies = 5;
m_numSecondarySpecies = 9;
m_numKineticReactions = 1;
m_solventDensity = mediumCaseParams.getSolventDensity();
break;

default:
m_numPrimarySpecies = 5;
m_numSecondarySpecies = 7;
m_numKineticReactions = 0;
m_solventDensity = easyCaseParams.getSolventDensity();
break;
}
}
Expand Down
4 changes: 4 additions & 0 deletions src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,8 @@ class ReactiveSinglePhaseFluid : public BASE

integer numKineticReactions() const { return m_numKineticReactions; }

real64 solventDensity() const { return m_solventDensity; }

/**
* @brief Kernel wrapper class for ReactiveSinglePhaseFluid.
*/
Expand Down Expand Up @@ -365,6 +367,8 @@ class ReactiveSinglePhaseFluid : public BASE
array4d< real64, constitutive::reactivefluid::LAYOUT_SPECIES_DC > m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations;

ChemicalSystemType m_chemicalSystemType;

real64 m_solventDensity;
};

// these aliases are useful in constitutive dispatch
Expand Down
2 changes: 1 addition & 1 deletion src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -530,7 +530,7 @@ class SinglePhaseReactiveTransportTest : public ::testing::Test
}

static real64 constexpr time = 0.0;
static real64 constexpr dt = 1.0;
static real64 constexpr dt = 0.01;
static real64 constexpr eps = std::numeric_limits< real64 >::epsilon();

GeosxState state;
Expand Down
39 changes: 32 additions & 7 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -463,14 +463,33 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt,

forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
MeshLevel const & mesh,
string_array const & )
string_array const & regionNames )
{
NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName );

string const & dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );

real64 solventDensity = 1.0;
mesh.getElemManager().forElementSubRegions( regionNames,
[&]( localIndex const,
ElementSubRegionBase const & subRegion )
{
if( m_isThermal )
{
reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid const & fluid =
getConstitutiveModel< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) );
solventDensity = fluid.solventDensity();
}
else
{
reactivefluid::ReactiveCompressibleSinglePhaseFluid const & fluid =
getConstitutiveModel< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) );
solventDensity = fluid.solventDensity();
}
} );

fluxApprox.forAllStencils( mesh, [&] ( auto & stencil )
{
typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper();
Expand All @@ -481,6 +500,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt,
FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies,
m_hasDiffusion,
mobilePrimarySpeciesFlags.toViewConst(),
solventDensity,
dofManager.rankOffset(),
dofKey,
getName(),
Expand All @@ -496,6 +516,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt,
FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies,
m_hasDiffusion,
mobilePrimarySpeciesFlags.toViewConst(),
solventDensity,
dofManager.rankOffset(),
dofKey,
getName(),
Expand Down Expand Up @@ -552,15 +573,16 @@ void SinglePhaseReactiveTransport::updateSpeciesAmount( ElementSubRegionBase & s
getConstitutiveModel< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration = fluid.primarySpeciesAggregateConcentration();
arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration_n = fluid.primarySpeciesAggregateConcentration_n();
real64 const solventDensity = fluid.solventDensity();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
for( integer is = 0; is < numPrimarySpecies; ++is )
{
primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is];
primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventDensity;

if( isZero( primarySpeciesAggregateMole_n[ei][is] ) )
primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is];
primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventDensity;
}
} );
}
Expand All @@ -570,15 +592,16 @@ void SinglePhaseReactiveTransport::updateSpeciesAmount( ElementSubRegionBase & s
getConstitutiveModel< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration = fluid.primarySpeciesAggregateConcentration();
arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration_n = fluid.primarySpeciesAggregateConcentration_n();
real64 const solventDensity = fluid.solventDensity();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
for( integer is = 0; is < numPrimarySpecies; ++is )
{
primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is];
primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventDensity;

if( isZero( primarySpeciesAggregateMole_n[ei][is] ) )
primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is];
primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventDensity;
}
} );
}
Expand All @@ -598,12 +621,13 @@ void SinglePhaseReactiveTransport::updateKineticReactionMolarIncrements( real64
reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid & fluid =
getConstitutiveModel< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView3d< real64 const, reactivefluid::USD_SPECIES > const kineticReactionRates = fluid.kineticReactionRates();
real64 const solventDensity = fluid.solventDensity();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
for( integer r = 0; r < numKineticReactions; ++r )
{
kineticReactionMolarIncrements[ei][r] = dt* kineticReactionRates[ei][0][r];
kineticReactionMolarIncrements[ei][r] = dt * kineticReactionRates[ei][0][r] * solventDensity;
}
} );
}
Expand All @@ -612,12 +636,13 @@ void SinglePhaseReactiveTransport::updateKineticReactionMolarIncrements( real64
reactivefluid::ReactiveCompressibleSinglePhaseFluid & fluid =
getConstitutiveModel< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
arrayView3d< real64 const, reactivefluid::USD_SPECIES > const kineticReactionRates = fluid.kineticReactionRates();
real64 const solventDensity = fluid.solventDensity();

forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
{
for( integer r = 0; r < numKineticReactions; ++r )
{
kineticReactionMolarIncrements[ei][r] = dt* kineticReactionRates[ei][0][r];
kineticReactionMolarIncrements[ei][r] = dt * kineticReactionRates[ei][0][r] * solventDensity;
}
} );
}
Expand Down
6 changes: 3 additions & 3 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,23 +41,23 @@ DECLARE_FIELD( logPrimarySpeciesConcentration,
0,
LEVEL_0,
WRITE_AND_READ,
"Natural log of primary species concentration (molarity)" );
"Natural log of primary species concentration (molality)" );

DECLARE_FIELD( logPrimarySpeciesConcentration_n,
"logPrimarySpeciesConcentration_n",
array2dLayoutComp,
0,
LEVEL_0,
WRITE_AND_READ,
"Natural log of primary species concentration (molarity) at the previous converged time step" );
"Natural log of primary species concentration (molality) at the previous converged time step" );

DECLARE_FIELD( bcLogPrimarySpeciesConcentration,
"bcLogPrimarySpeciesConcentration",
array2dLayoutComp,
0,
LEVEL_0,
WRITE_AND_READ,
"Boundary condition for natural log of primary species concentration (molarity)" );
"Boundary condition for natural log of primary species concentration (molality)" );

DECLARE_FIELD( primarySpeciesAggregateMole,
"primarySpeciesAggregateMole",
Expand Down
26 changes: 12 additions & 14 deletions ...eComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
arrayView1d< real64 > const & localRhs )
: Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ),
m_dt( dt ),
m_solventDensity( fluid.solventDensity() ),
m_volume( subRegion.getElementVolume() ),
m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ),
m_porosity( solid.getPorosity() ),
Expand Down Expand Up @@ -160,8 +161,8 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
void computeAccumulation( localIndex const ei,
StackVariables & stack ) const
{
// Residual[is] += (primarySpeciesAggregateConcentration[is] * stack.poreVolume - primarySpeciesAggregateMole_n[is])
// - dt * m_volume * primarySpeciesKineticRate[is] // To Check: what's the unit of the kinetic rate
// Residual[is] += (primarySpeciesAggregateConcentration[is] * solventDensity * stack.poreVolume - primarySpeciesAggregateMole_n[is])
// - dt * m_volume * primarySpeciesKineticRate[is] * solventDensity

Base::computeAccumulation( ei, stack );

Expand All @@ -183,14 +184,14 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
// Step 2.1: residual
// Primary species mole amount in pore volume
stack.localResidual[is+numEqn-numSpecies] -= m_primarySpeciesAggregateMole_n[ei][is];
stack.localResidual[is+numEqn-numSpecies] += m_primarySpeciesAggregateConcentration[ei][0][is] * stack.poreVolume;
stack.localResidual[is+numEqn-numSpecies] += m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity * stack.poreVolume;

// Reaction term
stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * m_primarySpeciesAggregateKineticRate[ei][0][is];
stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * m_primarySpeciesAggregateKineticRate[ei][0][is] * m_solventDensity;

// Step 2.1: jacobian
// Drivative of primary species amount in pore volume wrt pressure
stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * m_primarySpeciesAggregateConcentration[ei][0][is]
stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity
/* + stack.poreVolume * m_dTotalPrimarySpeciesConcentration_dPres[ei][is] */;
// // Derivative of reaction term wrt pressure
// stack.localJacobian[is+numEqn-numSpecies][0] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) *
Expand All @@ -201,15 +202,9 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
{
stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] = /* stack.dPoreVolume_dLogPrimaryConc[js] *
m_primarySpeciesAggregateConcentration[ei][0][is]
+ */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[is][js]; // To
// check
// if
// the
// permutation
// is
// consistent

stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js];
+ */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[is][js] * m_solventDensity;

stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js] * m_solventDensity;
}
}
}
Expand Down Expand Up @@ -245,6 +240,9 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
/// Time step size
real64 const m_dt;

/// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³]
real64 const m_solventDensity;

/// View on the element volumes
arrayView1d< real64 const > const m_volume;
arrayView1d< real64 const > const m_deltaVolume;
Expand Down
Loading
Loading