diff --git a/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml b/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml index 3c032ca0754..03d19fbbb7f 100644 --- a/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml +++ b/inputFiles/singlePhaseReactiveTransport/MoMaSEasy_1DDiffusive_base.xml @@ -19,7 +19,7 @@ lineSearchAction="None" timeStepIncreaseFactor="1.2" logLevel="2" - maxAllowedResidualNorm="1e100"/> + maxAllowedResidualNorm="1e8"/> diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_base.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_base.xml new file mode 100644 index 00000000000..625440f19c2 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_base.xml @@ -0,0 +1,155 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_smoke.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_smoke.xml new file mode 100644 index 00000000000..67eede6b0f9 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_1D_smoke.xml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_base.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_base.xml new file mode 100644 index 00000000000..2aca6139a30 --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_base.xml @@ -0,0 +1,133 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_smoke.xml b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_smoke.xml new file mode 100644 index 00000000000..9db70cb152d --- /dev/null +++ b/inputFiles/singlePhaseReactiveTransport/mixedReactionCarbonateSystemBdot_2D_smoke.xml @@ -0,0 +1,49 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats b/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats index 88bff398cfe..591e341a851 100644 --- a/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats +++ b/inputFiles/singlePhaseReactiveTransport/singlePhaseReactiveTransport.ats @@ -40,6 +40,23 @@ decks = [ restart_step=10, check_step=20, restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="mixedReactionCarbonateSystemBdot_1D_smoke", + description= + '1D Carbonate mixed reaction with the B-dot activity model (1D, mixedReaction, bdot)', + partitions=((1, 1, 1), (2, 1, 1)), + restart_step=10, + check_step=20, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + + TestDeck( + name="mixedReactionCarbonateSystemBdot_2D_smoke", + description= + '2D Carbonate mixed reaction with the B-dot activity model (2D, mixedReaction, bdot)', + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=20, + check_step=40, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), ] generate_geos_tests(decks) diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 39ddb39771b..90cd5405359 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -363,6 +363,7 @@ if (ENABLE_HPCREACT) set( constitutive_headers ${constitutive_headers} fluid/reactivefluid/ReactiveFluidSelector.hpp + fluid/reactivefluid/ReactiveFluidSystemSelector.hpp ) set( constitutive_headers diff --git a/src/coreComponents/constitutive/HPCReact b/src/coreComponents/constitutive/HPCReact index cee75f0309b..ebd67f2f4fc 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit cee75f0309b113845c5f7c48a46f4d14f376c844 +Subproject commit ebd67f2f4fcb4efe6a1fbb48a8fa83022488aaa0 diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp new file mode 100644 index 00000000000..a70e31407ed --- /dev/null +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp @@ -0,0 +1,267 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ReactiveFluidSystemSelector.hpp + * @brief The chemical systems ReactiveSinglePhaseFluid can instantiate, and the activity model each + * may be paired with. + * + * Adding a system, or an activity model for one, is one entry below plus its name in + * ReactionSystemList. The kernel wrapper variant, the dispatch, and the check that rejects an + * unsupported pairing all follow from that list. + */ +#ifndef GEOS_CONSTITUTIVE_FLUID_REACTIVEFLUID_REACTIVEFLUIDSYSTEMSELECTOR_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_REACTIVEFLUID_REACTIVEFLUIDSYSTEMSELECTOR_HPP_ + +#include "common/DataTypes.hpp" +#include "common/format/EnumStrings.hpp" + +#include "constitutive/HPCReact/src/reactions/geochemistry/GeochemicalSystems.hpp" +#include "constitutive/HPCReact/src/reactions/exampleSystems/BulkGeneric.hpp" +#include "constitutive/HPCReact/src/reactions/exampleSystems/ChainGeneric.hpp" +#include "constitutive/HPCReact/src/reactions/exampleSystems/MoMasBenchmark.hpp" + +#include +#include + +namespace geos +{ + +namespace constitutive +{ + +namespace reactivefluid +{ + +/** + * @brief The chemical system a ReactiveSinglePhaseFluid solves. + */ +enum class ChemicalSystemType : integer +{ + carbonate, + carbonateAllEquilibrium, + ultramafic, + momasEasy, + momasMedium, + chainSerialAllKinetic +}; + +/** + * @brief The activity model applied to it. + * @details ``identity`` is the identity activity model, gamma = 1 and a_w = 1. ``bdot`` needs ion + * size and b-dot parameters, which only the geochemical systems carry. + */ +enum class ActivityModelType : integer +{ + identity, + bdot +}; + +ENUM_STRINGS( ChemicalSystemType, + "carbonate", + "carbonateAllEquilibrium", + "ultramafic", + "momasEasy", + "momasMedium", + "chainSerialAllKinetic" ); + +ENUM_STRINGS( ActivityModelType, + "identity", + "bdot" ); + +/** + * @brief What every entry below shares: the two HPCReact types, and the two enum values the XML + * input selects it by. + * @tparam SYSTEM the ChemicalSystemType it answers to + * @tparam MODEL the ActivityModelType it answers to + * @tparam REACTION_PARAMS the HPCReact reaction parameters type + * @tparam ACTIVITY_MODEL the HPCReact activity model type + * + * Each entry adds the two parameter objects, as accessors returning by value. They cannot be + * template arguments: the HPCReact objects are namespace-scope constexpr and so have internal + * linkage, which would make the entry a different type in every translation unit. Returning a + * copy also matches how the kernel wrapper stores them. + */ +template< ChemicalSystemType SYSTEM, + ActivityModelType MODEL, + typename REACTION_PARAMS, + typename ACTIVITY_MODEL > +struct ReactionSystemEntry +{ + /// the HPCReact reaction parameters type + using ReactionParamsType = REACTION_PARAMS; + + /// the HPCReact activity model type + using ActivityType = ACTIVITY_MODEL; + + /// the chemical system this answers to + static constexpr ChemicalSystemType chemicalSystem = SYSTEM; + + /// the activity model this answers to + static constexpr ActivityModelType activityModel = MODEL; +}; + +/// Carbonate, 16 species with calcite a pure phase, identity activity model. +struct CarbonateIdentityEntry : ReactionSystemEntry< ChemicalSystemType::carbonate, + ActivityModelType::identity, + hpcReact::geochemistry::carbonateSystemType, + hpcReact::geochemistry::carbonateNosolidIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateNosolidIdentityActivityParams; } +}; + +/// Carbonate, 16 species with calcite a pure phase, B-dot activity model on EQ3/6 parameters. +struct CarbonateBdotEntry : ReactionSystemEntry< ChemicalSystemType::carbonate, + ActivityModelType::bdot, + hpcReact::geochemistry::carbonateSystemType, + hpcReact::geochemistry::carbonateNosolidActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateNosolidActivityParamsEQ36; } +}; + +/// Carbonate, 17 species with calcite among them, identity activity model. +struct CarbonateAllEquilibriumIdentityEntry : ReactionSystemEntry< ChemicalSystemType::carbonateAllEquilibrium, + ActivityModelType::identity, + hpcReact::geochemistry::carbonateSystemAllEquilibriumType, + hpcReact::geochemistry::carbonateIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystemAllEquilibrium; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateIdentityActivityParams; } +}; + +/// Carbonate, 17 species with calcite among them, B-dot activity model on EQ3/6 parameters. +struct CarbonateAllEquilibriumBdotEntry : ReactionSystemEntry< ChemicalSystemType::carbonateAllEquilibrium, + ActivityModelType::bdot, + hpcReact::geochemistry::carbonateSystemAllEquilibriumType, + hpcReact::geochemistry::carbonateActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::carbonateSystemAllEquilibrium; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::carbonateActivityParamsEQ36; } +}; + +/// Ultramafic, 20 species with the five minerals pure phases, identity activity model. +struct UltramaficIdentityEntry : ReactionSystemEntry< ChemicalSystemType::ultramafic, + ActivityModelType::identity, + hpcReact::geochemistry::ultramaficSystemType, + hpcReact::geochemistry::ultramaficIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::ultramaficSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::ultramaficIdentityActivityParams; } +}; + +/// Ultramafic, 20 species with the five minerals pure phases, B-dot activity model on EQ3/6 parameters. +struct UltramaficBdotEntry : ReactionSystemEntry< ChemicalSystemType::ultramafic, + ActivityModelType::bdot, + hpcReact::geochemistry::ultramaficSystemType, + hpcReact::geochemistry::ultramaficActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::geochemistry::ultramaficSystem; } + static ActivityType::Params activityParams() { return hpcReact::geochemistry::ultramaficActivityParamsEQ36; } +}; + +/// MoMaS easy benchmark, 12 abstract species, identity activity model. +struct MomasEasyIdentityEntry : ReactionSystemEntry< ChemicalSystemType::momasEasy, + ActivityModelType::identity, + hpcReact::MoMasBenchmark::easyCaseType, + hpcReact::MoMasBenchmark::easyCaseIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::MoMasBenchmark::easyCaseParams; } + static ActivityType::Params activityParams() { return hpcReact::MoMasBenchmark::easyCaseIdentityActivityParams; } +}; + +/// MoMaS medium benchmark, 14 abstract species, identity activity model. +struct MomasMediumIdentityEntry : ReactionSystemEntry< ChemicalSystemType::momasMedium, + ActivityModelType::identity, + hpcReact::MoMasBenchmark::mediumCaseType, + hpcReact::MoMasBenchmark::mediumCaseIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::MoMasBenchmark::mediumCaseParams; } + static ActivityType::Params activityParams() { return hpcReact::MoMasBenchmark::mediumCaseIdentityActivityParams; } +}; + +/// Generic serial chain, 3 species all kinetic, identity activity model. +struct ChainSerialIdentityEntry : ReactionSystemEntry< ChemicalSystemType::chainSerialAllKinetic, + ActivityModelType::identity, + hpcReact::ChainGeneric::serialAllKineticType, + hpcReact::ChainGeneric::serialAllKineticIdentityActivityType > +{ + static ReactionParamsType reactionParams() { return hpcReact::ChainGeneric::serialAllKineticParams; } + static ActivityType::Params activityParams() { return hpcReact::ChainGeneric::serialAllKineticIdentityActivityParams; } +}; + +/** + * @brief Every pairing the fluid model can build. B-dot appears only for the geochemical systems, + * the ones carrying ion size and b-dot parameters. + */ +using ReactionSystemList = std::variant< CarbonateIdentityEntry, + CarbonateBdotEntry, + CarbonateAllEquilibriumIdentityEntry, + CarbonateAllEquilibriumBdotEntry, + UltramaficIdentityEntry, + UltramaficBdotEntry, + MomasEasyIdentityEntry, + MomasMediumIdentityEntry, + ChainSerialIdentityEntry >; + +/** + * @brief Apply @p func to one instance of every system in ReactionSystemList. + * @tparam FUNC the callable type + * @tparam INDICES the alternatives of ReactionSystemList + * @param func a generic callable taking one entry + */ +template< typename FUNC, std::size_t ... INDICES > +void forEachReactionSystem( FUNC && func, std::index_sequence< INDICES... > ) +{ + ( func( std::variant_alternative_t< INDICES, ReactionSystemList >{} ), ... ); +} + +/** + * @copydoc forEachReactionSystem + */ +template< typename FUNC > +void forEachReactionSystem( FUNC && func ) +{ + forEachReactionSystem( std::forward< FUNC >( func ), + std::make_index_sequence< std::variant_size_v< ReactionSystemList > >{} ); +} + +/** + * @brief Whether @p chemicalSystem may be paired with @p activityModel. + * @param chemicalSystem the chemical system + * @param activityModel the activity model + * @return true if ReactionSystemList holds that pairing + */ +inline bool isSupportedReactionSystem( ChemicalSystemType const chemicalSystem, + ActivityModelType const activityModel ) +{ + bool supported = false; + forEachReactionSystem( [&]( auto system ) + { + using System = decltype( system ); + supported = supported || ( System::chemicalSystem == chemicalSystem && + System::activityModel == activityModel ); + } ); + return supported; +} + +} // namespace reactivefluid + +} // namespace constitutive + +} // namespace geos + +#endif // GEOS_CONSTITUTIVE_FLUID_REACTIVEFLUID_REACTIVEFLUIDSYSTEMSELECTOR_HPP_ diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index b73ffae31ae..b0a159297f3 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -41,6 +41,14 @@ ReactiveSinglePhaseFluid( string const & name, Group * const parent ): setDescription( "Chemical System type. Available options are: " "``" + EnumStrings< ChemicalSystemType >::concat( "|" ) + "``" ); + this->registerWrapper( viewKeyStruct::activityModelNameString(), &m_activityModelType ). + setApplyDefaultValue( ActivityModelType::identity ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Activity model applied to the chemical system. Available options are: " + "``" + EnumStrings< ActivityModelType >::concat( "|" ) + + ". ``bdot`` requires ion size and b-dot parameters, which only the geochemical " + "systems carry." ); + this->registerWrapper( viewKeyStruct::solventMassFractionString(), &m_solventMassFraction ). setApplyDefaultValue( 1.0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -71,6 +79,7 @@ deliverClone( string const & name, Group * const parent ) const ReactiveSinglePhaseFluid & newConstitutiveRelation = dynamicCast< ReactiveSinglePhaseFluid & >( *clone ); newConstitutiveRelation.m_chemicalSystemType = m_chemicalSystemType; + newConstitutiveRelation.m_activityModelType = m_activityModelType; newConstitutiveRelation.m_numPrimarySpecies = m_numPrimarySpecies; newConstitutiveRelation.m_numSecondarySpecies = m_numSecondarySpecies; newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions; @@ -84,44 +93,25 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() { BASE::postInputInitialization(); - switch( m_chemicalSystemType ) + GEOS_THROW_IF( !isSupportedReactionSystem( m_chemicalSystemType, m_activityModelType ), + GEOS_FMT( "{}: chemical system '{}' has no '{}' activity model.", + this->getDataContext(), + EnumStrings< ChemicalSystemType >::toString( m_chemicalSystemType ), + EnumStrings< ActivityModelType >::toString( m_activityModelType ) ), + InputError ); + + // taken from the reaction system itself + forEachReactionSystem( [&]( auto system ) { - case ChemicalSystemType::ultramafic: - m_numPrimarySpecies = 9; - m_numSecondarySpecies = 16; - m_numKineticReactions = 5; - break; - - case ChemicalSystemType::carbonate: - m_numPrimarySpecies = 7; - m_numSecondarySpecies = 10; - m_numKineticReactions = 1; - break; - - case ChemicalSystemType::carbonateAllEquilibrium: - m_numPrimarySpecies = 7; - m_numSecondarySpecies = 11; - m_numKineticReactions = 0; - break; - - case ChemicalSystemType::chainSerialAllKinetic: - m_numPrimarySpecies = 3; - m_numSecondarySpecies = 0; - m_numKineticReactions = 3; - break; - - case ChemicalSystemType::momasMedium: - m_numPrimarySpecies = 5; - m_numSecondarySpecies = 9; - m_numKineticReactions = 1; - break; - - default: - m_numPrimarySpecies = 5; - m_numSecondarySpecies = 7; - m_numKineticReactions = 0; - break; - } + using System = decltype( system ); + if( System::chemicalSystem == m_chemicalSystemType && System::activityModel == m_activityModelType ) + { + using ReactionParams = typename System::ReactionParamsType; + m_numPrimarySpecies = ReactionParams::numPrimarySpecies(); + m_numSecondarySpecies = ReactionParams::numSecondarySpecies(); + m_numKineticReactions = ReactionParams::numKineticReactions(); + } + } ); GEOS_THROW_IF( m_solventMassFraction <= 0.0 || m_solventMassFraction > 1.0, GEOS_FMT( "{}: invalid value of attribute '{}', expected a value in (0, 1]", diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index ee846698e6c..11a488f5c20 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -27,14 +27,11 @@ #include "constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp" #include "constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp" -#include "constitutive/HPCReact/src/reactions/geochemistry/GeochemicalSystems.hpp" -#include "constitutive/HPCReact/src/reactions/exampleSystems/BulkGeneric.hpp" -#include "constitutive/HPCReact/src/reactions/exampleSystems/ChainGeneric.hpp" -#include "constitutive/HPCReact/src/reactions/exampleSystems/MoMasBenchmark.hpp" +#include "constitutive/fluid/reactivefluid/ReactiveFluidSystemSelector.hpp" #include "constitutive/HPCReact/src/reactions/reactionsSystems/EquilibriumReactions.hpp" #include "constitutive/HPCReact/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp" -#include "constitutive/HPCReact/src/reactions/massActions/MassActions.hpp" #include +#include namespace geos { @@ -47,16 +44,6 @@ namespace reactivefluid using namespace hpcReact::reactionsSystems; -enum class ChemicalSystemType : integer -{ - carbonate, - carbonateAllEquilibrium, - ultramafic, - momasEasy, - momasMedium, - chainSerialAllKinetic -}; - template< typename BASE > class ReactiveSinglePhaseFluid : public BASE { @@ -125,7 +112,7 @@ class ReactiveSinglePhaseFluid : public BASE /** * @brief Kernel wrapper class for ReactiveSinglePhaseFluid. */ - template< typename REACTION_PARAMS_TYPE > + template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > class ReactionKernelWrapper { @@ -143,7 +130,8 @@ class ReactiveSinglePhaseFluid : public BASE integer const numPrimarySpecies, integer const numSecondarySpecies, integer const numKineticReactions, - REACTION_PARAMS_TYPE params ): + REACTION_PARAMS_TYPE params, + typename ACTIVITY_MODEL::Params activityParams ): m_numPrimarySpecies( numPrimarySpecies ), m_numSecondarySpecies( numSecondarySpecies ), m_numKineticReactions( numKineticReactions ), @@ -156,10 +144,11 @@ class ReactiveSinglePhaseFluid : public BASE m_kineticReactionRates( kineticReactionRates ), m_aggregateSpeciesRates( aggregateSpeciesRates ), m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations( dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ), - m_params( params ) + m_params( params ), + m_activityParams( activityParams ) {} - using EquilibriumReactionsType = hpcReact::reactionsSystems::EquilibriumReactions< real64, integer, localIndex >; + using EquilibriumReactionsType = hpcReact::reactionsSystems::EquilibriumReactions< real64, integer, localIndex, ACTIVITY_MODEL >; /** * @brief Get number of elements in this wrapper. @@ -168,14 +157,22 @@ class ReactiveSinglePhaseFluid : public BASE GEOS_HOST_DEVICE localIndex numElems() const { return m_secondarySpeciesConcentration.size( 0 ); } + /** + * @brief Speciate cell @p k at equilibrium. + * @return whether the equilibrium solve converged + */ GEOS_HOST_DEVICE - void updateEquilibriumReaction( localIndex const k, + bool updateEquilibriumReaction( localIndex const k, real64 const pressure, real64 const temperature, arraySlice1d< real64, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration ) const; + /** + * @brief Solve for the primary and secondary concentrations at the target aggregates. + * @return whether the solve converged + */ GEOS_HOST_DEVICE - void enforceEquilibrium( real64 const pressure, + bool enforceEquilibrium( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & targetPrimarySpeciesAggregateConcentration, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & initialPrimarySpeciesConcentration, @@ -231,113 +228,55 @@ class ReactiveSinglePhaseFluid : public BASE arrayView4d< real64, reactivefluid::USD_SPECIES_DC > m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations; REACTION_PARAMS_TYPE m_params; + + typename ACTIVITY_MODEL::Params m_activityParams; }; - std::variant< - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::ultramaficSystemType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::carbonateSystemType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::carbonateSystemAllEquilibriumType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::ChainGeneric::serialAllKineticType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::mediumCaseType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::easyCaseType > > - createReactionKernelWrapper() const + /// The kernel wrapper for one system of reactivefluid::ReactionSystemList. + template< typename SYSTEM > + using WrapperFor = ReactionKernelWrapper< typename SYSTEM::ReactionParamsType, typename SYSTEM::ActivityType >; + + /// @cond DO_NOT_DOCUMENT + template< typename LIST > + struct WrapperVariantHelper; + + template< typename ... SYSTEMS > + struct WrapperVariantHelper< std::variant< SYSTEMS... > > { - using namespace hpcReact::geochemistry; - using namespace hpcReact::MoMasBenchmark; - using namespace hpcReact::ChainGeneric; - switch( m_chemicalSystemType ) + using type = std::variant< WrapperFor< SYSTEMS > ... >; + }; + /// @endcond + + /// One alternative per system of reactivefluid::ReactionSystemList. + using ReactionKernelWrapperVariant = typename WrapperVariantHelper< reactivefluid::ReactionSystemList >::type; + + /** + * @brief Build the kernel wrapper for the chemical system and activity model this fluid was given. + * @return the wrapper, as the alternative of ReactionKernelWrapperVariant matching that pairing + * + * postInputInitialization has already rejected a pairing ReactionSystemList does not hold. + */ + ReactionKernelWrapperVariant createReactionKernelWrapper() const + { + std::optional< ReactionKernelWrapperVariant > wrapper; + + reactivefluid::forEachReactionSystem( [&]( auto system ) { - case ChemicalSystemType::ultramafic: - return ReactionKernelWrapper< ultramaficSystemType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - ultramaficSystem ); - - case ChemicalSystemType::carbonate: - return ReactionKernelWrapper< carbonateSystemType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - carbonateSystem ); - case ChemicalSystemType::carbonateAllEquilibrium: - return ReactionKernelWrapper< carbonateSystemAllEquilibriumType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - carbonateSystemAllEquilibrium ); - case ChemicalSystemType::chainSerialAllKinetic: - return ReactionKernelWrapper< serialAllKineticType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - serialAllKineticParams ); - case ChemicalSystemType::momasMedium: - return ReactionKernelWrapper< mediumCaseType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - mediumCaseParams ); - default: - return ReactionKernelWrapper< easyCaseType >( m_primarySpeciesAggregateConcentration, - m_primarySpeciesMobileAggregateConcentration, - m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, - m_initialPrimarySpeciesConcentration, - m_secondarySpeciesConcentration, - m_kineticReactionRates, - m_aggregateSpeciesRates, - m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, - m_numPrimarySpecies, - m_numSecondarySpecies, - m_numKineticReactions, - easyCaseParams ); - } + using System = decltype( system ); + if( System::chemicalSystem == m_chemicalSystemType && System::activityModel == m_activityModelType ) + { + wrapper.emplace( makeReactionKernelWrapper< WrapperFor< System > >( System::reactionParams(), + System::activityParams() ) ); + } + } ); + + return std::move( wrapper.value() ); } struct viewKeyStruct : ConstitutiveBase::viewKeyStruct { static constexpr char const * chemicalSystemNameString() { return "chemicalSystemType"; } + static constexpr char const * activityModelNameString() { return "activityModelType"; } static constexpr char const * solventMassFractionString() { return "solventMassFraction"; } }; @@ -347,6 +286,29 @@ class ReactiveSinglePhaseFluid : public BASE virtual void resizeFields( localIndex const size, localIndex const numPts ); + /** + * @brief Build one kernel wrapper for the given reaction system and activity model. + */ + template< typename WRAPPER_TYPE, typename REACTION_PARAMS_TYPE, typename ACTIVITY_PARAMS_TYPE > + WRAPPER_TYPE makeReactionKernelWrapper( REACTION_PARAMS_TYPE const & params, + ACTIVITY_PARAMS_TYPE const & activityParams ) const + { + return WRAPPER_TYPE( m_primarySpeciesAggregateConcentration, + m_primarySpeciesMobileAggregateConcentration, + m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, + m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, + m_initialPrimarySpeciesConcentration, + m_secondarySpeciesConcentration, + m_kineticReactionRates, + m_aggregateSpeciesRates, + m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, + m_numPrimarySpecies, + m_numSecondarySpecies, + m_numKineticReactions, + params, + activityParams ); + } + integer m_numPrimarySpecies; integer m_numSecondarySpecies; @@ -375,6 +337,8 @@ class ReactiveSinglePhaseFluid : public BASE ChemicalSystemType m_chemicalSystemType; + ActivityModelType m_activityModelType; + /// 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; @@ -386,41 +350,45 @@ using ReactiveCompressibleSinglePhaseFluid = ReactiveSinglePhaseFluid< Compressi using ReactiveThermalCompressibleSinglePhaseFluid = ReactiveSinglePhaseFluid< ThermalCompressibleSinglePhaseFluid >; template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE -inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +inline bool +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: updateEquilibriumReaction( localIndex const k, real64 const pressure, real64 const temperature, arraySlice1d< real64, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration ) const { - integer const numSecondarySpecies = m_numSecondarySpecies; + constexpr integer numSecondarySpecies = REACTION_PARAMS_TYPE::numSecondarySpecies(); - if( numSecondarySpecies > 0 ) + if constexpr ( numSecondarySpecies > 0 ) { - stackArray1d< real64, MAX_NUM_SPECIES > logSecondarySpeciesConcentration( numSecondarySpecies ); + stackArray1d< real64, numSecondarySpecies > logSecondarySpeciesConcentration( numSecondarySpecies ); - enforceEquilibrium( pressure, temperature, m_primarySpeciesAggregateConcentration[k][0], m_initialPrimarySpeciesConcentration[k][0], logPrimarySpeciesConcentration, - logSecondarySpeciesConcentration.toSlice() ); + bool const converged = enforceEquilibrium( pressure, temperature, m_primarySpeciesAggregateConcentration[k][0], + m_initialPrimarySpeciesConcentration[k][0], logPrimarySpeciesConcentration, + logSecondarySpeciesConcentration.toSlice() ); for( integer i=0; i < numSecondarySpecies; ++i ) { m_secondarySpeciesConcentration[k][0][i] = LvArray::math::exp( logSecondarySpeciesConcentration[i] ); } + + return converged; } else { GEOS_UNUSED_VAR( k, pressure, temperature, logPrimarySpeciesConcentration ); + return true; } } template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE -inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +inline bool +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: enforceEquilibrium( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & targetPrimarySpeciesAggregateConcentration, @@ -441,20 +409,21 @@ enforceEquilibrium( real64 const pressure, logPrimarySpeciesConcentration0[i] = LvArray::math::log( initialPrimarySpeciesConcentration[i] ); } - // 1. We enforce equilibrium - EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature, m_params, targetPrimarySpeciesAggregateConc, logPrimarySpeciesConcentration0, logPrimarySpeciesConcentration ); - - // 2. We calculate the secondary species concentration - hpcReact::massActions::calculateLogSecondarySpeciesConcentration< real64, - localIndex, - localIndex >( m_params, logPrimarySpeciesConcentration, logSecondarySpeciesConcentration ); + // Solve for the primary and secondary concentrations with equilibrium enforced at the target aggregates + return EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature, + m_params, + m_activityParams, + targetPrimarySpeciesAggregateConc, + logPrimarySpeciesConcentration0, + logPrimarySpeciesConcentration, + logSecondarySpeciesConcentration ); } template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: updateMixedReactionSystem( localIndex const k, real64 const pressure, real64 const temperature, @@ -489,10 +458,10 @@ updateMixedReactionSystem( localIndex const k, } template< typename BASE > -template< typename REACTION_PARAMS_TYPE > +template< typename REACTION_PARAMS_TYPE, typename ACTIVITY_MODEL > GEOS_HOST_DEVICE inline void -ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE >:: +ReactiveSinglePhaseFluid< BASE >::ReactionKernelWrapper< REACTION_PARAMS_TYPE, ACTIVITY_MODEL >:: computeAggregateConcentrationsAndRates( real64 const pressure, real64 const temperature, arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & logPrimarySpeciesConcentration, @@ -509,9 +478,10 @@ computeAggregateConcentrationsAndRates( real64 const pressure, { GEOS_UNUSED_VAR( pressure ); - MixedEquilibriumKineticReactions< real64, localIndex, localIndex, true >:: + MixedEquilibriumKineticReactions< real64, localIndex, localIndex, ACTIVITY_MODEL, true >:: updateMixedSystem( temperature, m_params, + m_activityParams, logPrimarySpeciesConcentration, surfaceArea, logSecondarySpeciesConcentration, @@ -525,14 +495,6 @@ computeAggregateConcentrationsAndRates( real64 const pressure, dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); } -ENUM_STRINGS( ChemicalSystemType, - "carbonate", - "carbonateAllEquilibrium", - "ultramafic", - "momasEasy", - "momasMedium", - "chainSerialAllKinetic" ); - } // namespace reactivefluid } // namespace constitutive diff --git a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp index c4b50dfab30..0875378f4c8 100644 --- a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp @@ -562,13 +562,15 @@ TEST_F( SinglePhaseReactiveTransportTest, jacobianNumericalCheck_accumulationBal real64 const perturb = std::sqrt( eps ); real64 const tol = 1e-6; // 1% error margin + real64 const accumulationDt = 1e-4; + DomainPartition & domain = state.getProblemManager().getDomainPartition(); testNumericalJacobian( *solver, domain, perturb, tol, [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - solver->assembleAccumulationTermsInMassBalanceAndSpeciesAmountEqs( dt, domain, solver->getDofManager(), localMatrix, localRhs ); + solver->assembleAccumulationTermsInMassBalanceAndSpeciesAmountEqs( accumulationDt, domain, solver->getDofManager(), localMatrix, localRhs ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 0d2632b065a..dbea76beb0b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -19,6 +19,7 @@ #include "SinglePhaseReactiveTransport.hpp" +#include "common/format/LogPart.hpp" #include "constitutive/ConstitutiveManager.hpp" #include "constitutive/ConstitutivePassThru.hpp" #include "constitutive/diffusion/DiffusionFields.hpp" @@ -106,6 +107,13 @@ SinglePhaseReactiveTransport::SinglePhaseReactiveTransport( const string & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Array to store the indices of immobile species. Default is {}, which indicates no immobile species." ); + this->registerWrapper( viewKeyStruct::maxAbsoluteLogConcChangeString(), &m_maxAbsoluteLogConcChange ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 2.0 * 2.302585092994046 ). // two ln10 units + setDescription( "Maximum (absolute) change in the natural log of a primary species concentration " + "in a Newton iteration. Zero or less disables the scaling." ); + addLogLevel< logInfo::BoundaryConditions >(); } @@ -773,6 +781,10 @@ void SinglePhaseReactiveTransport::updateSurfaceArea( ElementSubRegionBase & sub void SinglePhaseReactiveTransport::initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { + LogPart equilibriumLog( "Initial Chemical Equilibrium Enforcement", + MpiWrapper::commRank() == 0 && isLogLevelActive< logInfo::Convergence >( getLogLevel() ) ); + equilibriumLog.begin(); + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) { @@ -813,6 +825,8 @@ void SinglePhaseReactiveTransport::initializeFluidState( MeshLevel & mesh, strin diffusionMaterial.initializeTemperatureState( temperature ); } } ); + + equilibriumLog.end(); } void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegionBase & subRegion ) const @@ -823,6 +837,8 @@ void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegi arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >(); arrayView2d< real64, compflow::USD_COMP > const logPrimaryConc = subRegion.getField< fields::flow::logPrimarySpeciesConcentration >(); + bool converged = true; + if( m_isThermal ) { reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid & fluid = @@ -830,7 +846,7 @@ void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegi constitutive::constitutiveUpdatePassThru( fluid, [&]( auto & castedFluid ) { - singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); + converged = singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); } ); fluid.saveConvergedState(); @@ -842,11 +858,24 @@ void SinglePhaseReactiveTransport::initializeEquilibriumReaction( ElementSubRegi constitutive::constitutiveUpdatePassThru( fluid, [&]( auto & castedFluid ) { - singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); + converged = singlePhaseReactiveBaseKernels::EquilibriumReactionUpdateKernel::launch( castedFluid, pres, temp, logPrimaryConc ); } ); fluid.saveConvergedState(); } + + // Report the cell count and the status over all ranks, not just the one doing the logging. + globalIndex const numCells = MpiWrapper::sum< globalIndex >( subRegion.getNumberOfLocalIndices() ); + integer const allConverged = MpiWrapper::min< integer >( converged ? 1 : 0 ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Convergence, + GEOS_FMT( "{}: initial equilibrium speciation on {} ({} cells): {}", + getName(), subRegion.getName(), numCells, + allConverged ? "converged" : "NOT converged" ) ); + + GEOS_ERROR_IF( !converged, + GEOS_FMT( "{}: the initial equilibrium speciation did not converge.", + subRegion.getDataContext() ) ); } void SinglePhaseReactiveTransport::initializePostInitialConditionsPreSubGroups() @@ -1357,6 +1386,56 @@ real64 SinglePhaseReactiveTransport::calculateResidualNorm( real64 const & GEOS_ return residualNorm; } +real64 SinglePhaseReactiveTransport::scalingForSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) +{ + GEOS_MARK_FUNCTION; + + // Pressure, and temperature when thermal, are scaled by the base. + real64 scalingFactor = SinglePhaseBase::scalingForSystemSolution( domain, dofManager, localSolution ); + + string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + integer const speciesOffset = m_isThermal ? 2 : 1; + real64 maxDeltaLogConc = 0.0; + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + auto const subRegionData = + singlePhaseReactiveBaseKernels::SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + dofNumber, + ghostRank, + m_numDofPerCell, + speciesOffset, + m_numPrimarySpecies, + m_maxAbsoluteLogConcChange ); + + scalingFactor = std::min( scalingFactor, subRegionData.first ); + maxDeltaLogConc = std::max( maxDeltaLogConc, subRegionData.second ); + } ); + } ); + + scalingFactor = MpiWrapper::min( scalingFactor ); + maxDeltaLogConc = MpiWrapper::max( maxDeltaLogConc ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, + GEOS_FMT( " {}: Max log concentration change = {:.4g} (before scaling)", + getName(), maxDeltaLogConc ) ); + + return scalingFactor; +} + void SinglePhaseReactiveTransport::applySystemSolution( DofManager const & dofManager, arrayView1d< real64 const > const & localSolution, real64 const scalingFactor, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp index 544b0f612f3..4fe5341ae1b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp @@ -28,6 +28,7 @@ #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluidUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp" @@ -127,6 +128,18 @@ class SinglePhaseReactiveTransport : public SinglePhaseBase DofManager const & dofManager, arrayView1d< real64 const > const & localRhs ) override; + /** + * @brief Compute the factor the Newton update is scaled by before it is applied. + * @param domain the physical domain + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localSolution the solution vector + * @return the smaller of the pressure scaling from the base and the log concentration scaling + */ + virtual real64 + scalingForSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) override; + virtual void applySystemSolution( DofManager const & dofManager, arrayView1d< real64 const > const & localSolution, @@ -271,10 +284,15 @@ class SinglePhaseReactiveTransport : public SinglePhaseBase static constexpr char const * isUpdateReactivePorosityString() { return "isUpdateReactivePorosity"; } static constexpr char const * isUpdateSurfaceAreaString() { return "isUpdateSurfaceArea"; } static constexpr char const * immobilePrimarySpeciesIndicesString() { return "immobilePrimarySpeciesIndices"; } + + static constexpr char const * maxAbsoluteLogConcChangeString() { return "maxAbsoluteLogConcentrationChange"; } }; protected: + /// Maximum (absolute) change in log primary species concentration in a Newton iteration, in ln units + real64 m_maxAbsoluteLogConcChange; + /// the number of primary species in the fluid integer m_numPrimarySpecies; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp index 79496d888fe..8d7ea5efb36 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ReactionUpdateKernel.hpp @@ -36,28 +36,45 @@ namespace singlePhaseReactiveBaseKernels struct EquilibriumReactionUpdateKernel { + /** + * @return whether the equilibrium solve converged in every cell + */ template< typename REACTION_WRAPPER_TYPE > - static void helper( REACTION_WRAPPER_TYPE const & reactionWrapper, + static bool helper( REACTION_WRAPPER_TYPE const & reactionWrapper, arrayView1d< real64 const > const & pres, arrayView1d< real64 const > const & temp, arrayView2d< real64, compflow::USD_COMP > const logPrimaryConc ) { + RAJA::ReduceMax< parallelDeviceReduce, integer > anyNotConverged( 0 ); + forAll< parallelDevicePolicy<> >( reactionWrapper.numElems(), [=] GEOS_HOST_DEVICE ( localIndex const k ) { - reactionWrapper.updateEquilibriumReaction( k, pres[k], temp[k], logPrimaryConc[k] ); + if( !reactionWrapper.updateEquilibriumReaction( k, pres[k], temp[k], logPrimaryConc[k] ) ) + { + anyNotConverged.max( 1 ); + } } ); + + return anyNotConverged.get() == 0; } + /** + * @return whether the equilibrium solve converged in every cell, for the caller to report + */ template< typename REACTIVE_FLUID > - static void launch( REACTIVE_FLUID const & fluid, + static bool launch( REACTIVE_FLUID const & fluid, arrayView1d< real64 const > const & pres, arrayView1d< real64 const > const & temp, arrayView2d< real64, compflow::USD_COMP > const logPrimaryConc ) { + bool converged = true; + std::visit( [&]( auto const & reactionWrapper ) { - helper( reactionWrapper, pres, temp, logPrimaryConc ); + converged = helper( reactionWrapper, pres, temp, logPrimaryConc ); }, fluid.createReactionKernelWrapper()); + + return converged; } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp index c9fc35acce1..9eccb735f01 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp @@ -77,7 +77,7 @@ class IsothermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNo // step 2: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); real64 const valAmount = LvArray::math::abs( m_localResidual[stack.localRow + idof + 1] ) / speciesAmountNormalizer; if( valAmount > stack.localValue[1] ) { @@ -99,7 +99,7 @@ class IsothermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNo // step 2: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); stack.localValue[1] += m_localResidual[stack.localRow + idof + 1] * m_localResidual[stack.localRow + idof + 1]; stack.localNormalizer[1] += speciesAmountNormalizer; @@ -185,7 +185,7 @@ class ThermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNormK // step 3: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); real64 const valAmount = LvArray::math::abs( m_localResidual[stack.localRow + idof + 2] ) / speciesAmountNormalizer; if( valAmount > stack.localValue[2] ) { @@ -214,7 +214,7 @@ class ThermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNormK // step 3: species amount residuals for( integer idof = 0; idof < m_numPrimarySpecies; ++idof ) { - real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] ); + real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_primarySpeciesAggregateMole_n[ei][idof] ) ); stack.localValue[2] += m_localResidual[stack.localRow + idof + 2] * m_localResidual[stack.localRow + idof + 2]; stack.localNormalizer[2] += speciesAmountNormalizer; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp new file mode 100644 index 00000000000..241b035df39 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SolutionScalingKernel.hpp @@ -0,0 +1,94 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SolutionScalingKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_SOLUTIONSCALINGKERNEL_HPP +#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_SOLUTIONSCALINGKERNEL_HPP + +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" + +namespace geos +{ + +namespace singlePhaseReactiveBaseKernels +{ + +/** + * @brief Bound the change in log primary species concentration over a Newton iteration. + */ +struct SolutionScalingKernel +{ + /** + * @brief Compute the factor that keeps the largest log concentration change within the bound. + * @tparam POLICY the execution policy + * @param localSolution the solution vector + * @param rankOffset the offset of this rank in the solution vector + * @param dofNumber the first degree of freedom of each element + * @param ghostRank the ghost rank of each element + * @param numDofPerCell the number of degrees of freedom per cell + * @param speciesOffset the offset of the first species within a cell, 1 or 2 when thermal + * @param numPrimarySpecies the number of primary species + * @param maxAbsoluteLogConcChange the bound, in ln units. Zero or less disables the scaling. + * @return the scaling factor and the largest log concentration change seen, before scaling + */ + template< typename POLICY > + static std::pair< real64, real64 > launch( arrayView1d< real64 const > const & localSolution, + globalIndex const rankOffset, + arrayView1d< globalIndex const > const & dofNumber, + arrayView1d< integer const > const & ghostRank, + integer const numDofPerCell, + integer const speciesOffset, + integer const numPrimarySpecies, + real64 const maxAbsoluteLogConcChange ) + { + GEOS_UNUSED_VAR( numDofPerCell ); + + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > scalingFactor( 1.0 ); + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaLogConc( 0.0 ); + + forAll< POLICY >( dofNumber.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + if( ghostRank[ei] < 0 && dofNumber[ei] >= 0 ) + { + localIndex const lid = dofNumber[ei] - rankOffset; + + for( integer ic = 0; ic < numPrimarySpecies; ++ic ) + { + real64 const absLogConcChange = LvArray::math::abs( localSolution[lid + speciesOffset + ic] ); + maxDeltaLogConc.max( absLogConcChange ); + + if( maxAbsoluteLogConcChange > 0.0 && absLogConcChange > maxAbsoluteLogConcChange ) + { + scalingFactor.min( maxAbsoluteLogConcChange / absLogConcChange ); + } + } + } + + } ); + + return { scalingFactor.get(), maxDeltaLogConc.get() }; + } + +}; + +} // namespace singlePhaseReactiveBaseKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_SOLUTIONSCALINGKERNEL_HPP diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index c719a937ad1..35550fb8ce0 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -5355,6 +5355,8 @@ Information output from lower logLevels is added with the desired log level 2 - The summary of declared fields and coupling--> + + @@ -8811,6 +8813,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + @@ -8830,8 +8834,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia - - + + + + + + + @@ -8896,6 +8905,8 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + +