From 0a5f7985116e325c1078017386537599eacf6c4f Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sat, 14 Mar 2026 10:34:04 -0700 Subject: [PATCH 01/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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 56111728045a41278ee1a0d2ad862a128776a510 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 4 Sep 2026 22:03:01 -0700 Subject: [PATCH 07/12] 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 08/12] 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 9fb8b70c541329bda74b4932811c273d7f0936ca Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Mon, 14 Sep 2026 21:56:52 -0700 Subject: [PATCH 09/12] 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 22a8a6f78c441769f65f9a3086fb69baea24c2ed Mon Sep 17 00:00:00 2001 From: Omar Duran Date: Wed, 16 Sep 2026 22:27:09 -0700 Subject: [PATCH 10/12] wip: fix some typos --- .../fluid/reactivefluid/ReactiveFluidFields.hpp | 6 +++--- .../kernels/singlePhase/reactive/AccumulationKernels.hpp | 2 +- .../singlePhase/reactive/ThermalAccumulationKernels.hpp | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp index 81f2d6fb3d3..0ff28339da9 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp @@ -73,7 +73,7 @@ DECLARE_FIELD( dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentra 0, NOPLOT, WRITE_AND_READ, - "Deivatives of primarySpeciesAggregateConcentration w.r.t log primary species concentration" ); + "Derivatives of primarySpeciesAggregateConcentration w.r.t log primary species concentration" ); DECLARE_FIELD( dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, "dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations", @@ -81,7 +81,7 @@ DECLARE_FIELD( dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesCon 0, NOPLOT, WRITE_AND_READ, - "Deivatives of primarySpeciesMobileAggregateConcentration w.r.t log primary species concentration" ); + "Derivatives of primarySpeciesMobileAggregateConcentration w.r.t log primary species concentration" ); DECLARE_FIELD( secondarySpeciesConcentration, "secondarySpeciesConcentration", @@ -113,7 +113,7 @@ DECLARE_FIELD( dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, 0, NOPLOT, WRITE_AND_READ, - "Deivatives of aggregate concentration rates w.r.t log primary species concentration" ); + "Derivatives of aggregate concentration rates w.r.t log primary species concentration" ); } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index 983673c90d3..91bc14322cd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -191,7 +191,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU 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 + // Derivative of primary species amount in pore volume wrt pressure 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 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp index b22932096a8..8ab3567b8ba 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp @@ -161,7 +161,7 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe // Step 3: assemble the derivatives of the species amount balance equation w.r.t temperature for( integer is = 0; is < numSpecies; ++is ) { - // Drivative of primary species amount in pore volume wrt temperature + // Derivative 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_solventMassPerSolutionVolume /* + stack.poreVolume * m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] * From 718fc69954b98b42dd17dbc9eb4387dc00994156 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Thu, 17 Sep 2026 14:44:04 -0700 Subject: [PATCH 11/12] applied Omar's suggestions --- .../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 +- .../ReactiveSinglePhaseFluid.cpp | 25 +++----- .../ReactiveSinglePhaseFluid.hpp | 26 +++----- .../SinglePhaseReactiveTransport.cpp | 30 +++++---- .../reactive/AccumulationKernels.hpp | 23 ++++--- .../reactive/DirichletFluxComputeKernel.hpp | 25 ++++---- .../reactive/FluxComputeKernel.hpp | 62 +++++++++---------- .../reactive/SourceFluxComputeKernel.hpp | 15 ++--- .../reactive/ThermalAccumulationKernels.hpp | 13 ++-- .../ThermalDirichletFluxComputeKernel.hpp | 21 +++++-- .../reactive/ThermalFluxComputeKernel.hpp | 34 +++++----- .../ThermalSourceFluxComputeKernel.hpp | 11 ---- src/coreComponents/schema/schema.xsd | 4 +- src/coreComponents/schema/schema.xsd.other | 12 ++-- 18 files changed, 159 insertions(+), 157 deletions(-) diff --git a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml index 3982548312e..694fecc4544 100644 --- a/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml +++ b/inputFiles/singlePhaseReactiveTransport/ChainSerial_1DAllKinetic_base.xml @@ -37,6 +37,7 @@ + + solventMassFraction="0.001"/> + + solventMassFraction="0.001"/> + + solventMassFraction="0.001"/> + + solventMassFraction="0.001"/> + + solventMassFraction="0.001"/> ReactiveSinglePhaseFluid< BASE >:: ReactiveSinglePhaseFluid( string const & name, Group * const parent ): @@ -46,12 +41,12 @@ ReactiveSinglePhaseFluid( string const & name, Group * const parent ): setDescription( "Chemical System type. Available options are: " "``" + EnumStrings< ChemicalSystemType >::concat( "|" ) + "``" ); - this->registerWrapper( viewKeyStruct::solventMassPerSolutionVolumeString(), &m_solventMassPerSolutionVolume ). - setApplyDefaultValue( 1000.0 ). + this->registerWrapper( viewKeyStruct::solventMassFractionString(), &m_solventMassFraction ). + setApplyDefaultValue( 1.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." ); + setDescription( "Mass fraction of solvent in the solution [-]. Species molality [mol/kg solvent] " + "times this fraction is the amount per kg of solution, and times the fluid density " + "the amount per m^3 of solution. The default of 1 is a dilute aqueous solution." ); this->template registerField< fields::reactivefluid::initialPrimarySpeciesConcentration >( &m_initialPrimarySpeciesConcentration ); this->template registerField< fields::reactivefluid::secondarySpeciesConcentration >( &m_secondarySpeciesConcentration ); @@ -79,7 +74,7 @@ deliverClone( string const & name, Group * const parent ) const newConstitutiveRelation.m_numPrimarySpecies = m_numPrimarySpecies; newConstitutiveRelation.m_numSecondarySpecies = m_numSecondarySpecies; newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions; - newConstitutiveRelation.m_solventMassPerSolutionVolume = m_solventMassPerSolutionVolume; + newConstitutiveRelation.m_solventMassFraction = m_solventMassFraction; return clone; } @@ -128,10 +123,10 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() break; } - GEOS_THROW_IF_LE_MSG( m_solventMassPerSolutionVolume, 0.0, - GEOS_FMT( "invalid value of attribute '{}'", - viewKeyStruct::solventMassPerSolutionVolumeString() ), - InputError, this->getDataContext() ); + GEOS_THROW_IF( m_solventMassFraction <= 0.0 || m_solventMassFraction > 1.0, + GEOS_FMT( "{}: invalid value of attribute '{}', expected a value in (0, 1]", + this->getDataContext(), viewKeyStruct::solventMassFractionString() ), + InputError ); } template< typename BASE > diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index 50ed70db6f1..ee846698e6c 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -114,13 +114,13 @@ class ReactiveSinglePhaseFluid : public BASE integer numKineticReactions() const { return m_numKineticReactions; } /** - * @brief Mass of solvent per unit volume of solution [kg/m^3]. + * @brief Mass fraction of solvent in the solution [-]. * - * 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. + * HPCReact is a molality based library: concentrations, equilibrium constants and mass-action + * quotients are all on the molal scale [mol/kg solvent]. With w this fraction and rho the solution + * density, moles per kg of solution are m w and moles per m^3 of solution are m w rho. */ - real64 solventMassPerSolutionVolume() const { return m_solventMassPerSolutionVolume; } + real64 solventMassFraction() const { return m_solventMassFraction; } /** * @brief Kernel wrapper class for ReactiveSinglePhaseFluid. @@ -244,7 +244,6 @@ class ReactiveSinglePhaseFluid : public BASE { using namespace hpcReact::geochemistry; using namespace hpcReact::MoMasBenchmark; - using namespace hpcReact::bulkGeneric; using namespace hpcReact::ChainGeneric; switch( m_chemicalSystemType ) { @@ -339,7 +338,7 @@ class ReactiveSinglePhaseFluid : public BASE struct viewKeyStruct : ConstitutiveBase::viewKeyStruct { static constexpr char const * chemicalSystemNameString() { return "chemicalSystemType"; } - static constexpr char const * solventMassPerSolutionVolumeString() { return "solventMassPerSolutionVolume"; } + static constexpr char const * solventMassFractionString() { return "solventMassFraction"; } }; protected: @@ -376,16 +375,9 @@ class ReactiveSinglePhaseFluid : public BASE ChemicalSystemType m_chemicalSystemType; - /// 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; + /// Mass fraction of solvent in the solution [-], constant; for the carbonate brine EQ3/6 gives 0.898 + /// TODO: depends on concentration, w = 1 / (1 + sum_i m_i M_i) + real64 m_solventMassFraction; }; // these aliases are useful in constitutive dispatch diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 97732290591..0d2632b065a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -472,9 +472,9 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, } ConstitutiveManager const & cm = domain.getConstitutiveManager(); - real64 const solventMassPerSolutionVolume = - m_isThermal ? cm.getConstitutiveRelation< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventMassPerSolutionVolume() - : cm.getConstitutiveRelation< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventMassPerSolutionVolume(); + real64 const solventMassFraction = + m_isThermal ? cm.getConstitutiveRelation< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventMassFraction() + : cm.getConstitutiveRelation< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( m_reactiveFluidModelName ).solventMassFraction(); 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(), - solventMassPerSolutionVolume, + solventMassFraction, dofManager.rankOffset(), dofKey, getName(), @@ -512,7 +512,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies, m_hasDiffusion, mobilePrimarySpeciesFlags.toViewConst(), - solventMassPerSolutionVolume, + solventMassFraction, dofManager.rankOffset(), dofKey, getName(), @@ -569,16 +569,20 @@ 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 solventMassPerSolutionVolume = fluid.solventMassPerSolutionVolume(); + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const density = fluid.density(); + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const density_n = fluid.density_n(); + real64 const solventMassFraction = fluid.solventMassFraction(); + // moles in the pore volume: molality * solvent mass fraction * density * pore volume 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] * solventMassPerSolutionVolume; + primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventMassFraction * density[ei][0]; + // the previous amount is not written to restart files and is rebuilt here after a restart if( isZero( primarySpeciesAggregateMole_n[ei][is] ) ) - primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventMassPerSolutionVolume; + primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventMassFraction * density_n[ei][0]; } } ); } @@ -588,16 +592,20 @@ 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 solventMassPerSolutionVolume = fluid.solventMassPerSolutionVolume(); + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const density = fluid.density(); + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const density_n = fluid.density_n(); + real64 const solventMassFraction = fluid.solventMassFraction(); + // moles in the pore volume: molality * solvent mass fraction * density * pore volume 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] * solventMassPerSolutionVolume; + primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventMassFraction * density[ei][0]; + // the previous amount is not written to restart files and is rebuilt here after a restart if( isZero( primarySpeciesAggregateMole_n[ei][is] ) ) - primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventMassPerSolutionVolume; + primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventMassFraction * density_n[ei][0]; } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index 91bc14322cd..117ea50f94a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -82,7 +82,9 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU arrayView1d< real64 > const & localRhs ) : Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ), m_dt( dt ), - m_solventMassPerSolutionVolume( fluid.solventMassPerSolutionVolume() ), + m_solventMassFraction( fluid.solventMassFraction() ), + m_density( fluid.density() ), + m_dDensity( fluid.dDensity() ), m_volume( subRegion.getElementVolume() ), m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ), m_porosity( solid.getPorosity() ), @@ -161,7 +163,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU void computeAccumulation( localIndex const ei, StackVariables & stack ) const { - // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventMassPerSolutionVolume * stack.poreVolume - + // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventMassFraction * density * stack.poreVolume - // primarySpeciesAggregateMole_n[is]) // - dt * m_volume * primarySpeciesKineticRate[is] @@ -185,14 +187,17 @@ 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_solventMassPerSolutionVolume * stack.poreVolume; + // moles of species per m^3 of solution: molality * solvent mass fraction * density + real64 const aggregateConcMolarity = m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassFraction * m_density[ei][0]; + stack.localResidual[is+numEqn-numSpecies] += aggregateConcMolarity * 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 // Derivative of primary species amount in pore volume wrt pressure - stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassPerSolutionVolume + stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * aggregateConcMolarity + + stack.poreVolume * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassFraction * m_dDensity[ei][0][DerivOffset::dP] /* + 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] ) * @@ -204,7 +209,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_solventMassPerSolutionVolume; + m_solventMassFraction * m_density[ei][0]; stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js]; } @@ -242,8 +247,12 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU /// Time step size real64 const m_dt; - /// Mass of solvent per unit volume of solution [kg/m³], converting molality [mol/kg] to molarity [mol/m³] - real64 const m_solventMassPerSolutionVolume; + /// Mass fraction of solvent in the solution [-]; molality times this fraction times density is the molarity + real64 const m_solventMassFraction; + + /// View on the fluid density and its derivatives + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_density; + arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > const m_dDensity; /// 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 2ef73fddea2..97eae3adeda 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 & solventMassPerSolutionVolume, + real64 const & solventMassFraction, 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_solventMassPerSolutionVolume( solventMassPerSolutionVolume ) + m_solventMassFraction( solventMassFraction ) {} /** @@ -188,25 +188,24 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo real64 const & dens_up, real64 const & dDens_dP_up ) { - GEOS_UNUSED_VAR( kf ); + GEOS_UNUSED_VAR( kf, dens_up, dDens_dP_up ); real64 speciesFlux[numSpecies]{}; real64 dSpeciesFlux_dP[numSpecies]{}; real64 dSpeciesFlux_dLogConc[numSpecies][numSpecies]{}; + // the mass flux carries molality * solvent mass fraction moles per kg of solution for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is] * m_solventMassPerSolutionVolume; - speciesFlux[is] = aggregateConcMolarity_i / dens_up * fluxVal * mobility_up; + real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is] * m_solventMassFraction; + speciesFlux[is] = aggregateConcPerMass_i * fluxVal * mobility_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); + dSpeciesFlux_dP[is] = aggregateConcPerMass_i * ( dFlux_dP * mobility_up + fluxVal * dMobility_dP_up ); for( integer js = 0; js < numSpecies; ++js ) { - 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; + real64 const dAggregateConcPerMass_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[seri][sesri][sei][0][is][js] * m_solventMassFraction; + dSpeciesFlux_dLogConc[is][js] += dAggregateConcPerMass_i_dLogConc_j * fluxVal * mobility_up; } } @@ -269,8 +268,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; - /// Mass of solvent per unit volume of solution [kg/m³], converting molality [mol/kg] to molarity [mol/m³] - real64 const m_solventMassPerSolutionVolume; + /// Mass fraction of solvent in the solution [-]; molality times this fraction is the amount per kg of solution + real64 const m_solventMassFraction; }; @@ -349,7 +348,7 @@ class DirichletFluxComputeKernelFactory reactiveFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, - reactiveFluid.solventMassPerSolutionVolume(), + reactiveFluid.solventMassFraction(), 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 dcc870ff8c5..c2ad0d8a66e 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] solventMassPerSolutionVolume mass of solvent per unit volume of solution [kg/m3] + * @param[in] solventMassFraction mass fraction of solvent in the solution [-] * @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 & solventMassPerSolutionVolume, + real64 const & solventMassFraction, 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_solventMassPerSolutionVolume( solventMassPerSolutionVolume ) + m_solventMassFraction( solventMassFraction ) {} /** @@ -211,7 +211,6 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E StackVariables & stack, FUNC && kernelOp = NoOpFunc{} ) const { - using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; // *********************************************** // First, we call the base computeFlux to compute: // 1) massFlux and its derivatives, @@ -246,28 +245,22 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E localIndex const esr_up = sesri[k_up]; localIndex const ei_up = sei[k_up]; - real64 const fluidDens_up = m_dens[er_up][esr_up][ei_up][0]; - real64 const dDens_dPres = m_dDens[er_up][esr_up][ei_up][0][DerivOffset::dP]; - - // compute species fluxes and derivatives using upstream cell concentration + // compute species fluxes and derivatives using upstream cell concentration: the mass flux carries + // molality * solvent mass fraction moles per kg of solution, so the density does not enter for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] - * m_solventMassPerSolutionVolume; // convert from mol/kg to mol/m3 - speciesFlux[is] = aggregateConcMolarity_i / fluidDens_up * fluxVal; + real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] * m_solventMassFraction; + speciesFlux[is] = aggregateConcPerMass_i * fluxVal; for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { - dSpeciesFlux_dP[ke][is] += aggregateConcMolarity_i / fluidDens_up * dFlux_dP[ke]; + dSpeciesFlux_dP[ke][is] += aggregateConcPerMass_i * dFlux_dP[ke]; } - dSpeciesFlux_dP[k_up][is] += -aggregateConcMolarity_i * fluxVal * dDens_dPres / (fluidDens_up * fluidDens_up); - 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_solventMassPerSolutionVolume; // convert from mol/kg to mol/m3 - dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConcMolarity_i_dLogConc_j / fluidDens_up * fluxVal; + real64 const dAggregateConcPerMass_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js] * m_solventMassFraction; + dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConcPerMass_i_dLogConc_j * fluxVal; } } @@ -296,7 +289,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E } // Customize the kernel with this lambda - kernelOp( k, seri, sesri, sei, connectionIndex, alpha, mobility, potGrad, fluxVal, dFlux_dP, fluidDens_up ); + kernelOp( k, seri, sesri, sei, connectionIndex, alpha, mobility, potGrad, fluxVal, dFlux_dP ); } ); } @@ -313,6 +306,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E StackVariables & stack, FUNC && kernelOp = NoOpFunc{} ) const { + using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; if( m_hasDiffusion ) { // ***************************************************** @@ -339,7 +333,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E // clear working arrays real64 diffusionFlux[numSpecies]{}; real64 speciesGrad[numSpecies]{}; - // real64 dDiffusionFlux_dP[numFluxSupportPoints][numSpecies]{}; // Turn on if diffusionFlux is pressure-dependent + real64 dDiffusionFlux_dP[numFluxSupportPoints][numSpecies]{}; real64 dDiffusionFlux_dLogConc[numFluxSupportPoints][numSpecies][numSpecies]{}; real64 const diffusionTrans[numFluxSupportPoints] = { stack.diffusionTransmissibility[connectionIndex][0], @@ -349,26 +343,27 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E // loop over primary species for( integer is = 0; is < numSpecies; ++is ) { - // real64 dSpeciesGrad_i_dP[numFluxSupportPoints]{}; // Turn on if speciesGrad is pressure-dependent + real64 dSpeciesGrad_i_dP[numFluxSupportPoints]{}; real64 dSpeciesGrad_i_dLogConc[numFluxSupportPoints][numSpecies]{}; - // Step 2: compute species gradient + // Step 2: compute the gradient of the molarity, molality * solvent mass fraction * density for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { localIndex const er = seri[ke]; localIndex const esr = sesri[ke]; localIndex const ei = sei[ke]; - real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] - * m_solventMassPerSolutionVolume; // convert from mol/kg to mol/m3 + real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventMassFraction; + real64 const dens = m_dens[er][esr][ei][0]; - speciesGrad[is] += diffusionTrans[ke] * aggregateConcMolarity_i; + speciesGrad[is] += diffusionTrans[ke] * aggregateConcPerMass_i * dens; + dSpeciesGrad_i_dP[ke] += diffusionTrans[ke] * aggregateConcPerMass_i * m_dDens[er][esr][ei][0][DerivOffset::dP]; for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js] * m_solventMassPerSolutionVolume; + real64 const dAggregateConcPerMass_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js] * m_solventMassFraction; - dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConcMolarity_i_dLogConc_j; + dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConcPerMass_i_dLogConc_j * dens; } } @@ -385,6 +380,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E // add contributions of the derivatives of component fractions wrt pressure/component fractions for( integer ke = 0; ke < numFluxSupportPoints; ke++ ) { + dDiffusionFlux_dP[ke][is] += m_referencePorosity[er_up][esr_up][ei_up] * dSpeciesGrad_i_dP[ke]; for( integer js = 0; js < numSpecies; ++js ) { dDiffusionFlux_dLogConc[ke][is][js] += m_referencePorosity[er_up][esr_up][ei_up] * dSpeciesGrad_i_dLogConc[ke][js]; @@ -401,8 +397,8 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { localIndex const localDofIndexPres = k[ke] * numDof; - // stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dDiffusionFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is]; - // stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dDiffusionFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is]; + stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dDiffusionFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is]; + stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dDiffusionFlux_dP[ke][is] * m_mobilePrimarySpeciesFlags[is]; for( integer js = 0; js < numSpecies; ++js ) { @@ -507,8 +503,8 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; - /// Mass of solvent per unit volume of solution [kg/m3], converting molality [mol/kg] to molarity [mol/m3] - real64 const m_solventMassPerSolutionVolume; + /// Mass fraction of solvent in the solution [-]; molality times this fraction is the amount per kg of solution + real64 const m_solventMassFraction; }; /** @@ -525,7 +521,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] solventMassPerSolutionVolume mass of solvent per unit volume of solution [kg/m3] + * @param[in] solventMassFraction mass fraction of solvent in the solution [-] * @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 +536,7 @@ class FluxComputeKernelFactory createAndLaunch( integer const numSpecies, integer const hasDiffusion, arrayView1d< integer const > const mobilePrimarySpeciesFlags, - real64 const solventMassPerSolutionVolume, + real64 const solventMassFraction, globalIndex const rankOffset, string const & dofKey, string const & solverName, @@ -572,7 +568,7 @@ class FluxComputeKernelFactory KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor, flowAccessors, reactiveFlowAccessors, fluidAccessors, reactiveFluidAccessors, permAccessors, diffusionAccessors, porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags, - solventMassPerSolutionVolume, dt, localMatrix, localRhs ); + solventMassFraction, 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 a113f1771a9..94a41eeea30 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_solventMassPerSolutionVolume( fluid.solventMassPerSolutionVolume() ), + m_solventMassFraction( fluid.solventMassFraction() ), m_primarySpeciesAggregateConcentration( fluid.primarySpeciesAggregateConcentration() ), m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations( fluid.dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations() ), m_density( fluid.density() ), @@ -150,17 +150,14 @@ class SourceFluxComputeKernel { real64 const scaledInflowMass = stack.totalInflowMass / m_sizeScalingFactor; + // the inflow mass carries molality * solvent mass fraction moles per kg of solution 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]) * - scaledInflowMass; + stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventMassFraction * 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_solventMassFraction * scaledInflowMass; } } @@ -242,8 +239,8 @@ class SourceFluxComputeKernel /// size scaling factor real64 const m_sizeScalingFactor; - /// Mass of solvent per unit volume of solution [kg/m³], converting molality [mol/kg] to molarity [mol/m³] - real64 const m_solventMassPerSolutionVolume; + /// Mass fraction of solvent in the solution [-]; molality times this fraction is the amount per kg of solution + real64 const m_solventMassFraction; // 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 8ab3567b8ba..92ecf6d6c20 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp @@ -53,7 +53,9 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe using Base::m_volume; using Base::m_deltaVolume; using Base::m_primarySpeciesAggregateConcentration; - using Base::m_solventMassPerSolutionVolume; + using Base::m_solventMassFraction; + using Base::m_density; + using Base::m_dDensity; /// Note: Derivative lineup only supports dP & dT, not component terms using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; @@ -162,10 +164,11 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe for( integer is = 0; is < numSpecies; ++is ) { // Derivative 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_solventMassPerSolutionVolume - /* + stack.poreVolume * - m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] * - m_solventMassPerSolutionVolume */; + real64 const aggregateConcMolality = m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventMassFraction; + stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += ( stack.dPoreVolume_dTemp * m_density[ei][0] + + stack.poreVolume * m_dDensity[ei][0][DerivOffset::dT] ) * aggregateConcMolality + /* + stack.poreVolume * m_density[ei][0] * m_solventMassFraction * + m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] */; // // 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 740277c402f..72e43913d9c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp @@ -86,7 +86,9 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl using Base::m_sei; using Base::m_facePres; using Base::m_faceGravCoef; - using Base::m_solventMassPerSolutionVolume; + using Base::m_solventMassFraction; + using Base::m_primarySpeciesMobileAggregateConc; + using Base::m_mobilePrimarySpeciesFlags; using ReactiveSinglePhaseFlowAccessors = typename Base::ReactiveSinglePhaseFlowAccessors; using ReactiveSinglePhaseFluidAccessors = typename Base::ReactiveSinglePhaseFluidAccessors; @@ -138,7 +140,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl PermeabilityAccessors const & permeabilityAccessors, ThermalConductivityAccessors const & thermalConductivityAccessors, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, - real64 const & solventMassPerSolutionVolume, + real64 const & solventMassFraction, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -154,7 +156,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl reactiveSinglePhaseFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, - solventMassPerSolutionVolume, + solventMassFraction, dt, localMatrix, localRhs ), @@ -233,10 +235,11 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl // Compute the derivatives of the (upwinded) energy flux wrt pressure and temperature + real64 dFlux_dT = mobility_up * dF_dT; if( f >= 0 ) // the element is upstream { real64 const dFlux_dP = mobility_up * dF_dP + dMobility_dP_up * f; - real64 const dFlux_dT = mobility_up * dF_dT + m_dMob[er][esr][ei][DerivOffset::dT] * f; + dFlux_dT += m_dMob[er][esr][ei][DerivOffset::dT] * f; stack.dEnergyFlux_dP += dFlux_dP * enthalpy + flux * m_dEnthalpy[er][esr][ei][0][DerivOffset::dP]; stack.dEnergyFlux_dT += dFlux_dT * enthalpy + flux * m_dEnthalpy[er][esr][ei][0][DerivOffset::dT]; @@ -244,12 +247,18 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl else { real64 const dFlux_dP = mobility_up * dF_dP; - real64 const dFlux_dT = mobility_up * dF_dT; stack.dEnergyFlux_dP += dFlux_dP * enthalpy; stack.dEnergyFlux_dT += dFlux_dT * enthalpy; } + // derivative of the species flux, molality * solvent mass fraction * mass flux, wrt temperature + for( integer is = 0; is < numSpecies; ++is ) + { + real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventMassFraction; + stack.localFluxJacobian[numEqn - numSpecies + is][numDof - numSpecies - 1] += m_dt * aggregateConcPerMass_i * dFlux_dT * m_mobilePrimarySpeciesFlags[is]; + } + // Contribution of energy conduction through the solid phase real64 thermalTrans = 0.0; real64 dThermalTrans_dThermalCond[3]{}; @@ -395,7 +404,7 @@ class DirichletFluxComputeKernelFactory permeabilityAccessors, thermalConductivityAccessors, mobilePrimarySpeciesFlags, - reactiveFluid.solventMassPerSolutionVolume(), + reactiveFluid.solventMassFraction(), 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 6ed2df03064..cf07b134b11 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_solventMassPerSolutionVolume; + using Base::m_solventMassFraction; 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] solventMassPerSolutionVolume mass of solvent per unit volume of solution [kg/m3] + * @param[in] solventMassFraction mass fraction of solvent in the solution [-] * @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 & solventMassPerSolutionVolume, + real64 const & solventMassFraction, 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, - solventMassPerSolutionVolume, + solventMassFraction, dt, localMatrix, localRhs ), @@ -237,8 +237,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne real64 const mobility, real64 const & potGrad, real64 const & fluxVal, - real64 const (&dFlux_dP)[2], - real64 const fluidDens_up ) + real64 const (&dFlux_dP)[2] ) { // Step 1: compute the derivatives of the (upwinded) massFlux wrt temperature // -------------------------------------------------------------------------- @@ -319,28 +318,25 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne real64 dSpeciesFlux_dT[numFluxSupportPoints][numSpecies]{}; { - // Step 2.1: compute the derivatives of the upstream density wrt temperature - // choose upstream cell + // Step 2.1: choose upstream cell; the mass flux carries molality * solvent mass fraction moles per kg + // of solution, so the density does not enter and only the mass flux depends on temperature localIndex const k_up = (potGrad >= 0) ? 0 : 1; localIndex const er_up = seri[k_up]; localIndex const esr_up = sesri[k_up]; localIndex const ei_up = sei[k_up]; - real64 const dDens_dTemp = m_dDens[er_up][esr_up][ei_up][0][DerivOffset::dT]; - // 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_solventMassPerSolutionVolume; + real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] * m_solventMassFraction; // real64 const dAggregateConc_i_dTemp = m_dPrimarySpeciesMobileAggregateConcentration_dTemp[er_up][esr_up][ei_up][is]; - // 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); + // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * m_solventMassFraction * fluxVal; for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { - dSpeciesFlux_dT[ke][is] += aggregateConcMolarity_i / fluidDens_up * dFlux_dT[ke]; + dSpeciesFlux_dT[ke][is] += aggregateConcPerMass_i * dFlux_dT[ke]; } } } @@ -472,6 +468,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne void computeDiffusion( localIndex const iconn, StackVariables & stack ) const { + using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; Base::computeDiffusion( iconn, stack, [&] ( integer const is, localIndex const (&k)[2], localIndex const (&seri)[2], @@ -493,7 +490,10 @@ 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_solventMassPerSolutionVolume; + // molarity is molality * solvent mass fraction * density, so both the transmissibility and the density vary with temperature + real64 const aggregateConcPerMass_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventMassFraction; + dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * aggregateConcPerMass_i * m_dens[er][esr][ei][0] + + stack.diffusionTransmissibility[connectionIndex][ke] * aggregateConcPerMass_i * m_dDens[er][esr][ei][0][DerivOffset::dT]; } 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 solventMassPerSolutionVolume, + real64 const solventMassFraction, 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, solventMassPerSolutionVolume, dt, localMatrix, localRhs ); + hasDiffusion, mobilePrimarySpeciesFlags, solventMassFraction, 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 739b6936085..4184f7de27b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp @@ -45,10 +45,6 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu using Base::numDof; using Base::numEqn; using Base::m_sizeScalingFactor; - using Base::m_solventMassPerSolutionVolume; - using Base::m_primarySpeciesAggregateConcentration; - using Base::m_density; - using Base::m_dDensity; using Base::m_localMatrix; using Base::m_localRhs; @@ -120,13 +116,6 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu stack.localEnergyRhs += m_enthalpy[ei][0] * scaledInflowMass; stack.localEnergyJacobian[0] = scaledInflowMass * m_dEnthalpy[ei][0][DerivOffset::dP]; stack.localEnergyJacobian[numDof-numSpecies-1] = scaledInflowMass * m_dEnthalpy[ei][0][DerivOffset::dT]; - - for( integer i = 0; i < numSpecies; ++i ) - { - 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 5fcbfa97ae3..c719a937ad1 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -8924,8 +8924,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia - - + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index cf0ac4ad56b..7082b1875de 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -3153,7 +3153,7 @@ A field can represent a physical variable. (pressure, temperature, global compos - + @@ -3161,9 +3161,9 @@ A field can represent a physical variable. (pressure, temperature, global compos - + - + @@ -3218,7 +3218,7 @@ A field can represent a physical variable. (pressure, temperature, global compos - + @@ -3226,9 +3226,9 @@ A field can represent a physical variable. (pressure, temperature, global compos - + - + From 83fe1cd18f46ea823105f0572d8eeda5ed589384 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Sun, 20 Sep 2026 09:54:17 -0700 Subject: [PATCH 12/12] rebaseline --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index b6441b60bc5..74cea56702e 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr4108-17561-eb62184 + baseline: integratedTests/baseline_integratedTests-pr3996-17639-718fc69 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 1fbab508e88..e700a96f185 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -4,6 +4,11 @@ Notes This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3996 (2026-09-20) +===================== +Switch the single-phase reactive transport solver to a molality-based primary +species concentration by introducing solventMassFraction. + 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.