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/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/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..e3e3dee58c8 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() ), @@ -160,8 +161,8 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU void computeAccumulation( localIndex const ei, StackVariables & stack ) const { - // Residual[is] += (primarySpeciesAggregateConcentration[is] * stack.poreVolume - primarySpeciesAggregateMole_n[is]) - // - dt * m_volume * primarySpeciesKineticRate[is] // To Check: what's the unit of the kinetic rate + // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventDensity * stack.poreVolume - primarySpeciesAggregateMole_n[is]) + // - dt * m_volume * primarySpeciesKineticRate[is] * solventDensity Base::computeAccumulation( ei, stack ); @@ -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; } }