From 0a5f7985116e325c1078017386537599eacf6c4f Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sat, 14 Mar 2026 10:34:04 -0700 Subject: [PATCH 01/15] switch from molarity concentration to molality --- src/coreComponents/constitutive/HPCReact | 2 +- .../ReactiveSinglePhaseFluid.cpp | 10 +++++ .../ReactiveSinglePhaseFluid.hpp | 4 ++ .../SinglePhaseReactiveTransport.cpp | 39 +++++++++++++++---- .../SinglePhaseReactiveTransportFields.hpp | 6 +-- .../reactive/AccumulationKernels.hpp | 22 +++++------ .../reactive/DirichletFluxComputeKernel.hpp | 22 +++++++---- .../reactive/FluxComputeKernel.hpp | 35 +++++++++++------ .../reactive/SourceFluxComputeKernel.hpp | 10 +++-- .../reactive/ThermalAccumulationKernels.hpp | 5 ++- .../ThermalDirichletFluxComputeKernel.hpp | 4 ++ .../reactive/ThermalFluxComputeKernel.hpp | 17 +++++--- .../ThermalSourceFluxComputeKernel.hpp | 3 +- 13 files changed, 124 insertions(+), 55 deletions(-) diff --git a/src/coreComponents/constitutive/HPCReact b/src/coreComponents/constitutive/HPCReact index 5268dc5c2a5..d844e5851f8 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit 5268dc5c2a59e9dae23f435fba9aecf5a5c5de33 +Subproject commit d844e5851f83e6724dd4f8719e2e8fc20c7b7be0 diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index 69b7b750934..62120595ca1 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -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 >:: @@ -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; } @@ -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; } } diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index 6346d030c4f..e246f0a0155 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -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. */ @@ -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 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 5d543342aaa..3bb9edda8db 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -463,7 +463,7 @@ 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(); @@ -471,6 +471,25 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, 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(); @@ -481,6 +500,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies, m_hasDiffusion, mobilePrimarySpeciesFlags.toViewConst(), + solventDensity, dofManager.rankOffset(), dofKey, getName(), @@ -496,6 +516,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies, m_hasDiffusion, mobilePrimarySpeciesFlags.toViewConst(), + solventDensity, dofManager.rankOffset(), dofKey, getName(), @@ -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; } } ); } @@ -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; } } ); } @@ -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; } } ); } @@ -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; } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp index 28c72f5c170..0d6a47880a0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp @@ -41,7 +41,7 @@ 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", @@ -49,7 +49,7 @@ DECLARE_FIELD( logPrimarySpeciesConcentration_n, 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", @@ -57,7 +57,7 @@ DECLARE_FIELD( bcLogPrimarySpeciesConcentration, 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", diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index 28c3b1cc181..5df0803302f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -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() ), @@ -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] ) * @@ -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; } } } @@ -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; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp index 80ff9bce9ba..fac64cc70da 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp @@ -117,6 +117,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo ReactiveSinglePhaseFluidAccessors const & reactiveSinglePhaseFluidAccessors, PermeabilityAccessors const & permeabilityAccessors, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -135,7 +136,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo m_primarySpeciesMobileAggregateConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::primarySpeciesMobileAggregateConcentration {} ) ), m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations {} ) ), - m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ) + m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ), + m_solventDensity( solventDensity ) {} /** @@ -194,17 +196,17 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is]; - speciesFlux[is] = aggregateConc_i / dens_up * fluxVal * mobility_up; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is] * m_solventDensity; + speciesFlux[is] = aggregateConcMolarity_i / dens_up * fluxVal * mobility_up; - dSpeciesFlux_dP[is] = aggregateConc_i / dens_up * dFlux_dP * mobility_up - + aggregateConc_i / dens_up * fluxVal * dMobility_dP_up - - aggregateConc_i * fluxVal * mobility_up * dDens_dP_up / (dens_up * dens_up); + dSpeciesFlux_dP[is] = aggregateConcMolarity_i / dens_up * dFlux_dP * mobility_up + + aggregateConcMolarity_i / dens_up * fluxVal * dMobility_dP_up + - aggregateConcMolarity_i * fluxVal * mobility_up * dDens_dP_up / (dens_up * dens_up); for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[seri][sesri][sei][0][is][js]; - dSpeciesFlux_dLogConc[is][js] += dAggregateConc_i_dLogConc_j / dens_up * fluxVal * mobility_up; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[seri][sesri][sei][0][is][js] * m_solventDensity; + dSpeciesFlux_dLogConc[is][js] += dAggregateConcMolarity_i_dLogConc_j / dens_up * fluxVal * mobility_up; } } @@ -267,6 +269,9 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; + /// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³] + real64 const m_solventDensity; + }; @@ -344,6 +349,7 @@ class DirichletFluxComputeKernelFactory reactiveFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, + reactiveFluid.solventDensity(), dt, localMatrix, localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp index a408c29aa0d..2fdcc47d656 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp @@ -122,6 +122,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E * @param[in] porosityAccessors * @param[in] hasDiffusion the flag to turn on diffusion calculation * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species + * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector @@ -138,6 +139,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E PorosityAccessors const & porosityAccessors, integer const & hasDiffusion, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -159,7 +161,8 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E m_dDiffusivity_dTemp( diffusionAccessors.get( fields::diffusion::dDiffusivity_dTemperature {} ) ), m_referencePorosity( porosityAccessors.get( fields::porosity::referencePorosity {} ) ), m_hasDiffusion( hasDiffusion ), - m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ) + m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ), + m_solventDensity( solventDensity ) {} /** @@ -249,20 +252,22 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E // compute species fluxes and derivatives using upstream cell concentration for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is]; - speciesFlux[is] = aggregateConc_i / fluidDens_up * fluxVal; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] + * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density + speciesFlux[is] = aggregateConcMolarity_i / fluidDens_up * fluxVal; for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { - dSpeciesFlux_dP[ke][is] += aggregateConc_i / fluidDens_up * dFlux_dP[ke]; + dSpeciesFlux_dP[ke][is] += aggregateConcMolarity_i / fluidDens_up * dFlux_dP[ke]; } - dSpeciesFlux_dP[k_up][is] += -aggregateConc_i * fluxVal * dDens_dPres / (fluidDens_up * fluidDens_up); + dSpeciesFlux_dP[k_up][is] += -aggregateConcMolarity_i * fluxVal * dDens_dPres / (fluidDens_up * fluidDens_up); for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js]; - dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConc_i_dLogConc_j / fluidDens_up * fluxVal; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js] + * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density + dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConcMolarity_i_dLogConc_j / fluidDens_up * fluxVal; } } @@ -354,15 +359,16 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E localIndex const esr = sesri[ke]; localIndex const ei = sei[ke]; - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is]; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] + * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density - speciesGrad[is] += diffusionTrans[ke] * aggregateConc_i; + speciesGrad[is] += diffusionTrans[ke] * aggregateConcMolarity_i; for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js]; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js] * m_solventDensity; - dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConc_i_dLogConc_j; + dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConcMolarity_i_dLogConc_j; } } @@ -500,6 +506,9 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; + + /// Density of the solvent (e.g., water) [kg/m3] + real64 const m_solventDensity; }; /** @@ -516,6 +525,7 @@ class FluxComputeKernelFactory * @param[in] numSpecies the number of primary species * @param[in] hasDiffusion the flag of adding diffusion term * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species + * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] * @param[in] rankOffset the offset of my MPI rank * @param[in] dofKey string to get the element degrees of freedom numbers * @param[in] solverName name of the solver (to name accessors) @@ -530,6 +540,7 @@ class FluxComputeKernelFactory createAndLaunch( integer const numSpecies, integer const hasDiffusion, arrayView1d< integer const > const mobilePrimarySpeciesFlags, + real64 const solventDensity, globalIndex const rankOffset, string const & dofKey, string const & solverName, @@ -561,7 +572,7 @@ class FluxComputeKernelFactory KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor, flowAccessors, reactiveFlowAccessors, fluidAccessors, reactiveFluidAccessors, permAccessors, diffusionAccessors, porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags, - dt, localMatrix, localRhs ); + solventDensity, dt, localMatrix, localRhs ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp index 1ef58f942aa..8f2897a451f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp @@ -74,6 +74,7 @@ class SourceFluxComputeKernel m_elemGhostRank( elemGhostRank ), m_rhsContributionArrayView( rhsContributionArrayView ), m_sizeScalingFactor( sizeScalingFactor ), + m_solventDensity( fluid.solventDensity() ), m_primarySpeciesAggregateConcentration( fluid.primarySpeciesAggregateConcentration() ), m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations( fluid.dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations() ), m_density( fluid.density() ), @@ -151,12 +152,12 @@ class SourceFluxComputeKernel for( integer i = 0; i < numSpecies; ++i ) { - stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] / m_density[ei][0] * scaledInflowMass; - stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; + stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity / m_density[ei][0] * scaledInflowMass; + stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; for( integer j = 0; j < numSpecies; ++j ) { - stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] / m_density[ei][0] * scaledInflowMass; + stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventDensity / m_density[ei][0] * scaledInflowMass; } } } @@ -237,6 +238,9 @@ class SourceFluxComputeKernel /// size scaling factor real64 const m_sizeScalingFactor; + /// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³] + real64 const m_solventDensity; + // View on the total concentration of ions that contain the primary species arrayView3d< real64 const, constitutive::reactivefluid::USD_SPECIES > const m_primarySpeciesAggregateConcentration; // View on the derivatives of total ion concentration for the primary species wrt log of primary species concentration diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp index 33fa2ac0f3b..276ea4bea5d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp @@ -53,6 +53,7 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe using Base::m_volume; using Base::m_deltaVolume; using Base::m_primarySpeciesAggregateConcentration; + using Base::m_solventDensity; /// Note: Derivative lineup only supports dP & dT, not component terms using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; @@ -161,9 +162,9 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe for( integer is = 0; is < numSpecies; ++is ) { // Drivative of primary species amount in pore volume wrt temperature - stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][0][is] + stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity /* + stack.poreVolume * - m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] */; + m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] * m_solventDensity */; // // Derivative of reaction term wrt temperature // stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * // m_dPrimarySpeciesTotalKineticRate_dTemp[is]; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp index 7a87bdfbc50..73d7135862c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp @@ -86,6 +86,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl using Base::m_sei; using Base::m_facePres; using Base::m_faceGravCoef; + using Base::m_solventDensity; using ReactiveSinglePhaseFlowAccessors = typename Base::ReactiveSinglePhaseFlowAccessors; using ReactiveSinglePhaseFluidAccessors = typename Base::ReactiveSinglePhaseFluidAccessors; @@ -137,6 +138,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl PermeabilityAccessors const & permeabilityAccessors, ThermalConductivityAccessors const & thermalConductivityAccessors, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -152,6 +154,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl reactiveSinglePhaseFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, + solventDensity, dt, localMatrix, localRhs ), @@ -392,6 +395,7 @@ class DirichletFluxComputeKernelFactory permeabilityAccessors, thermalConductivityAccessors, mobilePrimarySpeciesFlags, + reactiveFluid.solventDensity(), dt, localMatrix, localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp index d436b261628..80c262efcd4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp @@ -89,6 +89,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne using Base::m_primarySpeciesMobileAggregateConc; using Base::m_referencePorosity; using Base::m_mobilePrimarySpeciesFlags; + using Base::m_solventDensity; using ThermalSinglePhaseFlowAccessors = StencilAccessors< fields::flow::temperature >; @@ -121,6 +122,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne * @param[in] thermalConductivityAccessors accessor for wrappers registered by the thermal conductivity model * @param[in] hasDiffusion the flag to turn on diffusion calculation * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species + * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector @@ -140,6 +142,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne ThermalConductivityAccessors const & thermalConductivityAccessors, integer const & hasDiffusion, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -155,6 +158,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags, + solventDensity, dt, localMatrix, localRhs ), @@ -328,15 +332,15 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne // Step 2.2: compute speciesFlux derivative wrt temperature for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is]; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] * m_solventDensity; // real64 const dAggregateConc_i_dTemp = m_dPrimarySpeciesMobileAggregateConcentration_dTemp[er_up][esr_up][ei_up][is]; - // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * fluxVal / fluidDens_up; - dSpeciesFlux_dT[k_up][is] += -aggregateConc_i * fluxVal * dDens_dTemp / (fluidDens_up * fluidDens_up); + // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * m_solventDensity * fluxVal / fluidDens_up; + dSpeciesFlux_dT[k_up][is] += -aggregateConcMolarity_i * fluxVal * dDens_dTemp / (fluidDens_up * fluidDens_up); for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { - dSpeciesFlux_dT[ke][is] += aggregateConc_i / fluidDens_up * dFlux_dT[ke]; + dSpeciesFlux_dT[ke][is] += aggregateConcMolarity_i / fluidDens_up * dFlux_dT[ke]; } } } @@ -489,7 +493,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne // dSpeciesGrad_dT[ke] += stack.diffusionTransmissibility[connectionIndex][ke] // * m_dPrimarySpeciesMobileAggregateConcentration_dTemp[er][esr][ei][is]; - dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is]; + dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventDensity; } for( integer ke = 0; ke < numFluxSupportPoints; ke++ ) @@ -588,6 +592,7 @@ class FluxComputeKernelFactory createAndLaunch( integer const numSpecies, integer const hasDiffusion, arrayView1d< integer const > const mobilePrimarySpeciesFlags, + real64 const solventDensity, globalIndex const rankOffset, string const & dofKey, string const & solverName, @@ -622,7 +627,7 @@ class FluxComputeKernelFactory KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor, flowAccessors, reactiveFlowAccessors, thermalFlowAccessors, fluidAccessors, reactiveFluidAccessors, thermalFluidAccessors, permAccessors, diffusionAccessors, porosityAccessors, thermalConductivityAccessors, - hasDiffusion, mobilePrimarySpeciesFlags, dt, localMatrix, localRhs ); + hasDiffusion, mobilePrimarySpeciesFlags, solventDensity, dt, localMatrix, localRhs ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp index 5fa7f2c201c..21f542c262a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp @@ -45,6 +45,7 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu using Base::numDof; using Base::numEqn; using Base::m_sizeScalingFactor; + using Base::m_solventDensity; using Base::m_primarySpeciesAggregateConcentration; using Base::m_density; using Base::m_dDensity; @@ -122,7 +123,7 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu for( integer i = 0; i < numSpecies; ++i ) { - stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_dDensity[ei][0][DerivOffset::dT] / (m_density[ei][0] * m_density[ei][0]) * + stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dT] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; } } From befd8a996638e3576c506334708c71e3d585376f Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Mon, 16 Mar 2026 13:10:33 -0700 Subject: [PATCH 02/15] make the timestep in unit test smaller as the jacobian is very sensitive to that --- .../fluidFlowTests/testSinglePhaseReactiveTransport.cpp | 2 +- .../kernels/singlePhase/reactive/AccumulationKernels.hpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp index c4b50dfab30..20e3980fe89 100644 --- a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp @@ -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; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index 5df0803302f..e3e3dee58c8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -161,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 ); From 2e211f4ee7f02b37e8351d0341b91b35de3418f0 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Tue, 1 Sep 2026 09:40:58 -0700 Subject: [PATCH 03/15] single reaction system across all regions to ensure consistent dof layout --- .../SinglePhaseReactiveTransport.cpp | 40 +++++++++---------- 1 file changed, 18 insertions(+), 22 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index c2359370af9..9ba3bd003b8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -127,12 +127,22 @@ void SinglePhaseReactiveTransport::registerDataOnMesh( Group & meshBodies ) [&]( localIndex const, ElementSubRegionBase & subRegion ) { + string const reactiveFluidModelName = m_isThermal? getConstitutiveName< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion ): + getConstitutiveName< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion ); + if( m_reactiveFluidModelName.empty() ) { - m_reactiveFluidModelName = m_isThermal? getConstitutiveName< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion ): - getConstitutiveName< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion ); + m_reactiveFluidModelName = reactiveFluidModelName; } + // The number of species, the dof layout and the solvent density are all taken from a single reactive + // fluid model, and the flux kernels are launched per stencil rather than per region, so every subregion + // must share that model. + GEOS_THROW_IF_NE_MSG( reactiveFluidModelName, m_reactiveFluidModelName, + GEOS_FMT( "SinglePhaseReactiveTransport {}: all regions must use the same reactive fluid model, but {} uses a different one", + getDataContext(), subRegion.getDataContext() ), + InputError ); + // If at least one region has a diffusion model, consider it enabled for all string const diffusionName = getConstitutiveName< DiffusionBase >( subRegion ); if( !diffusionName.empty() ) @@ -461,9 +471,14 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, } } + ConstitutiveManager const & cm = domain.getConstitutiveManager(); + real64 const solventDensity = + m_isThermal ? cm.getConstitutiveRelation< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventDensity() + : cm.getConstitutiveRelation< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventDensity(); + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, - string_array const & regionNames ) + string_array const & ) { NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); @@ -471,25 +486,6 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, 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(); From c366f3556d5f3091c071ceca2308298a3321ea84 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Tue, 1 Sep 2026 09:43:25 -0700 Subject: [PATCH 04/15] fix kinetic reaction rate with molality, uncrustify --- .../fluidFlowTests/testSinglePhaseReactiveTransport.cpp | 2 +- .../fluidFlow/SinglePhaseReactiveTransport.cpp | 6 ++---- .../kernels/singlePhase/reactive/AccumulationKernels.hpp | 6 +++--- .../kernels/singlePhase/reactive/FluxComputeKernel.hpp | 2 +- .../singlePhase/reactive/SourceFluxComputeKernel.hpp | 6 ++++-- .../singlePhase/reactive/ThermalAccumulationKernels.hpp | 3 ++- .../singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp | 3 ++- 7 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp index 20e3980fe89..c4b50dfab30 100644 --- a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp @@ -530,7 +530,7 @@ class SinglePhaseReactiveTransportTest : public ::testing::Test } static real64 constexpr time = 0.0; - static real64 constexpr dt = 0.01; + static real64 constexpr dt = 1.0; static real64 constexpr eps = std::numeric_limits< real64 >::epsilon(); GeosxState state; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 9ba3bd003b8..4ffb0f06cd8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -617,13 +617,12 @@ 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] * solventDensity; + kineticReactionMolarIncrements[ei][r] = dt * kineticReactionRates[ei][0][r]; } } ); } @@ -632,13 +631,12 @@ 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] * solventDensity; + kineticReactionMolarIncrements[ei][r] = dt * kineticReactionRates[ei][0][r]; } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index e3e3dee58c8..26f607f4bf5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -162,7 +162,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU StackVariables & stack ) const { // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventDensity * stack.poreVolume - primarySpeciesAggregateMole_n[is]) - // - dt * m_volume * primarySpeciesKineticRate[is] * solventDensity + // - dt * m_volume * primarySpeciesKineticRate[is] Base::computeAccumulation( ei, stack ); @@ -187,7 +187,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU 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] * m_solventDensity; + stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * m_primarySpeciesAggregateKineticRate[ei][0][is]; // Step 2.1: jacobian // Drivative of primary species amount in pore volume wrt pressure @@ -204,7 +204,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU m_primarySpeciesAggregateConcentration[ei][0][is] + */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; + stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js]; } } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp index 2fdcc47d656..8beaff2ab23 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp @@ -265,7 +265,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js] + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js] * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConcMolarity_i_dLogConc_j / fluidDens_up * fluxVal; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp index 8f2897a451f..76212f9dc01 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp @@ -153,11 +153,13 @@ class SourceFluxComputeKernel for( integer i = 0; i < numSpecies; ++i ) { stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity / m_density[ei][0] * scaledInflowMass; - stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; + stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * + scaledInflowMass; for( integer j = 0; j < numSpecies; ++j ) { - stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventDensity / m_density[ei][0] * scaledInflowMass; + stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventDensity / m_density[ei][0] * + scaledInflowMass; } } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp index 276ea4bea5d..39cc407199a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp @@ -164,7 +164,8 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe // Drivative of primary species amount in pore volume wrt temperature stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity /* + stack.poreVolume * - m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] * m_solventDensity */; + m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] * + m_solventDensity */; // // Derivative of reaction term wrt temperature // stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * // m_dPrimarySpeciesTotalKineticRate_dTemp[is]; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp index 21f542c262a..b5280db1eb4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp @@ -123,7 +123,8 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu for( integer i = 0; i < numSpecies; ++i ) { - stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dT] / (m_density[ei][0] * m_density[ei][0]) * + stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dT] / + (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; } } From 816c04822afe9f6387e9602936d6f0640934339d Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Wed, 2 Sep 2026 13:33:09 -0700 Subject: [PATCH 05/15] rename to solventMassPerSolutionVolume to be more precise, and define it in GEOS rather than HPCReact for now --- .../ChainSerial_1DAllKinetic_base.xml | 3 ++- .../MoMaSEasy_1DAdvective_base.xml | 3 ++- .../MoMaSEasy_1DDiffusive_base.xml | 3 ++- .../MoMaSEasy_2DAdvective_base.xml | 3 ++- .../MoMaSMedium_1DAdvective_base.xml | 3 ++- src/coreComponents/constitutive/HPCReact | 2 +- .../ReactiveSinglePhaseFluid.cpp | 20 ++++++++++------ .../ReactiveSinglePhaseFluid.hpp | 21 ++++++++++++++-- .../SinglePhaseReactiveTransport.cpp | 24 +++++++++---------- .../reactive/AccumulationKernels.hpp | 14 +++++------ .../reactive/DirichletFluxComputeKernel.hpp | 14 +++++------ .../reactive/FluxComputeKernel.hpp | 24 +++++++++---------- .../reactive/SourceFluxComputeKernel.hpp | 12 +++++----- .../reactive/ThermalAccumulationKernels.hpp | 6 ++--- .../ThermalDirichletFluxComputeKernel.hpp | 8 +++---- .../reactive/ThermalFluxComputeKernel.hpp | 18 +++++++------- .../ThermalSourceFluxComputeKernel.hpp | 4 ++-- src/coreComponents/schema/schema.xsd | 4 ++++ 18 files changed, 109 insertions(+), 77 deletions(-) diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml index 27df60a6124..c58004bb4f4 100644 --- a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml @@ -44,7 +44,8 @@ referencePressure="0.0" compressibility="0.0" viscosibility="0.0" - chemicalSystemType="chainSerialAllKinetic"/> + chemicalSystemType="chainSerialAllKinetic" + solventMassPerSolutionVolume="1"/> + chemicalSystemType="momasEasy" + solventMassPerSolutionVolume="1"/> + chemicalSystemType="momasEasy" + solventMassPerSolutionVolume="1"/> + chemicalSystemType="momasEasy" + solventMassPerSolutionVolume="1"/> + chemicalSystemType="momasMedium" + solventMassPerSolutionVolume="1"/> ::concat( "|" ) + "``" ); + this->registerWrapper( viewKeyStruct::solventMassPerSolutionVolumeString(), &m_solventMassPerSolutionVolume ). + setApplyDefaultValue( 1000.0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Mass of solvent per unit volume of solution [kg/m^3], used to convert species " + "molality [mol/kg solvent] to molarity [mol/m^3 solution]. The default of 1000 " + "approximates an aqueous solution by the density of pure water." ); + this->template registerField< fields::reactivefluid::initialPrimarySpeciesConcentration >( &m_initialPrimarySpeciesConcentration ); this->template registerField< fields::reactivefluid::secondarySpeciesConcentration >( &m_secondarySpeciesConcentration ); this->template registerField< fields::reactivefluid::primarySpeciesAggregateConcentration >( &m_primarySpeciesAggregateConcentration ); @@ -72,7 +79,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; + newConstitutiveRelation.m_solventMassPerSolutionVolume = m_solventMassPerSolutionVolume; return clone; } @@ -88,44 +95,43 @@ 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; } + + GEOS_THROW_IF_LE_MSG( m_solventMassPerSolutionVolume, 0.0, + GEOS_FMT( "invalid value of attribute '{}'", + viewKeyStruct::solventMassPerSolutionVolumeString() ), + InputError, this->getDataContext() ); } template< typename BASE > diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index e246f0a0155..50ed70db6f1 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -113,7 +113,14 @@ class ReactiveSinglePhaseFluid : public BASE integer numKineticReactions() const { return m_numKineticReactions; } - real64 solventDensity() const { return m_solventDensity; } + /** + * @brief Mass of solvent per unit volume of solution [kg/m^3]. + * + * Converts species molality [mol/kg solvent] to molarity [mol/m^3 solution]. HPCReact is a + * molality-based library: concentrations, equilibrium constants and mass-action quotients are all + * on the molal scale. + */ + real64 solventMassPerSolutionVolume() const { return m_solventMassPerSolutionVolume; } /** * @brief Kernel wrapper class for ReactiveSinglePhaseFluid. @@ -332,6 +339,7 @@ class ReactiveSinglePhaseFluid : public BASE struct viewKeyStruct : ConstitutiveBase::viewKeyStruct { static constexpr char const * chemicalSystemNameString() { return "chemicalSystemType"; } + static constexpr char const * solventMassPerSolutionVolumeString() { return "solventMassPerSolutionVolume"; } }; protected: @@ -368,7 +376,16 @@ class ReactiveSinglePhaseFluid : public BASE ChemicalSystemType m_chemicalSystemType; - real64 m_solventDensity; + /// TODO: prescribed as a constant for now. The exact factor is + /// + /// rho_s = rho * w + /// + /// where rho_s is this quantity [kg/m^3], rho the solution density [kg/m^3] and w the + /// solvent mass fraction [-]. For the carbonate brine EQ3/6 gives 1070.9 * 0.898 = 961.6, + /// not the 1000 defaulted here. Ideally rho is a function of pressure, temperature and + /// species concentration, and w a function of concentration. The update methods and where + /// they should be launched are TBD. + real64 m_solventMassPerSolutionVolume; }; // these aliases are useful in constitutive dispatch diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 4ffb0f06cd8..97732290591 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -135,7 +135,7 @@ void SinglePhaseReactiveTransport::registerDataOnMesh( Group & meshBodies ) m_reactiveFluidModelName = reactiveFluidModelName; } - // The number of species, the dof layout and the solvent density are all taken from a single reactive + // The number of species, the dof layout and the molality-to-molarity conversion factor are all taken from a single reactive // fluid model, and the flux kernels are launched per stencil rather than per region, so every subregion // must share that model. GEOS_THROW_IF_NE_MSG( reactiveFluidModelName, m_reactiveFluidModelName, @@ -472,9 +472,9 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, } ConstitutiveManager const & cm = domain.getConstitutiveManager(); - real64 const solventDensity = - m_isThermal ? cm.getConstitutiveRelation< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventDensity() - : cm.getConstitutiveRelation< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventDensity(); + real64 const solventMassPerSolutionVolume = + m_isThermal ? cm.getConstitutiveRelation< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventMassPerSolutionVolume() + : cm.getConstitutiveRelation< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventMassPerSolutionVolume(); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, @@ -496,7 +496,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies, m_hasDiffusion, mobilePrimarySpeciesFlags.toViewConst(), - solventDensity, + solventMassPerSolutionVolume, dofManager.rankOffset(), dofKey, getName(), @@ -512,7 +512,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies, m_hasDiffusion, mobilePrimarySpeciesFlags.toViewConst(), - solventDensity, + solventMassPerSolutionVolume, dofManager.rankOffset(), dofKey, getName(), @@ -569,16 +569,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(); + real64 const solventMassPerSolutionVolume = fluid.solventMassPerSolutionVolume(); 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] * solventDensity; + primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventMassPerSolutionVolume; if( isZero( primarySpeciesAggregateMole_n[ei][is] ) ) - primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventDensity; + primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventMassPerSolutionVolume; } } ); } @@ -588,16 +588,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(); + real64 const solventMassPerSolutionVolume = fluid.solventMassPerSolutionVolume(); 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] * solventDensity; + primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventMassPerSolutionVolume; if( isZero( primarySpeciesAggregateMole_n[ei][is] ) ) - primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventDensity; + primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventMassPerSolutionVolume; } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index 26f607f4bf5..4effb6c1cdf 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -82,7 +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_solventMassPerSolutionVolume( fluid.solventMassPerSolutionVolume() ), m_volume( subRegion.getElementVolume() ), m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ), m_porosity( solid.getPorosity() ), @@ -161,7 +161,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU void computeAccumulation( localIndex const ei, StackVariables & stack ) const { - // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventDensity * stack.poreVolume - primarySpeciesAggregateMole_n[is]) + // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventMassPerSolutionVolume * stack.poreVolume - primarySpeciesAggregateMole_n[is]) // - dt * m_volume * primarySpeciesKineticRate[is] Base::computeAccumulation( ei, stack ); @@ -184,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] * m_solventDensity * stack.poreVolume; + stack.localResidual[is+numEqn-numSpecies] += m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassPerSolutionVolume * stack.poreVolume; // Reaction term stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * m_primarySpeciesAggregateKineticRate[ei][0][is]; // 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] * m_solventDensity + stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassPerSolutionVolume /* + 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] ) * @@ -202,7 +202,7 @@ 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] * m_solventDensity; + + */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[is][js] * m_solventMassPerSolutionVolume; stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js]; } @@ -240,8 +240,8 @@ 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; + /// Mass of solvent per unit volume of solution [kg/m³], converting molality [mol/kg] to molarity [mol/m³] + real64 const m_solventMassPerSolutionVolume; /// View on the element volumes arrayView1d< real64 const > const m_volume; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp index fac64cc70da..2ef73fddea2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp @@ -117,7 +117,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo ReactiveSinglePhaseFluidAccessors const & reactiveSinglePhaseFluidAccessors, PermeabilityAccessors const & permeabilityAccessors, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, - real64 const & solventDensity, + real64 const & solventMassPerSolutionVolume, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -137,7 +137,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations {} ) ), m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ), - m_solventDensity( solventDensity ) + m_solventMassPerSolutionVolume( solventMassPerSolutionVolume ) {} /** @@ -196,7 +196,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is] * m_solventDensity; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is] * m_solventMassPerSolutionVolume; speciesFlux[is] = aggregateConcMolarity_i / dens_up * fluxVal * mobility_up; dSpeciesFlux_dP[is] = aggregateConcMolarity_i / dens_up * dFlux_dP * mobility_up @@ -205,7 +205,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[seri][sesri][sei][0][is][js] * m_solventDensity; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[seri][sesri][sei][0][is][js] * m_solventMassPerSolutionVolume; dSpeciesFlux_dLogConc[is][js] += dAggregateConcMolarity_i_dLogConc_j / dens_up * fluxVal * mobility_up; } } @@ -269,8 +269,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; - /// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³] - real64 const m_solventDensity; + /// Mass of solvent per unit volume of solution [kg/m³], converting molality [mol/kg] to molarity [mol/m³] + real64 const m_solventMassPerSolutionVolume; }; @@ -349,7 +349,7 @@ class DirichletFluxComputeKernelFactory reactiveFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, - reactiveFluid.solventDensity(), + reactiveFluid.solventMassPerSolutionVolume(), dt, localMatrix, localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp index 8beaff2ab23..dcc870ff8c5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp @@ -122,7 +122,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E * @param[in] porosityAccessors * @param[in] hasDiffusion the flag to turn on diffusion calculation * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species - * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] + * @param[in] solventMassPerSolutionVolume mass of solvent per unit volume of solution [kg/m3] * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector @@ -139,7 +139,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E PorosityAccessors const & porosityAccessors, integer const & hasDiffusion, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, - real64 const & solventDensity, + real64 const & solventMassPerSolutionVolume, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -162,7 +162,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E m_referencePorosity( porosityAccessors.get( fields::porosity::referencePorosity {} ) ), m_hasDiffusion( hasDiffusion ), m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ), - m_solventDensity( solventDensity ) + m_solventMassPerSolutionVolume( solventMassPerSolutionVolume ) {} /** @@ -253,7 +253,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E for( integer is = 0; is < numSpecies; ++is ) { real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] - * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density + * m_solventMassPerSolutionVolume; // convert from mol/kg to mol/m3 speciesFlux[is] = aggregateConcMolarity_i / fluidDens_up * fluxVal; for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) @@ -266,7 +266,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E for( integer js = 0; js < numSpecies; ++js ) { real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js] - * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density + * m_solventMassPerSolutionVolume; // convert from mol/kg to mol/m3 dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConcMolarity_i_dLogConc_j / fluidDens_up * fluxVal; } } @@ -360,13 +360,13 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E localIndex const ei = sei[ke]; real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] - * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density + * m_solventMassPerSolutionVolume; // convert from mol/kg to mol/m3 speciesGrad[is] += diffusionTrans[ke] * aggregateConcMolarity_i; for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js] * m_solventDensity; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js] * m_solventMassPerSolutionVolume; dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConcMolarity_i_dLogConc_j; } @@ -507,8 +507,8 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; - /// Density of the solvent (e.g., water) [kg/m3] - real64 const m_solventDensity; + /// Mass of solvent per unit volume of solution [kg/m3], converting molality [mol/kg] to molarity [mol/m3] + real64 const m_solventMassPerSolutionVolume; }; /** @@ -525,7 +525,7 @@ class FluxComputeKernelFactory * @param[in] numSpecies the number of primary species * @param[in] hasDiffusion the flag of adding diffusion term * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species - * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] + * @param[in] solventMassPerSolutionVolume mass of solvent per unit volume of solution [kg/m3] * @param[in] rankOffset the offset of my MPI rank * @param[in] dofKey string to get the element degrees of freedom numbers * @param[in] solverName name of the solver (to name accessors) @@ -540,7 +540,7 @@ class FluxComputeKernelFactory createAndLaunch( integer const numSpecies, integer const hasDiffusion, arrayView1d< integer const > const mobilePrimarySpeciesFlags, - real64 const solventDensity, + real64 const solventMassPerSolutionVolume, globalIndex const rankOffset, string const & dofKey, string const & solverName, @@ -572,7 +572,7 @@ class FluxComputeKernelFactory KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor, flowAccessors, reactiveFlowAccessors, fluidAccessors, reactiveFluidAccessors, permAccessors, diffusionAccessors, porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags, - solventDensity, dt, localMatrix, localRhs ); + solventMassPerSolutionVolume, dt, localMatrix, localRhs ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp index 76212f9dc01..d780e5b2b24 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp @@ -74,7 +74,7 @@ class SourceFluxComputeKernel m_elemGhostRank( elemGhostRank ), m_rhsContributionArrayView( rhsContributionArrayView ), m_sizeScalingFactor( sizeScalingFactor ), - m_solventDensity( fluid.solventDensity() ), + m_solventMassPerSolutionVolume( fluid.solventMassPerSolutionVolume() ), m_primarySpeciesAggregateConcentration( fluid.primarySpeciesAggregateConcentration() ), m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations( fluid.dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations() ), m_density( fluid.density() ), @@ -152,13 +152,13 @@ class SourceFluxComputeKernel for( integer i = 0; i < numSpecies; ++i ) { - stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity / m_density[ei][0] * scaledInflowMass; - stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * + stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventMassPerSolutionVolume / m_density[ei][0] * scaledInflowMass; + stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventMassPerSolutionVolume * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; for( integer j = 0; j < numSpecies; ++j ) { - stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventDensity / m_density[ei][0] * + stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventMassPerSolutionVolume / m_density[ei][0] * scaledInflowMass; } } @@ -240,8 +240,8 @@ class SourceFluxComputeKernel /// size scaling factor real64 const m_sizeScalingFactor; - /// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³] - real64 const m_solventDensity; + /// Mass of solvent per unit volume of solution [kg/m³], converting molality [mol/kg] to molarity [mol/m³] + real64 const m_solventMassPerSolutionVolume; // View on the total concentration of ions that contain the primary species arrayView3d< real64 const, constitutive::reactivefluid::USD_SPECIES > const m_primarySpeciesAggregateConcentration; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp index 39cc407199a..b22932096a8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp @@ -53,7 +53,7 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe using Base::m_volume; using Base::m_deltaVolume; using Base::m_primarySpeciesAggregateConcentration; - using Base::m_solventDensity; + using Base::m_solventMassPerSolutionVolume; /// Note: Derivative lineup only supports dP & dT, not component terms using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; @@ -162,10 +162,10 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe for( integer is = 0; is < numSpecies; ++is ) { // Drivative of primary species amount in pore volume wrt temperature - stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity + stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassPerSolutionVolume /* + stack.poreVolume * m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] * - m_solventDensity */; + m_solventMassPerSolutionVolume */; // // Derivative of reaction term wrt temperature // stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * // m_dPrimarySpeciesTotalKineticRate_dTemp[is]; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp index 73d7135862c..740277c402f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp @@ -86,7 +86,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl using Base::m_sei; using Base::m_facePres; using Base::m_faceGravCoef; - using Base::m_solventDensity; + using Base::m_solventMassPerSolutionVolume; using ReactiveSinglePhaseFlowAccessors = typename Base::ReactiveSinglePhaseFlowAccessors; using ReactiveSinglePhaseFluidAccessors = typename Base::ReactiveSinglePhaseFluidAccessors; @@ -138,7 +138,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl PermeabilityAccessors const & permeabilityAccessors, ThermalConductivityAccessors const & thermalConductivityAccessors, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, - real64 const & solventDensity, + real64 const & solventMassPerSolutionVolume, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -154,7 +154,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl reactiveSinglePhaseFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, - solventDensity, + solventMassPerSolutionVolume, dt, localMatrix, localRhs ), @@ -395,7 +395,7 @@ class DirichletFluxComputeKernelFactory permeabilityAccessors, thermalConductivityAccessors, mobilePrimarySpeciesFlags, - reactiveFluid.solventDensity(), + reactiveFluid.solventMassPerSolutionVolume(), dt, localMatrix, localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp index 80c262efcd4..6ed2df03064 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp @@ -89,7 +89,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne using Base::m_primarySpeciesMobileAggregateConc; using Base::m_referencePorosity; using Base::m_mobilePrimarySpeciesFlags; - using Base::m_solventDensity; + using Base::m_solventMassPerSolutionVolume; using ThermalSinglePhaseFlowAccessors = StencilAccessors< fields::flow::temperature >; @@ -122,7 +122,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne * @param[in] thermalConductivityAccessors accessor for wrappers registered by the thermal conductivity model * @param[in] hasDiffusion the flag to turn on diffusion calculation * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species - * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] + * @param[in] solventMassPerSolutionVolume mass of solvent per unit volume of solution [kg/m3] * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector @@ -142,7 +142,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne ThermalConductivityAccessors const & thermalConductivityAccessors, integer const & hasDiffusion, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, - real64 const & solventDensity, + real64 const & solventMassPerSolutionVolume, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -158,7 +158,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags, - solventDensity, + solventMassPerSolutionVolume, dt, localMatrix, localRhs ), @@ -332,10 +332,10 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne // Step 2.2: compute speciesFlux derivative wrt temperature for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] * m_solventDensity; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] * m_solventMassPerSolutionVolume; // real64 const dAggregateConc_i_dTemp = m_dPrimarySpeciesMobileAggregateConcentration_dTemp[er_up][esr_up][ei_up][is]; - // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * m_solventDensity * fluxVal / fluidDens_up; + // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * m_solventMassPerSolutionVolume * fluxVal / fluidDens_up; dSpeciesFlux_dT[k_up][is] += -aggregateConcMolarity_i * fluxVal * dDens_dTemp / (fluidDens_up * fluidDens_up); for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) @@ -493,7 +493,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne // dSpeciesGrad_dT[ke] += stack.diffusionTransmissibility[connectionIndex][ke] // * m_dPrimarySpeciesMobileAggregateConcentration_dTemp[er][esr][ei][is]; - dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventDensity; + dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventMassPerSolutionVolume; } for( integer ke = 0; ke < numFluxSupportPoints; ke++ ) @@ -592,7 +592,7 @@ class FluxComputeKernelFactory createAndLaunch( integer const numSpecies, integer const hasDiffusion, arrayView1d< integer const > const mobilePrimarySpeciesFlags, - real64 const solventDensity, + real64 const solventMassPerSolutionVolume, globalIndex const rankOffset, string const & dofKey, string const & solverName, @@ -627,7 +627,7 @@ class FluxComputeKernelFactory KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor, flowAccessors, reactiveFlowAccessors, thermalFlowAccessors, fluidAccessors, reactiveFluidAccessors, thermalFluidAccessors, permAccessors, diffusionAccessors, porosityAccessors, thermalConductivityAccessors, - hasDiffusion, mobilePrimarySpeciesFlags, solventDensity, dt, localMatrix, localRhs ); + hasDiffusion, mobilePrimarySpeciesFlags, solventMassPerSolutionVolume, dt, localMatrix, localRhs ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp index b5280db1eb4..739b6936085 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp @@ -45,7 +45,7 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu using Base::numDof; using Base::numEqn; using Base::m_sizeScalingFactor; - using Base::m_solventDensity; + using Base::m_solventMassPerSolutionVolume; using Base::m_primarySpeciesAggregateConcentration; using Base::m_density; using Base::m_dDensity; @@ -123,7 +123,7 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu for( integer i = 0; i < numSpecies; ++i ) { - stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dT] / + stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventMassPerSolutionVolume * m_dDensity[ei][0][DerivOffset::dT] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; } diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index c088a8aeb0e..7eede4c8c9d 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -8759,6 +8759,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + + + From 078d45f70f7941c7eacc4e9fd3b13797c777d0d2 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Wed, 2 Sep 2026 14:03:26 -0700 Subject: [PATCH 06/15] uncrustify --- .../kernels/singlePhase/reactive/AccumulationKernels.hpp | 6 ++++-- .../singlePhase/reactive/SourceFluxComputeKernel.hpp | 6 ++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index 4effb6c1cdf..983673c90d3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -161,7 +161,8 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU void computeAccumulation( localIndex const ei, StackVariables & stack ) const { - // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventMassPerSolutionVolume * stack.poreVolume - primarySpeciesAggregateMole_n[is]) + // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventMassPerSolutionVolume * stack.poreVolume - + // primarySpeciesAggregateMole_n[is]) // - dt * m_volume * primarySpeciesKineticRate[is] Base::computeAccumulation( ei, stack ); @@ -202,7 +203,8 @@ 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] * m_solventMassPerSolutionVolume; + + */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[is][js] * + m_solventMassPerSolutionVolume; stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js]; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp index d780e5b2b24..a113f1771a9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp @@ -153,12 +153,14 @@ class SourceFluxComputeKernel for( integer i = 0; i < numSpecies; ++i ) { stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventMassPerSolutionVolume / m_density[ei][0] * scaledInflowMass; - stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventMassPerSolutionVolume * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * + stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventMassPerSolutionVolume * m_dDensity[ei][0][DerivOffset::dP] / + (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; for( integer j = 0; j < numSpecies; ++j ) { - stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventMassPerSolutionVolume / m_density[ei][0] * + stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventMassPerSolutionVolume / + m_density[ei][0] * scaledInflowMass; } } From c6228bae6fe72ac8d722c52ea90108e4c27bda71 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Thu, 3 Sep 2026 21:01:44 -0700 Subject: [PATCH 07/15] plugged in activity model --- .../constitutive/CMakeLists.txt | 1 + src/coreComponents/constitutive/HPCReact | 2 +- .../ReactiveFluidSystemSelector.hpp | 267 ++++++++++++++++++ .../ReactiveSinglePhaseFluid.cpp | 64 ++--- .../ReactiveSinglePhaseFluid.hpp | 263 ++++++++--------- .../SinglePhaseReactiveTransport.cpp | 10 +- .../reactive/ReactionUpdateKernel.hpp | 25 +- src/coreComponents/schema/schema.xsd | 9 + 8 files changed, 446 insertions(+), 195 deletions(-) create mode 100644 src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 39ddb39771b..90cd5405359 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -363,6 +363,7 @@ if (ENABLE_HPCREACT) set( constitutive_headers ${constitutive_headers} fluid/reactivefluid/ReactiveFluidSelector.hpp + fluid/reactivefluid/ReactiveFluidSystemSelector.hpp ) set( constitutive_headers diff --git a/src/coreComponents/constitutive/HPCReact b/src/coreComponents/constitutive/HPCReact index 7423e3e2cd5..02afd5d435e 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit 7423e3e2cd55e959d9a823ab8241351f8b1087df +Subproject commit 02afd5d435edcabcd7ac226b837b61038e5da0d3 diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp new file mode 100644 index 00000000000..a70e31407ed --- /dev/null +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp @@ -0,0 +1,267 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ReactiveFluidSystemSelector.hpp + * @brief The chemical systems ReactiveSinglePhaseFluid can instantiate, and the activity model each + * may be paired with. + * + * Adding a system, or an activity model for one, is one entry below plus its name in + * ReactionSystemList. The kernel wrapper variant, the dispatch, and the check that rejects an + * unsupported pairing all follow from that list. + */ +#ifndef GEOS_CONSTITUTIVE_FLUID_REACTIVEFLUID_REACTIVEFLUIDSYSTEMSELECTOR_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_REACTIVEFLUID_REACTIVEFLUIDSYSTEMSELECTOR_HPP_ + +#include "common/DataTypes.hpp" +#include "common/format/EnumStrings.hpp" + +#include "constitutive/HPCReact/src/reactions/geochemistry/GeochemicalSystems.hpp" +#include "constitutive/HPCReact/src/reactions/exampleSystems/BulkGeneric.hpp" +#include "constitutive/HPCReact/src/reactions/exampleSystems/ChainGeneric.hpp" +#include "constitutive/HPCReact/src/reactions/exampleSystems/MoMasBenchmark.hpp" + +#include +#include + +namespace geos +{ + +namespace constitutive +{ + +namespace reactivefluid +{ + +/** + * @brief The chemical system a ReactiveSinglePhaseFluid solves. + */ +enum class ChemicalSystemType : integer +{ + carbonate, + carbonateAllEquilibrium, + ultramafic, + momasEasy, + momasMedium, + chainSerialAllKinetic +}; + +/** + * @brief The activity model applied to it. + * @details ``identity`` is the identity activity model, gamma = 1 and a_w = 1. ``bdot`` needs ion + * size and b-dot parameters, which only the geochemical systems carry. + */ +enum class ActivityModelType : integer +{ + identity, + bdot +}; + +ENUM_STRINGS( ChemicalSystemType, + "carbonate", + "carbonateAllEquilibrium", + "ultramafic", + "momasEasy", + "momasMedium", + "chainSerialAllKinetic" ); + +ENUM_STRINGS( ActivityModelType, + "identity", + "bdot" ); + +/** + * @brief What every entry below shares: the two HPCReact types, and the two enum values the XML + * input selects it by. + * @tparam SYSTEM the ChemicalSystemType it answers to + * @tparam MODEL the ActivityModelType it answers to + * @tparam REACTION_PARAMS the HPCReact reaction parameters type + * @tparam ACTIVITY_MODEL the HPCReact activity model type + * + * Each entry adds the two parameter objects, as accessors returning by value. They cannot be + * template arguments: the HPCReact objects are namespace-scope constexpr and so have internal + * linkage, which would make the entry a different type in every translation unit. Returning a + * copy also matches how the kernel wrapper stores them. + */ +template< ChemicalSystemType SYSTEM, + ActivityModelType MODEL, + typename REACTION_PARAMS, + typename ACTIVITY_MODEL > +struct ReactionSystemEntry +{ + /// the HPCReact reaction parameters type + using ReactionParamsType = REACTION_PARAMS; + + /// the HPCReact activity model type + using ActivityType = ACTIVITY_MODEL; + + /// the chemical system this answers to + static constexpr ChemicalSystemType chemicalSystem = SYSTEM; + + /// the activity model this answers to + static constexpr ActivityModelType activityModel = MODEL; +}; + +/// Carbonate, 16 species with calcite a pure phase, identity activity model. +struct CarbonateIdentityEntry : ReactionSystemEntry< ChemicalSystemType::carbonate, + ActivityModelType::identity, + hpcReact::geochemistry::carbonateSystemType, + hpcReact::geochemistry::carbonateNosolidIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateNosolidIdentityActivityParams; } +}; + +/// Carbonate, 16 species with calcite a pure phase, B-dot activity model on EQ3/6 parameters. +struct CarbonateBdotEntry : ReactionSystemEntry< ChemicalSystemType::carbonate, + ActivityModelType::bdot, + hpcReact::geochemistry::carbonateSystemType, + hpcReact::geochemistry::carbonateNosolidActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateNosolidActivityParamsEQ36; } +}; + +/// Carbonate, 17 species with calcite among them, identity activity model. +struct CarbonateAllEquilibriumIdentityEntry : ReactionSystemEntry< ChemicalSystemType::carbonateAllEquilibrium, + ActivityModelType::identity, + hpcReact::geochemistry::carbonateSystemAllEquilibriumType, + hpcReact::geochemistry::carbonateIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystemAllEquilibrium; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateIdentityActivityParams; } +}; + +/// Carbonate, 17 species with calcite among them, B-dot activity model on EQ3/6 parameters. +struct CarbonateAllEquilibriumBdotEntry : ReactionSystemEntry< ChemicalSystemType::carbonateAllEquilibrium, + ActivityModelType::bdot, + hpcReact::geochemistry::carbonateSystemAllEquilibriumType, + hpcReact::geochemistry::carbonateActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystemAllEquilibrium; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateActivityParamsEQ36; } +}; + +/// Ultramafic, 20 species with the five minerals pure phases, identity activity model. +struct UltramaficIdentityEntry : ReactionSystemEntry< ChemicalSystemType::ultramafic, + ActivityModelType::identity, + hpcReact::geochemistry::ultramaficSystemType, + hpcReact::geochemistry::ultramaficIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::ultramaficSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::ultramaficIdentityActivityParams; } +}; + +/// Ultramafic, 20 species with the five minerals pure phases, B-dot activity model on EQ3/6 parameters. +struct UltramaficBdotEntry : ReactionSystemEntry< ChemicalSystemType::ultramafic, + ActivityModelType::bdot, + hpcReact::geochemistry::ultramaficSystemType, + hpcReact::geochemistry::ultramaficActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::ultramaficSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::ultramaficActivityParamsEQ36; } +}; + +/// MoMaS easy benchmark, 12 abstract species, identity activity model. +struct MomasEasyIdentityEntry : ReactionSystemEntry< ChemicalSystemType::momasEasy, + ActivityModelType::identity, + hpcReact::MoMasBenchmark::easyCaseType, + hpcReact::MoMasBenchmark::easyCaseIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::MoMasBenchmark::easyCaseParams; } + static ActivityType::Params activityParams() { return hpcReact::MoMasBenchmark::easyCaseIdentityActivityParams; } +}; + +/// MoMaS medium benchmark, 14 abstract species, identity activity model. +struct MomasMediumIdentityEntry : ReactionSystemEntry< ChemicalSystemType::momasMedium, + ActivityModelType::identity, + hpcReact::MoMasBenchmark::mediumCaseType, + hpcReact::MoMasBenchmark::mediumCaseIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::MoMasBenchmark::mediumCaseParams; } + static ActivityType::Params activityParams() { return hpcReact::MoMasBenchmark::mediumCaseIdentityActivityParams; } +}; + +/// Generic serial chain, 3 species all kinetic, identity activity model. +struct ChainSerialIdentityEntry : ReactionSystemEntry< ChemicalSystemType::chainSerialAllKinetic, + ActivityModelType::identity, + hpcReact::ChainGeneric::serialAllKineticType, + hpcReact::ChainGeneric::serialAllKineticIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::ChainGeneric::serialAllKineticParams; } + static ActivityType::Params activityParams() { return hpcReact::ChainGeneric::serialAllKineticIdentityActivityParams; } +}; + +/** + * @brief Every pairing the fluid model can build. B-dot appears only for the geochemical systems, + * the ones carrying ion size and b-dot parameters. + */ +using ReactionSystemList = std::variant< CarbonateIdentityEntry, + CarbonateBdotEntry, + CarbonateAllEquilibriumIdentityEntry, + CarbonateAllEquilibriumBdotEntry, + UltramaficIdentityEntry, + UltramaficBdotEntry, + MomasEasyIdentityEntry, + MomasMediumIdentityEntry, + ChainSerialIdentityEntry >; + +/** + * @brief Apply @p func to one instance of every system in ReactionSystemList. + * @tparam FUNC the callable type + * @tparam INDICES the alternatives of ReactionSystemList + * @param func a generic callable taking one entry + */ +template< typename FUNC, std::size_t ... INDICES > +void forEachReactionSystem( FUNC && func, std::index_sequence< INDICES... > ) +{ + ( func( std::variant_alternative_t< INDICES, ReactionSystemList >{} ), ... ); +} + +/** + * @copydoc forEachReactionSystem + */ +template< typename FUNC > +void forEachReactionSystem( FUNC && func ) +{ + forEachReactionSystem( std::forward< FUNC >( func ), + std::make_index_sequence< std::variant_size_v< ReactionSystemList > >{} ); +} + +/** + * @brief Whether @p chemicalSystem may be paired with @p activityModel. + * @param chemicalSystem the chemical system + * @param activityModel the activity model + * @return true if ReactionSystemList holds that pairing + */ +inline bool isSupportedReactionSystem( ChemicalSystemType const chemicalSystem, + ActivityModelType const activityModel ) +{ + bool supported = false; + forEachReactionSystem( [&]( auto system ) + { + using System = decltype( system ); + supported = supported || ( System::chemicalSystem == chemicalSystem && + System::activityModel == activityModel ); + } ); + return supported; +} + +} // namespace reactivefluid + +} // namespace constitutive + +} // namespace geos + +#endif // GEOS_CONSTITUTIVE_FLUID_REACTIVEFLUID_REACTIVEFLUIDSYSTEMSELECTOR_HPP_ diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index f497488a650..8341637f490 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -46,6 +46,14 @@ ReactiveSinglePhaseFluid( string const & name, Group * const parent ): setDescription( "Chemical System type. Available options are: " "``" + EnumStrings< ChemicalSystemType >::concat( "|" ) + "``" ); + this->registerWrapper( viewKeyStruct::activityModelNameString(), &m_activityModelType ). + setApplyDefaultValue( ActivityModelType::identity ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Activity model applied to the chemical system. Available options are: " + "``" + EnumStrings< ActivityModelType >::concat( "|" ) + "``. " + "``bdot`` requires ion size and b-dot parameters, which only the geochemical " + "systems carry." ); + this->registerWrapper( viewKeyStruct::solventMassPerSolutionVolumeString(), &m_solventMassPerSolutionVolume ). setApplyDefaultValue( 1000.0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -76,6 +84,7 @@ deliverClone( string const & name, Group * const parent ) const ReactiveSinglePhaseFluid & newConstitutiveRelation = dynamicCast< ReactiveSinglePhaseFluid & >( *clone ); newConstitutiveRelation.m_chemicalSystemType = m_chemicalSystemType; + newConstitutiveRelation.m_activityModelType = m_activityModelType; newConstitutiveRelation.m_numPrimarySpecies = m_numPrimarySpecies; newConstitutiveRelation.m_numSecondarySpecies = m_numSecondarySpecies; newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions; @@ -89,44 +98,25 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() { BASE::postInputInitialization(); - switch( m_chemicalSystemType ) + GEOS_THROW_IF( !isSupportedReactionSystem( m_chemicalSystemType, m_activityModelType ), + GEOS_FMT( "{}: chemical system '{}' has no '{}' activity model.", + this->getDataContext(), + EnumStrings< ChemicalSystemType >::toString( m_chemicalSystemType ), + EnumStrings< ActivityModelType >::toString( m_activityModelType ) ), + InputError ); + + // taken from the reaction system itself + forEachReactionSystem( [&]( auto system ) { - case ChemicalSystemType::ultramafic: - m_numPrimarySpecies = 9; - m_numSecondarySpecies = 16; - m_numKineticReactions = 5; - break; - - case ChemicalSystemType::carbonate: - m_numPrimarySpecies = 7; - m_numSecondarySpecies = 10; - m_numKineticReactions = 1; - break; - - case ChemicalSystemType::carbonateAllEquilibrium: - m_numPrimarySpecies = 7; - m_numSecondarySpecies = 11; - m_numKineticReactions = 0; - break; - - case ChemicalSystemType::chainSerialAllKinetic: - m_numPrimarySpecies = 3; - m_numSecondarySpecies = 0; - m_numKineticReactions = 3; - break; - - case ChemicalSystemType::momasMedium: - m_numPrimarySpecies = 5; - m_numSecondarySpecies = 9; - m_numKineticReactions = 1; - break; - - default: - m_numPrimarySpecies = 5; - m_numSecondarySpecies = 7; - m_numKineticReactions = 0; - break; - } + using System = decltype( system ); + if( System::chemicalSystem == m_chemicalSystemType && System::activityModel == m_activityModelType ) + { + using ReactionParams = typename System::ReactionParamsType; + m_numPrimarySpecies = ReactionParams::numPrimarySpecies(); + m_numSecondarySpecies = ReactionParams::numSecondarySpecies(); + m_numKineticReactions = ReactionParams::numKineticReactions(); + } + } ); GEOS_THROW_IF_LE_MSG( m_solventMassPerSolutionVolume, 0.0, GEOS_FMT( "invalid value of attribute '{}'", diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index 50ed70db6f1..6dbf99de1db 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -27,14 +27,11 @@ #include "constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp" #include "constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp" -#include "constitutive/HPCReact/src/reactions/geochemistry/GeochemicalSystems.hpp" -#include "constitutive/HPCReact/src/reactions/exampleSystems/BulkGeneric.hpp" -#include "constitutive/HPCReact/src/reactions/exampleSystems/ChainGeneric.hpp" -#include "constitutive/HPCReact/src/reactions/exampleSystems/MoMasBenchmark.hpp" +#include "constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp" #include "constitutive/HPCReact/src/reactions/reactionsSystems/EquilibriumReactions.hpp" #include "constitutive/HPCReact/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp" -#include "constitutive/HPCReact/src/reactions/massActions/MassActions.hpp" #include +#include namespace geos { @@ -47,16 +44,6 @@ namespace reactivefluid using namespace hpcReact::reactionsSystems; -enum class ChemicalSystemType : integer -{ - carbonate, - carbonateAllEquilibrium, - ultramafic, - momasEasy, - momasMedium, - chainSerialAllKinetic -}; - template< typename BASE > class ReactiveSinglePhaseFluid : public BASE { @@ -125,7 +112,7 @@ class ReactiveSinglePhaseFluid : public BASE /** * @brief Kernel wrapper class for ReactiveSinglePhaseFluid. */ - template< typename REACTION_PARAMS_TYPE > + template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > class ReactionKernelWrapper { @@ -143,7 +130,8 @@ class ReactiveSinglePhaseFluid : public BASE integer const numPrimarySpecies, integer const numSecondarySpecies, integer const numKineticReactions, - REACTION_PARAMS_TYPE params ): + REACTION_PARAMS_TYPE params, + typename ACTIVITY_MODEL::Params activityParams ): m_numPrimarySpecies( numPrimarySpecies ), m_numSecondarySpecies( numSecondarySpecies ), m_numKineticReactions( numKineticReactions ), @@ -156,10 +144,11 @@ class ReactiveSinglePhaseFluid : public BASE m_kineticReactionRates( kineticReactionRates ), m_aggregateSpeciesRates( aggregateSpeciesRates ), m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations( dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ), - m_params( params ) + m_params( params ), + m_activityParams( activityParams ) {} - using EquilibriumReactionsType = hpcReact::reactionsSystems::EquilibriumReactions< real64, integer, localIndex >; + using EquilibriumReactionsType = hpcReact::reactionsSystems::EquilibriumReactions< real64, integer, localIndex, ACTIVITY_MODEL >; /** * @brief Get number of elements in this wrapper. @@ -168,14 +157,22 @@ class ReactiveSinglePhaseFluid : public BASE GEOS_HOST_DEVICE localIndex numElems() const { return m_secondarySpeciesConcentration.size( 0 ); } + /** + * @brief Speciate cell @p k at equilibrium. + * @return whether the equilibrium solve converged + */ GEOS_HOST_DEVICE - void updateEquilibriumReaction( localIndex const k, + bool updateEquilibriumReaction( localIndex const k, real64 const pressure, real64 const temperature, arraySlice1d< real64, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration ) const; + /** + * @brief Solve for the primary and secondary concentrations at the target aggregates. + * @return whether the solve converged + */ GEOS_HOST_DEVICE - void enforceEquilibrium( real64 const pressure, + bool enforceEquilibrium( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & targetPrimarySpeciesAggregateConcentration, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & initialPrimarySpeciesConcentration, @@ -231,114 +228,55 @@ class ReactiveSinglePhaseFluid : public BASE arrayView4d< real64, reactivefluid::USD_SPECIES_DC > m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations; REACTION_PARAMS_TYPE m_params; + + typename ACTIVITY_MODEL::Params m_activityParams; }; - std::variant< - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::ultramaficSystemType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::carbonateSystemType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::carbonateSystemAllEquilibriumType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::ChainGeneric::serialAllKineticType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::mediumCaseType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::easyCaseType > > - createReactionKernelWrapper() const + /// The kernel wrapper for one system of reactivefluid::ReactionSystemList. + template< typename SYSTEM > + using WrapperFor = ReactionKernelWrapper< typename SYSTEM::ReactionParamsType, typename SYSTEM::ActivityType >; + + /// @cond DO_NOT_DOCUMENT + template< typename LIST > + struct WrapperVariantHelper; + + template< typename ... SYSTEMS > + struct WrapperVariantHelper< std::variant< SYSTEMS... > > { - using namespace hpcReact::geochemistry; - using namespace hpcReact::MoMasBenchmark; - using namespace hpcReact::bulkGeneric; - using namespace hpcReact::ChainGeneric; - switch( m_chemicalSystemType ) + using type = std::variant< WrapperFor< SYSTEMS > ... >; + }; + /// @endcond + + /// One alternative per system of reactivefluid::ReactionSystemList. + using ReactionKernelWrapperVariant = typename WrapperVariantHelper< reactivefluid::ReactionSystemList >::type; + + /** + * @brief Build the kernel wrapper for the chemical system and activity model this fluid was given. + * @return the wrapper, as the alternative of ReactionKernelWrapperVariant matching that pairing + * + * postInputInitialization has already rejected a pairing ReactionSystemList does not hold. + */ + ReactionKernelWrapperVariant createReactionKernelWrapper() const + { + std::optional< ReactionKernelWrapperVariant > wrapper; + + reactivefluid::forEachReactionSystem( [&]( auto system ) { - case ChemicalSystemType::ultramafic: - return ReactionKernelWrapper< ultramaficSystemType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - ultramaficSystem ); - - case ChemicalSystemType::carbonate: - return ReactionKernelWrapper< carbonateSystemType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - carbonateSystem ); - case ChemicalSystemType::carbonateAllEquilibrium: - return ReactionKernelWrapper< carbonateSystemAllEquilibriumType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - carbonateSystemAllEquilibrium ); - case ChemicalSystemType::chainSerialAllKinetic: - return ReactionKernelWrapper< serialAllKineticType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - serialAllKineticParams ); - case ChemicalSystemType::momasMedium: - return ReactionKernelWrapper< mediumCaseType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - mediumCaseParams ); - default: - return ReactionKernelWrapper< easyCaseType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - easyCaseParams ); - } + using System = decltype( system ); + if( System::chemicalSystem == m_chemicalSystemType && System::activityModel == m_activityModelType ) + { + wrapper.emplace( makeReactionKernelWrapper< WrapperFor< System > >( System::reactionParams(), + System::activityParams() ) ); + } + } ); + + return std::move( wrapper.value() ); } struct viewKeyStruct : ConstitutiveBase::viewKeyStruct { static constexpr char const * chemicalSystemNameString() { return "chemicalSystemType"; } + static constexpr char const * activityModelNameString() { return "activityModelType"; } static constexpr char const * solventMassPerSolutionVolumeString() { return "solventMassPerSolutionVolume"; } }; @@ -348,6 +286,29 @@ class ReactiveSinglePhaseFluid : public BASE virtual void resizeFields( localIndex const size, localIndex const numPts ); + /** + * @brief Build one kernel wrapper for the given reaction system and activity model. + */ + template< typename WRAPPER_TYPE, typename REACTION_PARAMS_TYPE, typename ACTIVITY_PARAMS_TYPE > + WRAPPER_TYPE makeReactionKernelWrapper( REACTION_PARAMS_TYPE const & params, + ACTIVITY_PARAMS_TYPE const & activityParams ) const + { + return WRAPPER_TYPE( m_primarySpeciesAggregateConcentration, + m_primarySpeciesMobileAggregateConcentration, + m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, + m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, + m_initialPrimarySpeciesConcentration, + m_secondarySpeciesConcentration, + m_kineticReactionRates, + m_aggregateSpeciesRates, + m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, + m_numPrimarySpecies, + m_numSecondarySpecies, + m_numKineticReactions, + params, + activityParams ); + } + integer m_numPrimarySpecies; integer m_numSecondarySpecies; @@ -376,6 +337,8 @@ class ReactiveSinglePhaseFluid : public BASE ChemicalSystemType m_chemicalSystemType; + ActivityModelType m_activityModelType; + /// TODO: prescribed as a constant for now. The exact factor is /// /// rho_s = rho * w @@ -394,41 +357,45 @@ using ReactiveCompressibleSinglePhaseFluid = ReactiveSinglePhaseFluid< Compressi using ReactiveThermalCompressibleSinglePhaseFluid = ReactiveSinglePhaseFluid< ThermalCompressibleSinglePhaseFluid >; template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE -inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +inline bool +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: updateEquilibriumReaction( localIndex const k, real64 const pressure, real64 const temperature, arraySlice1d< real64, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration ) const { - integer const numSecondarySpecies = m_numSecondarySpecies; + constexpr integer numSecondarySpecies = REACTION_PARAMS_TYPE::numSecondarySpecies(); - if( numSecondarySpecies > 0 ) + if constexpr( numSecondarySpecies > 0 ) { - stackArray1d< real64, MAX_NUM_SPECIES > logSecondarySpeciesConcentration( numSecondarySpecies ); + stackArray1d< real64, numSecondarySpecies > logSecondarySpeciesConcentration( numSecondarySpecies ); - enforceEquilibrium( pressure, temperature, m_primarySpeciesAggregateConcentration[k][0], m_initialPrimarySpeciesConcentration[k][0], logPrimarySpeciesConcentration, - logSecondarySpeciesConcentration.toSlice() ); + bool const converged = enforceEquilibrium( pressure, temperature, m_primarySpeciesAggregateConcentration[k][0], + m_initialPrimarySpeciesConcentration[k][0], logPrimarySpeciesConcentration, + logSecondarySpeciesConcentration.toSlice() ); for( integer i=0; i < numSecondarySpecies; ++i ) { m_secondarySpeciesConcentration[k][0][i] = LvArray::math::exp( logSecondarySpeciesConcentration[i] ); } + + return converged; } else { GEOS_UNUSED_VAR( k, pressure, temperature, logPrimarySpeciesConcentration ); + return true; } } template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE -inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +inline bool +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: enforceEquilibrium( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & targetPrimarySpeciesAggregateConcentration, @@ -449,20 +416,21 @@ enforceEquilibrium( real64 const pressure, logPrimarySpeciesConcentration0[i] = LvArray::math::log( initialPrimarySpeciesConcentration[i] ); } - // 1. We enforce equilibrium - EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature, m_params, targetPrimarySpeciesAggregateConc, logPrimarySpeciesConcentration0, logPrimarySpeciesConcentration ); - - // 2. We calculate the secondary species concentration - hpcReact::massActions::calculateLogSecondarySpeciesConcentration< real64, - localIndex, - localIndex >( m_params, logPrimarySpeciesConcentration, logSecondarySpeciesConcentration ); + // Solve for the primary and secondary concentrations with equilibrium enforced at the target aggregates + return EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature, + m_params, + m_activityParams, + targetPrimarySpeciesAggregateConc, + logPrimarySpeciesConcentration0, + logPrimarySpeciesConcentration, + logSecondarySpeciesConcentration ); } template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: updateMixedReactionSystem( localIndex const k, real64 const pressure, real64 const temperature, @@ -497,10 +465,10 @@ updateMixedReactionSystem( localIndex const k, } template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: computeAggregateConcentrationsAndRates( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration, @@ -517,9 +485,10 @@ computeAggregateConcentrationsAndRates( real64 const pressure, { GEOS_UNUSED_VAR( pressure ); - MixedEquilibriumKineticReactions< real64, localIndex, localIndex, true >:: + MixedEquilibriumKineticReactions< real64, localIndex, localIndex, ACTIVITY_MODEL, true >:: updateMixedSystem( temperature, m_params, + m_activityParams, logPrimarySpeciesConcentration, surfaceArea, logSecondarySpeciesConcentration, @@ -533,14 +502,6 @@ computeAggregateConcentrationsAndRates( real64 const pressure, dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); } -ENUM_STRINGS( ChemicalSystemType, - "carbonate", - "carbonateAllEquilibrium", - "ultramafic", - "momasEasy", - "momasMedium", - "chainSerialAllKinetic" ); - } // namespace reactivefluid } // namespace constitutive diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 97732290591..88702b535a1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -815,6 +815,8 @@ void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegi arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >(); arrayView2d< real64, compflow::USD_COMP > const logPrimaryConc = subRegion.getField< fields::flow::logPrimarySpeciesConcentration >(); + bool converged = true; + if( m_isThermal ) { reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid & fluid = @@ -822,7 +824,7 @@ void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegi constitutive::constitutiveUpdatePassThru( fluid, [&]( auto & castedFluid ) { - singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); + converged = singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); } ); fluid.saveConvergedState(); @@ -834,11 +836,15 @@ void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegi constitutive::constitutiveUpdatePassThru( fluid, [&]( auto & castedFluid ) { - singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); + converged = singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); } ); fluid.saveConvergedState(); } + + GEOS_ERROR_IF( !converged, + GEOS_FMT( "{}: the initial equilibrium speciation did not converge.", + subRegion.getDataContext() ) ); } void SinglePhaseReactiveTransport::initializePostInitialConditionsPreSubGroups() diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp index 79496d888fe..8d7ea5efb36 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp @@ -36,28 +36,45 @@ namespace singlePhaseReactiveBaseKernels struct EquilibriumReactionUpdateKernel { + /** + * @return whether the equilibrium solve converged in every cell + */ template< typename REACTION_WRAPPER_TYPE > - static void helper( REACTION_WRAPPER_TYPE const & reactionWrapper, + static bool helper( REACTION_WRAPPER_TYPE const & reactionWrapper, arrayView1d< real64 const > const & pres, arrayView1d< real64 const > const & temp, arrayView2d< real64, compflow::USD_COMP > const logPrimaryConc ) { + RAJA::ReduceMax< parallelDeviceReduce, integer > anyNotConverged( 0 ); + forAll< parallelDevicePolicy<> >( reactionWrapper.numElems(), [=] GEOS_HOST_DEVICE ( localIndex const k ) { - reactionWrapper.updateEquilibriumReaction( k, pres[k], temp[k], logPrimaryConc[k] ); + if( !reactionWrapper.updateEquilibriumReaction( k, pres[k], temp[k], logPrimaryConc[k] ) ) + { + anyNotConverged.max( 1 ); + } } ); + + return anyNotConverged.get() == 0; } + /** + * @return whether the equilibrium solve converged in every cell, for the caller to report + */ template< typename REACTIVE_FLUID > - static void launch( REACTIVE_FLUID const & fluid, + static bool launch( REACTIVE_FLUID const & fluid, arrayView1d< real64 const > const & pres, arrayView1d< real64 const > const & temp, arrayView2d< real64, compflow::USD_COMP > const logPrimaryConc ) { + bool converged = true; + std::visit( [&]( auto const & reactionWrapper ) { - helper( reactionWrapper, pres, temp, logPrimaryConc ); + converged = helper( reactionWrapper, pres, temp, logPrimaryConc ); }, fluid.createReactionKernelWrapper()); + + return converged; } }; diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 54cadeb69ab..09011b85355 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -8809,6 +8809,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + @@ -8840,6 +8842,11 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + + + + @@ -8894,6 +8901,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + From 174e0322210a2885db32586034072970c9e35a30 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 4 Sep 2026 21:24:49 -0700 Subject: [PATCH 08/15] applied solution scaling and fixed the normalizer --- .../MoMaSEasy_1DDiffusive_base.xml | 2 +- .../ReactiveSinglePhaseFluid.hpp | 2 +- .../SinglePhaseReactiveTransport.cpp | 57 +++++++++++ .../SinglePhaseReactiveTransport.hpp | 18 ++++ .../reactive/ResidualNormKernel.hpp | 8 +- .../reactive/SolutionScalingKernel.hpp | 94 +++++++++++++++++++ src/coreComponents/schema/schema.xsd | 2 + 7 files changed, 177 insertions(+), 6 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp diff --git a/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml b/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml index 07bc255f1fb..7a1387e0d5d 100644 --- a/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml +++ b/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml @@ -19,7 +19,7 @@ lineSearchAction="None" timeStepIncreaseFactor="1.2" logLevel="2" - maxAllowedResidualNorm="1e100"/> + maxAllowedResidualNorm="1e8"/> diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index 6dbf99de1db..ba73e35046c 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -368,7 +368,7 @@ updateEquilibriumReaction( localIndex const k, { constexpr integer numSecondarySpecies = REACTION_PARAMS_TYPE::numSecondarySpecies(); - if constexpr( numSecondarySpecies > 0 ) + if constexpr ( numSecondarySpecies > 0 ) { stackArray1d< real64, numSecondarySpecies > logSecondarySpeciesConcentration( numSecondarySpecies ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 88702b535a1..6ae48ef12ed 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -106,6 +106,13 @@ SinglePhaseReactiveTransport::SinglePhaseReactiveTransport( const string & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Array to store the indices of immobile species. Default is {}, which indicates no immobile species." ); + this->registerWrapper( viewKeyStruct::maxAbsoluteLogConcChangeString(), &m_maxAbsoluteLogConcChange ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 2.0 * 2.302585092994046 ). // two ln10 units + setDescription( "Maximum (absolute) change in the natural log of a primary species concentration " + "in a Newton iteration. Zero or less disables the scaling." ); + addLogLevel< logInfo::BoundaryConditions >(); } @@ -1355,6 +1362,56 @@ real64 SinglePhaseReactiveTransport::calculateResidualNorm( real64 const & GEOS_ return residualNorm; } +real64 SinglePhaseReactiveTransport::scalingForSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) +{ + GEOS_MARK_FUNCTION; + + // Pressure, and temperature when thermal, are scaled by the base. + real64 scalingFactor = SinglePhaseBase::scalingForSystemSolution( domain, dofManager, localSolution ); + + string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + integer const speciesOffset = m_isThermal ? 2 : 1; + real64 maxDeltaLogConc = 0.0; + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + auto const subRegionData = + singlePhaseReactiveBaseKernels::SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + dofNumber, + ghostRank, + m_numDofPerCell, + speciesOffset, + m_numPrimarySpecies, + m_maxAbsoluteLogConcChange ); + + scalingFactor = std::min( scalingFactor, subRegionData.first ); + maxDeltaLogConc = std::max( maxDeltaLogConc, subRegionData.second ); + } ); + } ); + + scalingFactor = MpiWrapper::min( scalingFactor ); + maxDeltaLogConc = MpiWrapper::max( maxDeltaLogConc ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, + GEOS_FMT( " {}: Max log concentration change = {:.4g} (before scaling)", + getName(), maxDeltaLogConc ) ); + + return scalingFactor; +} + void SinglePhaseReactiveTransport::applySystemSolution( DofManager const & dofManager, arrayView1d< real64 const > const & localSolution, real64 const scalingFactor, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp index 544b0f612f3..4fe5341ae1b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp @@ -28,6 +28,7 @@ #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluidUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp" @@ -127,6 +128,18 @@ class SinglePhaseReactiveTransport : public SinglePhaseBase DofManager const & dofManager, arrayView1d< real64 const > const & localRhs ) override; + /** + * @brief Compute the factor the Newton update is scaled by before it is applied. + * @param domain the physical domain + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localSolution the solution vector + * @return the smaller of the pressure scaling from the base and the log concentration scaling + */ + virtual real64 + scalingForSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) override; + virtual void applySystemSolution( DofManager const & dofManager, arrayView1d< real64 const > const & localSolution, @@ -271,10 +284,15 @@ class SinglePhaseReactiveTransport : public SinglePhaseBase static constexpr char const * isUpdateReactivePorosityString() { return "isUpdateReactivePorosity"; } static constexpr char const * isUpdateSurfaceAreaString() { return "isUpdateSurfaceArea"; } static constexpr char const * immobilePrimarySpeciesIndicesString() { return "immobilePrimarySpeciesIndices"; } + + static constexpr char const * maxAbsoluteLogConcChangeString() { return "maxAbsoluteLogConcentrationChange"; } }; protected: + /// Maximum (absolute) change in log primary species concentration in a Newton iteration, in ln units + real64 m_maxAbsoluteLogConcChange; + /// the number of primary species in the fluid integer m_numPrimarySpecies; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp index c9fc35acce1..9eccb735f01 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp @@ -77,7 +77,7 @@ class IsothermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNo // step 2: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); real64 const valAmount = LvArray::math::abs( m_localResidual[stack.localRow + idof + 1] ) / speciesAmountNormalizer; if( valAmount > stack.localValue[1] ) { @@ -99,7 +99,7 @@ class IsothermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNo // step 2: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); stack.localValue[1] += m_localResidual[stack.localRow + idof + 1] * m_localResidual[stack.localRow + idof + 1]; stack.localNormalizer[1] += speciesAmountNormalizer; @@ -185,7 +185,7 @@ class ThermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNormK // step 3: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); real64 const valAmount = LvArray::math::abs( m_localResidual[stack.localRow + idof + 2] ) / speciesAmountNormalizer; if( valAmount > stack.localValue[2] ) { @@ -214,7 +214,7 @@ class ThermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNormK // step 3: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); stack.localValue[2] += m_localResidual[stack.localRow + idof + 2] * m_localResidual[stack.localRow + idof + 2]; stack.localNormalizer[2] += speciesAmountNormalizer; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp new file mode 100644 index 00000000000..241b035df39 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp @@ -0,0 +1,94 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SolutionScalingKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_SOLUTIONSCALINGKERNEL_HPP +#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_SOLUTIONSCALINGKERNEL_HPP + +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" + +namespace geos +{ + +namespace singlePhaseReactiveBaseKernels +{ + +/** + * @brief Bound the change in log primary species concentration over a Newton iteration. + */ +struct SolutionScalingKernel +{ + /** + * @brief Compute the factor that keeps the largest log concentration change within the bound. + * @tparam POLICY the execution policy + * @param localSolution the solution vector + * @param rankOffset the offset of this rank in the solution vector + * @param dofNumber the first degree of freedom of each element + * @param ghostRank the ghost rank of each element + * @param numDofPerCell the number of degrees of freedom per cell + * @param speciesOffset the offset of the first species within a cell, 1 or 2 when thermal + * @param numPrimarySpecies the number of primary species + * @param maxAbsoluteLogConcChange the bound, in ln units. Zero or less disables the scaling. + * @return the scaling factor and the largest log concentration change seen, before scaling + */ + template< typename POLICY > + static std::pair< real64, real64 > launch( arrayView1d< real64 const > const & localSolution, + globalIndex const rankOffset, + arrayView1d< globalIndex const > const & dofNumber, + arrayView1d< integer const > const & ghostRank, + integer const numDofPerCell, + integer const speciesOffset, + integer const numPrimarySpecies, + real64 const maxAbsoluteLogConcChange ) + { + GEOS_UNUSED_VAR( numDofPerCell ); + + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > scalingFactor( 1.0 ); + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaLogConc( 0.0 ); + + forAll< POLICY >( dofNumber.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + if( ghostRank[ei] < 0 && dofNumber[ei] >= 0 ) + { + localIndex const lid = dofNumber[ei] - rankOffset; + + for( integer ic = 0; ic < numPrimarySpecies; ++ic ) + { + real64 const absLogConcChange = LvArray::math::abs( localSolution[lid + speciesOffset + ic] ); + maxDeltaLogConc.max( absLogConcChange ); + + if( maxAbsoluteLogConcChange > 0.0 && absLogConcChange > maxAbsoluteLogConcChange ) + { + scalingFactor.min( maxAbsoluteLogConcChange / absLogConcChange ); + } + } + } + + } ); + + return { scalingFactor.get(), maxDeltaLogConc.get() }; + } + +}; + +} // namespace singlePhaseReactiveBaseKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_SOLUTIONSCALINGKERNEL_HPP diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 09011b85355..86a6d9c96c5 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -5355,6 +5355,8 @@ Information output from lower logLevels is added with the desired log level 2 - The summary of declared fields and coupling--> + + From bb7c918ac149b269ec43006876fe298a76715495 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 4 Sep 2026 21:26:18 -0700 Subject: [PATCH 09/15] added smoke tests for bdot model and added log for equil. reaction run --- ...xedReactionCarbonateSystemBdot_1D_base.xml | 155 ++++++++++++++++++ ...edReactionCarbonateSystemBdot_1D_smoke.xml | 76 +++++++++ ...xedReactionCarbonateSystemBdot_2D_base.xml | 133 +++++++++++++++ ...edReactionCarbonateSystemBdot_2D_smoke.xml | 49 ++++++ .../singlePhaseReactiveTransport.ats | 17 ++ .../SinglePhaseReactiveTransport.cpp | 16 ++ 6 files changed, 446 insertions(+) create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_base.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_smoke.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_base.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_smoke.xml diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_base.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_base.xml new file mode 100644 index 00000000000..625440f19c2 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_base.xml @@ -0,0 +1,155 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_smoke.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_smoke.xml new file mode 100644 index 00000000000..67eede6b0f9 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_smoke.xml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_base.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_base.xml new file mode 100644 index 00000000000..2aca6139a30 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_base.xml @@ -0,0 +1,133 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_smoke.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_smoke.xml new file mode 100644 index 00000000000..9db70cb152d --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_smoke.xml @@ -0,0 +1,49 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats b/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats index 4f73e9e6fa3..7b38bb742f9 100644 --- a/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats +++ b/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats @@ -31,6 +31,23 @@ decks = [ restart_step=10, check_step=20, restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="mixedReactionCarbonateSystemBdot_1D_smoke", + description= + '1D Carbonate mixed reaction with the B-dot activity model (1D, mixedReaction, bdot)', + partitions=((1, 1, 1), (2, 1, 1)), + restart_step=10, + check_step=20, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + + TestDeck( + name="mixedReactionCarbonateSystemBdot_2D_smoke", + description= + '2D Carbonate mixed reaction with the B-dot activity model (2D, mixedReaction, bdot)', + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=20, + check_step=40, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), ] generate_geos_tests(decks) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 6ae48ef12ed..d7f8da65439 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -19,6 +19,7 @@ #include "SinglePhaseReactiveTransport.hpp" +#include "common/format/LogPart.hpp" #include "constitutive/ConstitutiveManager.hpp" #include "constitutive/ConstitutivePassThru.hpp" #include "constitutive/diffusion/DiffusionFields.hpp" @@ -772,6 +773,10 @@ void SinglePhaseReactiveTransport::updateSurfaceArea( ElementSubRegionBase & sub void SinglePhaseReactiveTransport::initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { + LogPart equilibriumLog( "Initial Chemical Equilibrium Enforcement", + MpiWrapper::commRank() == 0 && isLogLevelActive< logInfo::Convergence >( getLogLevel() ) ); + equilibriumLog.begin(); + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) { @@ -812,6 +817,8 @@ void SinglePhaseReactiveTransport::initializeFluidState( MeshLevel & mesh, strin diffusionMaterial.initializeTemperatureState( temperature ); } } ); + + equilibriumLog.end(); } void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegionBase & subRegion ) const @@ -849,6 +856,15 @@ void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegi fluid.saveConvergedState(); } + // Report the cell count and the status over all ranks, not just the one doing the logging. + globalIndex const numCells = MpiWrapper::sum< globalIndex >( subRegion.getNumberOfLocalIndices() ); + integer const allConverged = MpiWrapper::min< integer >( converged ? 1 : 0 ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, + GEOS_FMT( "{}: initial equilibrium speciation on {} ({} cells): {}", + getName(), subRegion.getName(), numCells, + allConverged ? "converged" : "NOT converged" ) ); + GEOS_ERROR_IF( !converged, GEOS_FMT( "{}: the initial equilibrium speciation did not converge.", subRegion.getDataContext() ) ); From 56111728045a41278ee1a0d2ad862a128776a510 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 4 Sep 2026 22:03:01 -0700 Subject: [PATCH 10/15] No plot for some field variables --- .../fluid/reactivefluid/ReactiveFluidFields.hpp | 8 ++++---- .../fluidFlow/SinglePhaseReactiveTransportFields.hpp | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp index 54e1c9a1fc2..81f2d6fb3d3 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp +++ b/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,7 +71,7 @@ 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" ); @@ -79,7 +79,7 @@ DECLARE_FIELD( dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesCon "dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations", array4dLayoutSpecies_dC, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Deivatives of primarySpeciesMobileAggregateConcentration w.r.t log primary species concentration" ); @@ -111,7 +111,7 @@ 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" ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp index 0d6a47880a0..4a0115f36b5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp @@ -47,7 +47,7 @@ DECLARE_FIELD( logPrimarySpeciesConcentration_n, "logPrimarySpeciesConcentration_n", array2dLayoutComp, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Natural log of primary species concentration (molality) at the previous converged time step" ); @@ -55,7 +55,7 @@ DECLARE_FIELD( bcLogPrimarySpeciesConcentration, "bcLogPrimarySpeciesConcentration", array2dLayoutComp, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Boundary condition for natural log of primary species concentration (molality)" ); @@ -71,8 +71,8 @@ DECLARE_FIELD( primarySpeciesAggregateMole_n, "primarySpeciesAggregateMole_n", array2dLayoutComp, 0, - LEVEL_0, - WRITE_AND_READ, + NOPLOT, + NO_WRITE, "Aggregate amount of primary species in mole at the previous converged time step" ); DECLARE_FIELD( kineticReactionMolarIncrements, From ba33803ba0a1bb651aa1f79cc46a1a8dbf545126 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sat, 5 Sep 2026 21:40:36 -0700 Subject: [PATCH 11/15] rebaseline --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 5dced8164ef..9f2a511c1d7 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3994-17525-4ae3593 + baseline: integratedTests/baseline_integratedTests-pr3996-17548-5611172 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 812f1523bcf..533bda15eeb 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -5,6 +5,10 @@ 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-05) +===================== +Convert species concentration fed into HPCReact from molarity to molality + PR #3994 (2026-09-03) ===================== Introduce mimetic inner products for saddle-point operators and align face-centroid computation with exact volume integration on planar polyhedra. From 0c87a3e764fe8210532b4eecff5a8955fa802c38 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sat, 5 Sep 2026 22:40:11 -0700 Subject: [PATCH 12/15] fixed unused var --- src/coreComponents/constitutive/HPCReact | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/constitutive/HPCReact b/src/coreComponents/constitutive/HPCReact index 02afd5d435e..e3acfa22211 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit 02afd5d435edcabcd7ac226b837b61038e5da0d3 +Subproject commit e3acfa22211162e9b9f3e651e8859ac9729fddf9 From 0a7813ef5c5105f00e4d60cc3b49541c54caa5c5 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sun, 6 Sep 2026 12:01:57 -0700 Subject: [PATCH 13/15] fixed CI tests --- .../fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp | 4 ++-- .../fluidFlowTests/testSinglePhaseReactiveTransport.cpp | 4 +++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index 8341637f490..0bfe027d519 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -50,8 +50,8 @@ ReactiveSinglePhaseFluid( string const & name, Group * const parent ): setApplyDefaultValue( ActivityModelType::identity ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "Activity model applied to the chemical system. Available options are: " - "``" + EnumStrings< ActivityModelType >::concat( "|" ) + "``. " - "``bdot`` requires ion size and b-dot parameters, which only the geochemical " + "``" + EnumStrings< ActivityModelType >::concat( "|" ) + + ". ``bdot`` requires ion size and b-dot parameters, which only the geochemical " "systems carry." ); this->registerWrapper( viewKeyStruct::solventMassPerSolutionVolumeString(), &m_solventMassPerSolutionVolume ). diff --git a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp index c4b50dfab30..0875378f4c8 100644 --- a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp @@ -562,13 +562,15 @@ TEST_F( SinglePhaseReactiveTransportTest, jacobianNumericalCheck_accumulationBal real64 const perturb = std::sqrt( eps ); real64 const tol = 1e-6; // 1% error margin + real64 const accumulationDt = 1e-4; + DomainPartition & domain = state.getProblemManager().getDomainPartition(); testNumericalJacobian( *solver, domain, perturb, tol, [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - solver->assembleAccumulationTermsInMassBalanceAndSpeciesAmountEqs( dt, domain, solver->getDofManager(), localMatrix, localRhs ); + solver->assembleAccumulationTermsInMassBalanceAndSpeciesAmountEqs( accumulationDt, domain, solver->getDofManager(), localMatrix, localRhs ); } ); } From 9fb8b70c541329bda74b4932811c273d7f0936ca Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Mon, 14 Sep 2026 21:56:52 -0700 Subject: [PATCH 14/15] fixed 1D chain reaction run and added to integrated test --- .../ChainSerial_1DAllKinetic_base.xml | 4 +- .../ChainSerial_1DAllKinetic_benchmark.xml | 9 ++-- .../ChainSerial_1DAllKinetic_smoke.xml | 43 +++++++++++++++++++ .../ChainSerial_initialLogConc.xml | 33 ++++++++++++++ .../singlePhaseReactiveTransport.ats | 9 ++++ 5 files changed, 92 insertions(+), 6 deletions(-) create mode 100644 inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/ChainSerial_initialLogConc.xml diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml index c58004bb4f4..3982548312e 100644 --- a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml +++ b/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"> + maxAllowedResidualNorm="1e10"/> diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml index 355a4289718..e75583f52dd 100644 --- a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml @@ -5,6 +5,7 @@ + @@ -25,17 +26,17 @@ maxTime="400"> diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml new file mode 100644 index 00000000000..4ebc3d6c238 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml @@ -0,0 +1,43 @@ + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialLogConc.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialLogConc.xml new file mode 100644 index 00000000000..1368324cd18 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialLogConc.xml @@ -0,0 +1,33 @@ + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats b/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats index 4f73e9e6fa3..88bff398cfe 100644 --- a/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats +++ b/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= From 59c712959d2babf74775e6ba4be060d659096a51 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Tue, 15 Sep 2026 20:34:47 -0700 Subject: [PATCH 15/15] enable additional equilibrium enforce constraint options --- .../ChainSerial_1DAllKinetic_benchmark.xml | 2 +- .../ChainSerial_1DAllKinetic_smoke.xml | 2 +- .../ChainSerial_initialAggregate.xml | 12 +- .../ChainSerial_initialPrimaryConc.xml | 61 +++++++ .../MoMaSEasy_initialAggregate.xml | 30 ++-- .../MoMaSMedium_initialAggregate.xml | 30 ++-- ...rbonateSystemBdotChargeBalance_1D_base.xml | 156 ++++++++++++++++++ ...bonateSystemBdotChargeBalance_1D_smoke.xml | 76 +++++++++ ...rbonateSystemBdotChargeBalance_logConc.xml | 62 +++++++ ...dReactionCarbonateSystemBdotPX_1D_base.xml | 156 ++++++++++++++++++ ...ReactionCarbonateSystemBdotPX_1D_smoke.xml | 76 +++++++++ ...nCarbonateSystemBdotPX_constraintValue.xml | 133 +++++++++++++++ ...actionCarbonateSystem_initialAggregate.xml | 28 ++-- src/coreComponents/constitutive/HPCReact | 2 +- .../reactivefluid/ReactiveFluidFields.hpp | 8 + .../ReactiveFluidSystemSelector.hpp | 28 ++++ .../ReactiveSinglePhaseFluid.cpp | 90 ++++++++++ .../ReactiveSinglePhaseFluid.hpp | 109 ++++++++---- .../testSinglePhaseReactiveTransport.cpp | 28 ++-- 19 files changed, 987 insertions(+), 102 deletions(-) create mode 100644 inputFiles/singlePhaseReactiveTransport/ChainSerial_initialPrimaryConc.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_base.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_smoke.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_logConc.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_base.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_smoke.xml create mode 100644 inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_constraintValue.xml diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml index e75583f52dd..a49b3a96bdb 100644 --- a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_benchmark.xml @@ -5,7 +5,7 @@ - + diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml index 4ebc3d6c238..d6b6f52f387 100644 --- a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_smoke.xml @@ -5,7 +5,7 @@ - + diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialAggregate.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialAggregate.xml index a72597cb050..f60e6a75525 100644 --- a/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialAggregate.xml +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialAggregate.xml @@ -8,7 +8,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Medium" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1e-8"/> @@ -17,7 +17,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Medium" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="1e-8"/> @@ -26,7 +26,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Medium" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1e-8"/> @@ -36,7 +36,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1.0"/> @@ -45,7 +45,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="1e-8"/> @@ -54,7 +54,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1e-8"/> diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialPrimaryConc.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialPrimaryConc.xml new file mode 100644 index 00000000000..512dda6439a --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_initialPrimaryConc.xml @@ -0,0 +1,61 @@ + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_initialAggregate.xml b/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_initialAggregate.xml index d17545e2d17..82d23359027 100644 --- a/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_initialAggregate.xml +++ b/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_initialAggregate.xml @@ -8,7 +8,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1e-8"/> @@ -17,7 +17,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="-2.0"/> @@ -26,7 +26,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1e-8"/> @@ -35,7 +35,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="2.0"/> @@ -44,7 +44,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="1.0"/> @@ -54,7 +54,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1e-8"/> @@ -63,7 +63,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="-2.0"/> @@ -72,7 +72,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1e-8"/> @@ -81,7 +81,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="2.0"/> @@ -90,7 +90,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="10.0"/> @@ -100,7 +100,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="0.3"/> @@ -109,7 +109,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="0.3"/> @@ -118,7 +118,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="0.3"/> @@ -127,7 +127,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="1e-20"/> @@ -136,7 +136,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="1e-20"/> diff --git a/inputFiles/singlePhaseReactiveTransport/MoMaSMedium_initialAggregate.xml b/inputFiles/singlePhaseReactiveTransport/MoMaSMedium_initialAggregate.xml index 2fb76a9a089..950591742b6 100644 --- a/inputFiles/singlePhaseReactiveTransport/MoMaSMedium_initialAggregate.xml +++ b/inputFiles/singlePhaseReactiveTransport/MoMaSMedium_initialAggregate.xml @@ -8,7 +8,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1e-8"/> @@ -17,7 +17,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="-3.0"/> @@ -26,7 +26,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1e-8"/> @@ -35,7 +35,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="1.0"/> @@ -44,7 +44,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumA" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="1.0"/> @@ -54,7 +54,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1e-8"/> @@ -63,7 +63,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="-9.0"/> @@ -72,7 +72,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1e-8"/> @@ -81,7 +81,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="3.0"/> @@ -90,7 +90,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/MediumB" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="10.0"/> @@ -100,7 +100,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="0.3"/> @@ -109,7 +109,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="0.3"/> @@ -118,7 +118,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="0.3"/> @@ -127,7 +127,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="1e-8"/> @@ -136,7 +136,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Source" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="1e-8"/> diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_base.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_base.xml new file mode 100644 index 00000000000..a9fe333a564 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_base.xml @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_smoke.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_smoke.xml new file mode 100644 index 00000000000..1fa8b2f5deb --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_1D_smoke.xml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_logConc.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_logConc.xml new file mode 100644 index 00000000000..643220ff222 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotChargeBalance_logConc.xml @@ -0,0 +1,62 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_base.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_base.xml new file mode 100644 index 00000000000..40c35c7f027 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_base.xml @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_smoke.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_smoke.xml new file mode 100644 index 00000000000..6345c9ecb49 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_1D_smoke.xml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_constraintValue.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_constraintValue.xml new file mode 100644 index 00000000000..2163b6fe929 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdotPX_constraintValue.xml @@ -0,0 +1,133 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystem_initialAggregate.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystem_initialAggregate.xml index 767b9387978..34965b497e0 100644 --- a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystem_initialAggregate.xml +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystem_initialAggregate.xml @@ -8,7 +8,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1.585e-7"/> @@ -17,7 +17,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="8.293e-3"/> @@ -26,7 +26,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="2.171e-3"/> @@ -35,7 +35,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="1.666e-6"/> @@ -44,7 +44,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="2.821e-5"/> @@ -53,7 +53,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="5" scale="1.605e-3"/> @@ -62,7 +62,7 @@ initialCondition="1" setNames="{ all }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="6" scale="7.817e-4"/> @@ -72,7 +72,7 @@ initialCondition="1" setNames="{ source }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1.585e-7"/> @@ -81,7 +81,7 @@ initialCondition="1" setNames="{ source }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="7.317e-2"/> @@ -90,7 +90,7 @@ initialCondition="1" setNames="{ source }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1.517e-2"/> @@ -99,7 +99,7 @@ initialCondition="1" setNames="{ source }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="1.666e-6"/> @@ -108,7 +108,7 @@ initialCondition="1" setNames="{ source }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="2.821e-5"/> @@ -117,7 +117,7 @@ initialCondition="1" setNames="{ source }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="5" scale="1.605e-3"/> @@ -126,7 +126,7 @@ initialCondition="1" setNames="{ source }" objectPath="ElementRegions/Region1/block1" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="6" scale="7.817e-4"/> diff --git a/src/coreComponents/constitutive/HPCReact b/src/coreComponents/constitutive/HPCReact index e3acfa22211..be4f3b86354 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit e3acfa22211162e9b9f3e651e8859ac9729fddf9 +Subproject commit be4f3b863542840b271a679a92c4a7255db181b9 diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp index 81f2d6fb3d3..08c9e5f760a 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp @@ -43,6 +43,14 @@ DECLARE_FIELD( initialPrimarySpeciesConcentration, WRITE_AND_READ, "initialPrimarySpeciesConcentration" ); +DECLARE_FIELD( primarySpeciesConstraintValue, + "primarySpeciesConstraintValue", + array3dLayoutSpecies, + 1e-16, + NOPLOT, + WRITE_AND_READ, + "Value of the constraint closing each primary species' row of the initial equilibrium solve" ); + DECLARE_FIELD( primarySpeciesAggregateConcentration, "primarySpeciesAggregateConcentration", array3dLayoutSpecies, diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp index a70e31407ed..baaee174a42 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp @@ -36,6 +36,23 @@ #include #include +namespace hpcReact +{ +namespace reactionsSystems +{ + +using geos::EnumStrings; +using geos::string; + +ENUM_STRINGS( PrimarySpeciesConstraintType, + "aggregateConcentration", + "pX", + "chargeBalance", + "mineralEquilibrium" ); + +} // namespace reactionsSystems +} // namespace hpcReact + namespace geos { @@ -81,6 +98,17 @@ ENUM_STRINGS( ActivityModelType, "identity", "bdot" ); +/** + * @brief The constraint that closes one primary species' row of the initial equilibrium solve. + * @details XML names: + * - ``aggregateConcentration`` (default): the aggregate concentration supplied by the input. + * - ``pX``: the negative log of the species' activity, i.e. pH when the species is H+. + * - ``chargeBalance``: replaces the aggregate concentration with electroneutrality. At most + * one species may carry it. + * - ``mineralEquilibrium``: reserved and rejected, HPCReact does not implement it for now. + */ +using PrimarySpeciesConstraintType = hpcReact::reactionsSystems::PrimarySpeciesConstraintType; + /** * @brief What every entry below shares: the two HPCReact types, and the two enum values the XML * input selects it by. diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index 0bfe027d519..5f9b1d5cc1b 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -61,6 +61,18 @@ ReactiveSinglePhaseFluid( string const & name, Group * const parent ): "molality [mol/kg solvent] to molarity [mol/m^3 solution]. The default of 1000 " "approximates an aqueous solution by the density of pure water." ); + this->registerWrapper( viewKeyStruct::primarySpeciesConstraintTypesString(), &m_primarySpeciesConstraintTypeInput ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Constraint closing each primary species' row of the initial equilibrium solve, " + "one entry per primary species in declaration order. Available options are: " + "``" + EnumStrings< PrimarySpeciesConstraintType >::concat( "|" ) + + "``. Defaults to ``aggregateConcentration`` for every species when unset. Each " + "species reads its entry of the ``primarySpeciesConstraintValue`` field, set as " + "``_primarySpeciesConstraintValue``, in the units its constraint type " + "implies, and the solve reports the aggregate concentration of every species " + "not constrained by one." ); + + this->template registerField< fields::reactivefluid::primarySpeciesConstraintValue >( &m_primarySpeciesConstraintValue ); this->template registerField< fields::reactivefluid::initialPrimarySpeciesConcentration >( &m_initialPrimarySpeciesConcentration ); this->template registerField< fields::reactivefluid::secondarySpeciesConcentration >( &m_secondarySpeciesConcentration ); this->template registerField< fields::reactivefluid::primarySpeciesAggregateConcentration >( &m_primarySpeciesAggregateConcentration ); @@ -89,6 +101,8 @@ deliverClone( string const & name, Group * const parent ) const newConstitutiveRelation.m_numSecondarySpecies = m_numSecondarySpecies; newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions; newConstitutiveRelation.m_solventMassPerSolutionVolume = m_solventMassPerSolutionVolume; + newConstitutiveRelation.m_primarySpeciesConstraintTypeInput = m_primarySpeciesConstraintTypeInput; + newConstitutiveRelation.m_primarySpeciesConstraintType = m_primarySpeciesConstraintType; return clone; } @@ -122,6 +136,81 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() GEOS_FMT( "invalid value of attribute '{}'", viewKeyStruct::solventMassPerSolutionVolumeString() ), InputError, this->getDataContext() ); + + checkPrimarySpeciesConstraints(); +} + +template< typename BASE > +void ReactiveSinglePhaseFluid< BASE >::checkPrimarySpeciesConstraints() +{ + string const constraintKey = viewKeyStruct::primarySpeciesConstraintTypesString(); + + m_primarySpeciesConstraintType.resize( m_numPrimarySpecies ); + + // If no constraint types are given, constrain every species by its aggregate concentration. + if( m_primarySpeciesConstraintTypeInput.empty() ) + { + for( integer i = 0; i < m_numPrimarySpecies; ++i ) + { + m_primarySpeciesConstraintType[i] = static_cast< integer >( PrimarySpeciesConstraintType::AggregateConcentration ); + } + return; + } + + GEOS_THROW_IF_NE_MSG( LvArray::integerConversion< integer >( m_primarySpeciesConstraintTypeInput.size() ), + m_numPrimarySpecies, + GEOS_FMT( "{}: '{}' must have one entry per primary species. Chemical system " + "'{}' has {} of them.", + this->getDataContext(), constraintKey, + EnumStrings< ChemicalSystemType >::toString( m_chemicalSystemType ), + m_numPrimarySpecies ), + InputError ); + + integer chargeBalanceCount = 0; + + for( integer i = 0; i < m_numPrimarySpecies; ++i ) + { + PrimarySpeciesConstraintType const type = + EnumStrings< PrimarySpeciesConstraintType >::fromString( m_primarySpeciesConstraintTypeInput[i] ); + + GEOS_THROW_IF( type == PrimarySpeciesConstraintType::MineralEquilibrium, + GEOS_FMT( "{}: '{}' entry {} requests '{}', which is not implemented yet.", + this->getDataContext(), constraintKey, i, + EnumStrings< PrimarySpeciesConstraintType >::toString( type ) ), + InputError ); + + if( type == PrimarySpeciesConstraintType::ChargeBalance ) + { + // A neutral species does not appear in the charge balance, so it cannot be solved for by it. + real64 charge = 0.0; + forEachReactionSystem( [&]( auto system ) + { + using System = decltype( system ); + if( System::chemicalSystem == m_chemicalSystemType && System::activityModel == m_activityModelType ) + { + charge = System::activityParams().m_speciesCharge[i + m_numSecondarySpecies]; + } + } ); + + GEOS_THROW_IF( charge == 0.0, + GEOS_FMT( "{}: '{}' entry {} requests '{}' on a neutral species.", + this->getDataContext(), constraintKey, i, + EnumStrings< PrimarySpeciesConstraintType >::toString( type ) ), + InputError ); + + ++chargeBalanceCount; + } + + m_primarySpeciesConstraintType[i] = static_cast< integer >( type ); + } + + // Electroneutrality is one equation; more than one species carrying it makes the solve singular. + GEOS_THROW_IF_GT_MSG( chargeBalanceCount, 1, + GEOS_FMT( "{}: {} primary species request '{}'. At most one species may carry it.", + this->getDataContext(), chargeBalanceCount, + EnumStrings< PrimarySpeciesConstraintType >::toString( + PrimarySpeciesConstraintType::ChargeBalance ) ), + InputError ); } template< typename BASE > @@ -139,6 +228,7 @@ void ReactiveSinglePhaseFluid< BASE >::resizeFields( localIndex const size, loca integer const numSecondarySpecies = this->numSecondarySpecies(); integer const numKineticReactions = this->numKineticReactions(); + m_primarySpeciesConstraintValue.resize( size, numPts, numPrimarySpecies ); m_initialPrimarySpeciesConcentration.resize( size, numPts, numPrimarySpecies ); m_secondarySpeciesConcentration.resize( size, numPts, numSecondarySpecies ); m_primarySpeciesAggregateConcentration.resize( size, numPts, numPrimarySpecies ); diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index ba73e35046c..9b3884de3f4 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -67,6 +67,9 @@ class ReactiveSinglePhaseFluid : public BASE static constexpr integer MAX_NUM_SPECIES = 20; static constexpr integer MAX_NUM_KINETIC_REACTIONS = 10; + arrayView3d< real64 const, reactivefluid::USD_SPECIES > primarySpeciesConstraintValue() const + { return m_primarySpeciesConstraintValue; } + arrayView3d< real64 const, reactivefluid::USD_SPECIES > primarySpeciesAggregateConcentration() const { return m_primarySpeciesAggregateConcentration; } @@ -122,16 +125,19 @@ class ReactiveSinglePhaseFluid : public BASE arrayView3d< real64, reactivefluid::USD_SPECIES > const & primarySpeciesMobileAggregateConcentration, arrayView4d< real64, reactivefluid::USD_SPECIES_DC > const & dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, arrayView4d< real64, reactivefluid::USD_SPECIES_DC > const & dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & primarySpeciesConstraintValue, arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & initialPrimarySpeciesConcentration, arrayView3d< real64, reactivefluid::USD_SPECIES > const & secondarySpeciesConcentration, arrayView3d< real64, reactivefluid::USD_SPECIES > const & kineticReactionRates, arrayView3d< real64, reactivefluid::USD_SPECIES > const & aggregateSpeciesRates, arrayView4d< real64, reactivefluid::USD_SPECIES_DC > const & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, + arrayView1d< integer const > const & primarySpeciesConstraintType, integer const numPrimarySpecies, integer const numSecondarySpecies, integer const numKineticReactions, REACTION_PARAMS_TYPE params, typename ACTIVITY_MODEL::Params activityParams ): + m_primarySpeciesConstraintType( primarySpeciesConstraintType ), m_numPrimarySpecies( numPrimarySpecies ), m_numSecondarySpecies( numSecondarySpecies ), m_numKineticReactions( numKineticReactions ), @@ -139,6 +145,7 @@ class ReactiveSinglePhaseFluid : public BASE m_primarySpeciesMobileAggregateConcentration( primarySpeciesMobileAggregateConcentration ), m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations( dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations ), m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations( dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations ), + m_primarySpeciesConstraintValue( primarySpeciesConstraintValue ), m_initialPrimarySpeciesConcentration( initialPrimarySpeciesConcentration ), m_secondarySpeciesConcentration( secondarySpeciesConcentration ), m_kineticReactionRates( kineticReactionRates ), @@ -174,10 +181,11 @@ class ReactiveSinglePhaseFluid : public BASE GEOS_HOST_DEVICE bool enforceEquilibrium( real64 const pressure, real64 const temperature, - arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & targetPrimarySpeciesAggregateConcentration, + arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & primarySpeciesConstraintValue, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & initialPrimarySpeciesConcentration, arraySlice1d< real64, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration, - arraySlice1d< real64 > const & logSecondarySpeciesConcentration ) const; + arraySlice1d< real64 > const & logSecondarySpeciesConcentration, + arraySlice1d< real64, reactivefluid::USD_SPECIES - 2 > const & primarySpeciesAggregateConcentration ) const; GEOS_HOST_DEVICE void updateMixedReactionSystem( localIndex const k, @@ -203,6 +211,8 @@ class ReactiveSinglePhaseFluid : public BASE protected: + arrayView1d< integer const > m_primarySpeciesConstraintType; + integer m_numPrimarySpecies; integer m_numSecondarySpecies; @@ -217,6 +227,8 @@ class ReactiveSinglePhaseFluid : public BASE arrayView4d< real64, reactivefluid::USD_SPECIES_DC > m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations; + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const m_primarySpeciesConstraintValue; + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const m_initialPrimarySpeciesConcentration; arrayView3d< real64, reactivefluid::USD_SPECIES > m_secondarySpeciesConcentration; @@ -278,12 +290,20 @@ class ReactiveSinglePhaseFluid : public BASE static constexpr char const * chemicalSystemNameString() { return "chemicalSystemType"; } static constexpr char const * activityModelNameString() { return "activityModelType"; } static constexpr char const * solventMassPerSolutionVolumeString() { return "solventMassPerSolutionVolume"; } + static constexpr char const * primarySpeciesConstraintTypesString() { return "primarySpeciesConstraintTypes"; } }; protected: virtual void postInputInitialization() override; + /** + * @brief Check the input constraint types and store them in the array the solve takes. + * @details Rejects a length that does not match the species count, mineralEquilibrium, and more + * than one chargeBalance. + */ + void checkPrimarySpeciesConstraints(); + virtual void resizeFields( localIndex const size, localIndex const numPts ); /** @@ -297,11 +317,13 @@ class ReactiveSinglePhaseFluid : public BASE m_primarySpeciesMobileAggregateConcentration, m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, + m_primarySpeciesConstraintValue, m_initialPrimarySpeciesConcentration, m_secondarySpeciesConcentration, m_kineticReactionRates, m_aggregateSpeciesRates, m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, + m_primarySpeciesConstraintType.toViewConst(), m_numPrimarySpecies, m_numSecondarySpecies, m_numKineticReactions, @@ -317,6 +339,9 @@ class ReactiveSinglePhaseFluid : public BASE array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_initialPrimarySpeciesConcentration; + /// Value of the constraint each primary species carries, in the units its constraint type reads. + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_primarySpeciesConstraintValue; + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_secondarySpeciesConcentration; array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_primarySpeciesAggregateConcentration; @@ -339,6 +364,15 @@ class ReactiveSinglePhaseFluid : public BASE ActivityModelType m_activityModelType; + /// Constraint closing each primary species' row of the initial equilibrium solve, as named in the + /// input file, in the species order of the chemical system. Empty means every species is + /// constrained by its total concentration. + string_array m_primarySpeciesConstraintTypeInput; + + /// m_primarySpeciesConstraintTypeInput resolved to PrimarySpeciesConstraintType values, always + /// numPrimarySpecies long. Not an input. + array1d< integer > m_primarySpeciesConstraintType; + /// TODO: prescribed as a constant for now. The exact factor is /// /// rho_s = rho * w @@ -367,28 +401,23 @@ updateEquilibriumReaction( localIndex const k, arraySlice1d< real64, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration ) const { constexpr integer numSecondarySpecies = REACTION_PARAMS_TYPE::numSecondarySpecies(); + // A stack array needs a capacity of at least one, even when there are no secondary species. + constexpr integer numSecondarySpeciesStorage = numSecondarySpecies > 0 ? numSecondarySpecies : 1; - if constexpr ( numSecondarySpecies > 0 ) - { - stackArray1d< real64, numSecondarySpecies > logSecondarySpeciesConcentration( numSecondarySpecies ); + stackArray1d< real64, numSecondarySpeciesStorage > logSecondarySpeciesConcentration( numSecondarySpecies ); - bool const converged = enforceEquilibrium( pressure, temperature, m_primarySpeciesAggregateConcentration[k][0], - m_initialPrimarySpeciesConcentration[k][0], logPrimarySpeciesConcentration, - logSecondarySpeciesConcentration.toSlice() ); + bool const converged = enforceEquilibrium( pressure, temperature, m_primarySpeciesConstraintValue[k][0], + m_initialPrimarySpeciesConcentration[k][0], logPrimarySpeciesConcentration, + logSecondarySpeciesConcentration.toSlice(), + m_primarySpeciesAggregateConcentration[k][0] ); - for( integer i=0; i < numSecondarySpecies; ++i ) - { - m_secondarySpeciesConcentration[k][0][i] = LvArray::math::exp( logSecondarySpeciesConcentration[i] ); - } - - return converged; - } - else + for( integer i=0; i < numSecondarySpecies; ++i ) { - GEOS_UNUSED_VAR( k, pressure, temperature, logPrimarySpeciesConcentration ); - return true; + m_secondarySpeciesConcentration[k][0][i] = LvArray::math::exp( logSecondarySpeciesConcentration[i] ); } + return converged; + } template< typename BASE > @@ -398,32 +427,39 @@ inline bool ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: enforceEquilibrium( real64 const pressure, real64 const temperature, - arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & targetPrimarySpeciesAggregateConcentration, + arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & primarySpeciesConstraintValue, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & initialPrimarySpeciesConcentration, arraySlice1d< real64, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration, - arraySlice1d< real64 > const & logSecondarySpeciesConcentration ) const + arraySlice1d< real64 > const & logSecondarySpeciesConcentration, + arraySlice1d< real64, reactivefluid::USD_SPECIES - 2 > const & primarySpeciesAggregateConcentration ) const { GEOS_UNUSED_VAR( pressure ); integer const numPrimarySpecies = m_numPrimarySpecies; stackArray1d< real64, MAX_NUM_SPECIES > logPrimarySpeciesConcentration0( numPrimarySpecies ); - stackArray1d< real64, MAX_NUM_SPECIES > targetPrimarySpeciesAggregateConc( numPrimarySpecies ); + stackArray1d< real64, MAX_NUM_SPECIES > constraintValue( numPrimarySpecies ); + + using ConstraintType = PrimarySpeciesConstraintType; + stackArray1d< ConstraintType, MAX_NUM_SPECIES > constraintType( numPrimarySpecies ); for( integer i=0; i < numPrimarySpecies; ++i ) { - targetPrimarySpeciesAggregateConc[i] = targetPrimarySpeciesAggregateConcentration[i]; + constraintValue[i] = primarySpeciesConstraintValue[i]; logPrimarySpeciesConcentration0[i] = LvArray::math::log( initialPrimarySpeciesConcentration[i] ); + constraintType[i] = static_cast< ConstraintType >( m_primarySpeciesConstraintType[i] ); } - // Solve for the primary and secondary concentrations with equilibrium enforced at the target aggregates - return EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature, - m_params, - m_activityParams, - targetPrimarySpeciesAggregateConc, - logPrimarySpeciesConcentration0, - logPrimarySpeciesConcentration, - logSecondarySpeciesConcentration ); + // Solve for the primary and secondary concentrations under the constraint each species carries. + return EquilibriumReactionsType::enforceEquilibrium_PrimaryConcentrations( temperature, + m_params, + m_activityParams, + constraintType, + constraintValue, + logPrimarySpeciesConcentration0, + logPrimarySpeciesConcentration, + logSecondarySpeciesConcentration, + primarySpeciesAggregateConcentration ); } template< typename BASE > @@ -437,12 +473,15 @@ updateMixedReactionSystem( localIndex const k, arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration, arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & surfaceArea ) const { - integer const numPrimarySpecies = m_numPrimarySpecies; - integer const numSecondarySpecies = m_numSecondarySpecies; - integer const numKineticReactions = m_numKineticReactions; + constexpr integer numPrimarySpecies = REACTION_PARAMS_TYPE::numPrimarySpecies(); + constexpr integer numSecondarySpecies = REACTION_PARAMS_TYPE::numSecondarySpecies(); + constexpr integer numKineticReactions = REACTION_PARAMS_TYPE::numKineticReactions(); + // A stack array needs a capacity of at least one, even when the system has none of these. + constexpr integer numSecondarySpeciesStorage = numSecondarySpecies > 0 ? numSecondarySpecies : 1; + constexpr integer dReactionRatesStorage = numKineticReactions * numPrimarySpecies > 0 ? numKineticReactions * numPrimarySpecies : 1; - stackArray1d< real64, MAX_NUM_SPECIES > logSecondarySpeciesConcentration( numSecondarySpecies ); - stackArray2d< real64, MAX_NUM_KINETIC_REACTIONS * MAX_NUM_SPECIES > dReactionRates_dLogPrimarySpeciesConcentrations( numKineticReactions, numPrimarySpecies ); + stackArray1d< real64, numSecondarySpeciesStorage > logSecondarySpeciesConcentration( numSecondarySpecies ); + stackArray2d< real64, dReactionRatesStorage > dReactionRates_dLogPrimarySpeciesConcentrations( numKineticReactions, numPrimarySpecies ); computeAggregateConcentrationsAndRates( pressure, temperature, diff --git a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp index 0875378f4c8..04c44e24765 100644 --- a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp @@ -147,7 +147,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ all }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1.585e-7"/> @@ -156,7 +156,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ all }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="8.293e-3"/> @@ -165,7 +165,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ all }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="2.171e-3"/> @@ -174,7 +174,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ all }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="1.666e-6"/> @@ -183,7 +183,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ all }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="2.821e-5"/> @@ -192,7 +192,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ all }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="5" scale="1.605e-3"/> @@ -201,7 +201,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ all }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="6" scale="7.817e-4"/> @@ -210,7 +210,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ source }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="0" scale="1.585e-7"/> @@ -219,7 +219,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ source }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="1" scale="7.317e-2"/> @@ -228,7 +228,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ source }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="2" scale="1.517e-2"/> @@ -237,7 +237,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ source }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="3" scale="1.666e-6"/> @@ -246,7 +246,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ source }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="4" scale="2.821e-5"/> @@ -255,7 +255,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ source }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="5" scale="1.605e-3"/> @@ -264,7 +264,7 @@ char const * xmlInputCarbonate = initialCondition="1" setNames="{ source }" objectPath="ElementRegions/region/cb" - fieldName="water_primarySpeciesAggregateConcentration" + fieldName="water_primarySpeciesConstraintValue" component="6" scale="7.817e-4"/>