diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 746f223b649..2cb02daa192 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3495-9552-257c87e + baseline: integratedTests/baseline_integratedTests-pr3228-9676-61994fe allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index c5617ba8ed9..ea19f4ac5cd 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3228 (2024-01-15) +===================== +deltaVolume added in multiphase. + PR #3495 (2024-01-08) ===================== Add missing logic to support switching from fixed mass rate injection rate constraint to max injection pressure. diff --git a/examples/ObjectCatalog/Base.hpp b/examples/ObjectCatalog/Base.hpp index 2075058c05a..c8e320afe48 100644 --- a/examples/ObjectCatalog/Base.hpp +++ b/examples/ObjectCatalog/Base.hpp @@ -5,7 +5,7 @@ * 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) 2023-2024 Chevron * Copyright (c) 2019- GEOS/GEOSX Contributors * All rights reserved * @@ -46,7 +46,7 @@ class Parameter class Base { public: - Base( int junk, double const & junk2, Parameter& pbv ) + Base( int junk, double const & junk2, Parameter & pbv ) { GEOS_LOG( "calling Base constructor with arguments (" << junk << " " << junk2 << ")" ); } @@ -56,14 +56,20 @@ class Base GEOS_LOG( "calling Base destructor" ); } - using CatalogInterface = dataRepository::CatalogInterface< Base, int, double const &, Parameter& >; - static CatalogInterface::CatalogType& getCatalog() + using CatalogInterface = dataRepository::CatalogInterface< Base, int, double const &, Parameter & >; + static CatalogInterface::CatalogType & getCatalog() { static CatalogInterface::CatalogType catalog; return catalog; } virtual std::string const getName() const = 0; + + static DataContext const & getDataContext() + { + static DataFileContext const context = DataFileContext( "Base Test Class", __FILE__, __LINE__ ); + return context; + } }; #endif diff --git a/examples/ObjectCatalog/main.cpp b/examples/ObjectCatalog/main.cpp index 2a9ed5ba436..88c52c48d95 100644 --- a/examples/ObjectCatalog/main.cpp +++ b/examples/ObjectCatalog/main.cpp @@ -5,7 +5,7 @@ * 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) 2023-2024 Chevron * Copyright (c) 2019- GEOS/GEOSX Contributors * All rights reserved * @@ -31,10 +31,14 @@ int main( int argc, char *argv[] ) GEOS_LOG( "Attempting to create a Derived1 object" ); - std::unique_ptr derived1 = Base::CatalogInterface::Factory( "derived1", junk, junk2, param ); + std::unique_ptr< Base > derived1 = Base::CatalogInterface::factory( "derived1", + Base::getDataContext(), + junk, junk2, param ); GEOS_LOG( "Attempting to create a Derived2 object" ); - std::unique_ptr derived2 = Base::CatalogInterface::Factory( "derived2", junk, junk3, param ); + std::unique_ptr< Base > derived2 = Base::CatalogInterface::factory( "derived2", + Base::getDataContext(), + junk, junk3, param ); - Base::CatalogInterface::catalog_cast( *(derived2.get())); + Base::CatalogInterface::catalog_cast< Derived1 >( *(derived2.get())); GEOS_LOG( "EXITING MAIN" ); } diff --git a/inputFiles/poromechanicsFractures/co2flash.txt b/inputFiles/poromechanicsFractures/co2flash.txt new file mode 100644 index 00000000000..8de4352642b --- /dev/null +++ b/inputFiles/poromechanicsFractures/co2flash.txt @@ -0,0 +1 @@ +FlashModel CO2Solubility 1e6 10e7 5e4 367.15 369.15 1 0 diff --git a/inputFiles/poromechanicsFractures/elevation.txt b/inputFiles/poromechanicsFractures/elevation.txt new file mode 100644 index 00000000000..9f9ff9d4f67 --- /dev/null +++ b/inputFiles/poromechanicsFractures/elevation.txt @@ -0,0 +1,2 @@ +0 +1e4 \ No newline at end of file diff --git a/inputFiles/poromechanicsFractures/initCO2CompFrac.txt b/inputFiles/poromechanicsFractures/initCO2CompFrac.txt new file mode 100644 index 00000000000..f3c5e203d6f --- /dev/null +++ b/inputFiles/poromechanicsFractures/initCO2CompFrac.txt @@ -0,0 +1,2 @@ +0.001 +0.001 diff --git a/inputFiles/poromechanicsFractures/initTemp.txt b/inputFiles/poromechanicsFractures/initTemp.txt new file mode 100644 index 00000000000..da7b374fdb7 --- /dev/null +++ b/inputFiles/poromechanicsFractures/initTemp.txt @@ -0,0 +1,2 @@ +368.15 +368.15 diff --git a/inputFiles/poromechanicsFractures/initWaterCompFrac.txt b/inputFiles/poromechanicsFractures/initWaterCompFrac.txt new file mode 100644 index 00000000000..c88775171f6 --- /dev/null +++ b/inputFiles/poromechanicsFractures/initWaterCompFrac.txt @@ -0,0 +1,2 @@ +0.999 +0.999 diff --git a/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_base.xml b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_base.xml new file mode 100644 index 00000000000..fc1798c6851 --- /dev/null +++ b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_base.xml @@ -0,0 +1,207 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_smoke.xml b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_smoke.xml new file mode 100644 index 00000000000..cb87c0a50e3 --- /dev/null +++ b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_smoke.xml @@ -0,0 +1,130 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/poromechanicsFractures/phaseVolFraction_gas.txt b/inputFiles/poromechanicsFractures/phaseVolFraction_gas.txt new file mode 100644 index 00000000000..e2d1da1f128 --- /dev/null +++ b/inputFiles/poromechanicsFractures/phaseVolFraction_gas.txt @@ -0,0 +1,21 @@ +0 +0.05 +0.1 +0.15 +0.2 +0.25 +0.3 +0.35 +0.4 +0.45 +0.5 +0.55 +0.6 +0.65 +0.7 +0.75 +0.8 +0.85 +0.9 +0.95 +1.0 diff --git a/inputFiles/poromechanicsFractures/phaseVolFraction_water.txt b/inputFiles/poromechanicsFractures/phaseVolFraction_water.txt new file mode 100644 index 00000000000..e2d1da1f128 --- /dev/null +++ b/inputFiles/poromechanicsFractures/phaseVolFraction_water.txt @@ -0,0 +1,21 @@ +0 +0.05 +0.1 +0.15 +0.2 +0.25 +0.3 +0.35 +0.4 +0.45 +0.5 +0.55 +0.6 +0.65 +0.7 +0.75 +0.8 +0.85 +0.9 +0.95 +1.0 diff --git a/inputFiles/poromechanicsFractures/pvtgas.txt b/inputFiles/poromechanicsFractures/pvtgas.txt new file mode 100644 index 00000000000..9da1273477c --- /dev/null +++ b/inputFiles/poromechanicsFractures/pvtgas.txt @@ -0,0 +1,2 @@ +DensityFun SpanWagnerCO2Density 1e6 10e7 5e4 367.15 369.15 1 +ViscosityFun FenghourCO2Viscosity 1e6 10e7 5e4 367.15 369.15 1 diff --git a/inputFiles/poromechanicsFractures/pvtliquid.txt b/inputFiles/poromechanicsFractures/pvtliquid.txt new file mode 100644 index 00000000000..b0bfb5fc76e --- /dev/null +++ b/inputFiles/poromechanicsFractures/pvtliquid.txt @@ -0,0 +1,2 @@ +DensityFun PhillipsBrineDensity 1e6 10e7 5e4 367.15 369.15 1 0 +ViscosityFun PhillipsBrineViscosity 0 diff --git a/inputFiles/poromechanicsFractures/relPerm_gas.txt b/inputFiles/poromechanicsFractures/relPerm_gas.txt new file mode 100644 index 00000000000..d6c4e7c035e --- /dev/null +++ b/inputFiles/poromechanicsFractures/relPerm_gas.txt @@ -0,0 +1,21 @@ +0 +0 +0.0118 +0.0333 +0.0612 +0.0943 +0.1318 +0.1732 +0.2183 +0.2667 +0.3182 +0.3727 +0.4300 +0.4899 +0.5524 +0.6173 +0.6847 +0.7542 +0.8261 +0.9000 +0.9760 diff --git a/inputFiles/poromechanicsFractures/relPerm_water.txt b/inputFiles/poromechanicsFractures/relPerm_water.txt new file mode 100644 index 00000000000..d6c4e7c035e --- /dev/null +++ b/inputFiles/poromechanicsFractures/relPerm_water.txt @@ -0,0 +1,21 @@ +0 +0 +0.0118 +0.0333 +0.0612 +0.0943 +0.1318 +0.1732 +0.2183 +0.2667 +0.3182 +0.3727 +0.4300 +0.4899 +0.5524 +0.6173 +0.6847 +0.7542 +0.8261 +0.9000 +0.9760 diff --git a/src/coreComponents/common/format/StringUtilities.hpp b/src/coreComponents/common/format/StringUtilities.hpp index 79b0333d0a2..b79db4b1168 100644 --- a/src/coreComponents/common/format/StringUtilities.hpp +++ b/src/coreComponents/common/format/StringUtilities.hpp @@ -88,7 +88,7 @@ string join( CONTAINER const & container, S const & delim = S() ) * @return a string containing input values concatenated with a delimiter */ template< typename IT, typename S, typename LAMBDA > -string joinLamda( IT first, IT last, S const & delim, LAMBDA formattingFunc ) +string joinLambda( IT first, IT last, S const & delim, LAMBDA formattingFunc ) { if( first == last ) { @@ -114,9 +114,9 @@ string joinLamda( IT first, IT last, S const & delim, LAMBDA formattingFunc ) * @return a string containing input values concatenated with a delimiter */ template< typename CONTAINER, typename S, typename LAMBDA > -string joinLamda( CONTAINER const & container, S const & delim, LAMBDA formattingFunc ) +string joinLambda( CONTAINER const & container, S const & delim, LAMBDA formattingFunc ) { - return joinLamda( std::begin( container ), std::end( container ), delim, formattingFunc ); + return joinLambda( std::begin( container ), std::end( container ), delim, formattingFunc ); } /** diff --git a/src/coreComponents/constitutive/ConstitutiveBase.cpp b/src/coreComponents/constitutive/ConstitutiveBase.cpp index 03bf6852ecc..7323a649605 100644 --- a/src/coreComponents/constitutive/ConstitutiveBase.cpp +++ b/src/coreComponents/constitutive/ConstitutiveBase.cpp @@ -83,7 +83,8 @@ ConstitutiveBase::deliverClone( string const & name, Group * const parent ) const { std::unique_ptr< ConstitutiveBase > - newModel = ConstitutiveBase::CatalogInterface::factory( this->getCatalogName(), name, parent ); + newModel = ConstitutiveBase::CatalogInterface::factory( this->getCatalogName(), getDataContext(), + name, parent ); newModel->forWrappers( [&]( WrapperBase & wrapper ) { diff --git a/src/coreComponents/constitutive/ConstitutiveManager.cpp b/src/coreComponents/constitutive/ConstitutiveManager.cpp index 3babd88203e..913f0898971 100644 --- a/src/coreComponents/constitutive/ConstitutiveManager.cpp +++ b/src/coreComponents/constitutive/ConstitutiveManager.cpp @@ -44,7 +44,8 @@ ConstitutiveManager::~ConstitutiveManager() Group * ConstitutiveManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< ConstitutiveBase > material = ConstitutiveBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< ConstitutiveBase > material = + ConstitutiveBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return ®isterGroup< ConstitutiveBase >( childName, std::move( material ) ); } diff --git a/src/coreComponents/dataRepository/DataContext.cpp b/src/coreComponents/dataRepository/DataContext.cpp index 87e9f791750..d34ab4b72c1 100644 --- a/src/coreComponents/dataRepository/DataContext.cpp +++ b/src/coreComponents/dataRepository/DataContext.cpp @@ -25,7 +25,7 @@ namespace dataRepository { -DataContext::DataContext( string const & targetName ): +DataContext::DataContext( string_view targetName ): m_targetName( targetName ) {} @@ -35,12 +35,12 @@ std::ostream & operator<<( std::ostream & os, DataContext const & ctx ) return os; } -DataContext::ToStringInfo::ToStringInfo( string const & targetName, string const & filePath, size_t line ): +DataContext::ToStringInfo::ToStringInfo( string_view targetName, string_view filePath, size_t line ): m_targetName( targetName ), m_filePath( filePath ), m_line( line ) {} -DataContext::ToStringInfo::ToStringInfo( string const & targetName ): +DataContext::ToStringInfo::ToStringInfo( string_view targetName ): m_targetName( targetName ) {} @@ -83,6 +83,15 @@ DataFileContext::DataFileContext( xmlWrapper::xmlNode const & targetNode, m_offset( attPos.offset ) {} +DataFileContext::DataFileContext( string_view targetName, string_view file, size_t line ): + DataContext( targetName ), + m_typeName( "C++ Source File" ), + m_filePath( file ), + m_line( line ), + m_offsetInLine( 0 ), + m_offset( 0 ) +{} + string DataFileContext::toString() const { if( m_line != xmlWrapper::xmlDocument::npos ) @@ -95,7 +104,7 @@ string DataFileContext::toString() const } else { - return GEOS_FMT( "{} (Source file not found)", m_targetName ); + return GEOS_FMT( "{} ({})", m_targetName, ( m_filePath.empty() ? "Source file not found" : m_filePath ) ); } } diff --git a/src/coreComponents/dataRepository/DataContext.hpp b/src/coreComponents/dataRepository/DataContext.hpp index 46daeb38506..ee3abd0f1c1 100644 --- a/src/coreComponents/dataRepository/DataContext.hpp +++ b/src/coreComponents/dataRepository/DataContext.hpp @@ -47,7 +47,7 @@ class DataContext * @brief Construct a new DataContext object. * @param targetName the target object name */ - DataContext( string const & targetName ); + DataContext( string_view targetName ); /** * @brief Destroy the DataContext object @@ -95,12 +95,12 @@ class DataContext * @param filePath the input file path where the target is declared. * @param line the line in the file where the target is declared. */ - ToStringInfo( string const & targetName, string const & filePath, size_t line ); + ToStringInfo( string_view targetName, string_view filePath, size_t line ); /** * @brief Construct a new ToStringInfo object from a DataContext that has no input file info. * @param targetName the target name. */ - ToStringInfo( string const & targetName ); + ToStringInfo( string_view targetName ); /** * @return true if a location has been found to declare the target in an input file. */ @@ -142,6 +142,14 @@ class DataFileContext final : public DataContext DataFileContext( xmlWrapper::xmlNode const & targetNode, xmlWrapper::xmlAttribute const & att, xmlWrapper::xmlAttributePos const & attPos ); + /** + * @brief Constructs the file context of a Group from a C++ source file. + * @param targetName The name of the target Group. + * @param file The name of the source file. + * @param line The line number in the source file. + */ + DataFileContext( string_view targetName, string_view file, size_t line ); + /** * @return the target object name followed by the the file and line declaring it. */ diff --git a/src/coreComponents/dataRepository/Group.cpp b/src/coreComponents/dataRepository/Group.cpp index 440178c2a7d..8665b5aa9c7 100644 --- a/src/coreComponents/dataRepository/Group.cpp +++ b/src/coreComponents/dataRepository/Group.cpp @@ -135,18 +135,64 @@ string Group::getPath() const return noProblem.empty() ? "/" : noProblem; } +string Group::processInputName( xmlWrapper::xmlNode const & targetNode, + xmlWrapper::xmlNodePos const & targetNodePos, + string_view parentNodeName, + xmlWrapper::xmlNodePos const & parentNodePos, + std::set< string > & siblingNames ) +{ + GEOS_THROW_IF( targetNode.type() != xmlWrapper::xmlNodeType::node_element, + GEOS_FMT( "Error in node named \"{}\" ({}): GEOS XML nodes cannot contain " + "text data nor anything but XML nodes.\nErroneous content: \"{}\"\n", + parentNodeName, parentNodePos.toString(), + stringutilities::trimSpaces( targetNode.value() ) ), + InputError ); + + string targetNodeName; + try + { // read & validate the name attribute + xmlWrapper::readAttributeAsType( targetNodeName, "name", + rtTypes::getTypeRegex< string >( rtTypes::CustomTypes::groupName ), + targetNode, string( "" ) ); + } catch( std::exception const & ex ) + { + xmlWrapper::processInputException( ex, "name", targetNode, targetNodePos ); + } + + if( targetNodeName.empty() ) + { // if there is no name attribute, we use the node tag as a node name + targetNodeName = targetNode.name(); + } + else + { // Make sure names are not duplicated by checking all previous siblings + GEOS_THROW_IF( siblingNames.count( targetNodeName ) != 0, + GEOS_FMT( "Error at node named \"{}\" ({}): " + "An XML block cannot contain children with duplicated names.\n", + targetNodeName, targetNodePos.toString() ), + InputError ); + siblingNames.insert( targetNodeName ); + } + + return targetNodeName; +} + void Group::processInputFileRecursive( xmlWrapper::xmlDocument & xmlDocument, xmlWrapper::xmlNode & targetNode ) { - xmlWrapper::xmlNodePos nodePos = xmlDocument.getNodePosition( targetNode ); - processInputFileRecursive( xmlDocument, targetNode, nodePos ); + xmlWrapper::xmlNodePos targetNodePos = xmlDocument.getNodePosition( targetNode ); + processInputFileRecursive( xmlDocument, targetNode, targetNodePos ); } void Group::processInputFileRecursive( xmlWrapper::xmlDocument & xmlDocument, xmlWrapper::xmlNode & targetNode, - xmlWrapper::xmlNodePos const & nodePos ) + xmlWrapper::xmlNodePos const & targetNodePos ) { xmlDocument.addIncludedXML( targetNode ); + if( targetNodePos.isFound() ) + { + m_dataContext = std::make_unique< DataFileContext >( targetNode, targetNodePos ); + } + // Handle the case where the node was imported from a different input file // Set the path prefix to make sure all relative Path variables are interpreted correctly string const oldPrefix = std::string( Path::getPathPrefix() ); @@ -157,40 +203,12 @@ void Group::processInputFileRecursive( xmlWrapper::xmlDocument & xmlDocument, } // Loop over the child nodes of the targetNode - array1d< string > childNames; + std::set< string > childNames; for( xmlWrapper::xmlNode childNode : targetNode.children() ) { xmlWrapper::xmlNodePos childNodePos = xmlDocument.getNodePosition( childNode ); - - // Get the child tag and name - string childName; - - try - { - xmlWrapper::readAttributeAsType( childName, "name", - rtTypes::getTypeRegex< string >( rtTypes::CustomTypes::groupName ), - childNode, string( "" ) ); - } catch( std::exception const & ex ) - { - xmlWrapper::processInputException( ex, "name", childNode, childNodePos ); - } - - if( childName.empty() ) - { - childName = childNode.name(); - } - else - { - // Make sure child names are not duplicated - GEOS_ERROR_IF( std::find( childNames.begin(), childNames.end(), childName ) != childNames.end(), - GEOS_FMT( "Error: An XML block cannot contain children with duplicated names.\n" - "Error detected at node {} with name = {} ({}:l.{})", - childNode.path(), childName, xmlDocument.getFilePath(), - xmlDocument.getNodePosition( childNode ).line ) ); - - childNames.emplace_back( childName ); - } - + string const childName = processInputName( childNode, childNodePos, + getName(), targetNodePos, childNames ); // Create children Group * newChild = createChild( childNode.name(), childName ); if( newChild == nullptr ) @@ -203,7 +221,7 @@ void Group::processInputFileRecursive( xmlWrapper::xmlDocument & xmlDocument, } } - processInputFile( targetNode, nodePos ); + processInputFile( targetNode, targetNodePos ); // Restore original prefix once the node is processed Path::setPathPrefix( oldPrefix ); @@ -212,12 +230,6 @@ void Group::processInputFileRecursive( xmlWrapper::xmlDocument & xmlDocument, void Group::processInputFile( xmlWrapper::xmlNode const & targetNode, xmlWrapper::xmlNodePos const & nodePos ) { - if( nodePos.isFound() ) - { - m_dataContext = std::make_unique< DataFileContext >( targetNode, nodePos ); - } - - std::set< string > processedAttributes; for( std::pair< string const, WrapperBase * > & pair : m_wrappers ) { @@ -233,10 +245,10 @@ void Group::processInputFile( xmlWrapper::xmlNode const & targetNode, if( !xmlWrapper::isFileMetadataAttribute( attributeName ) ) { GEOS_THROW_IF( processedAttributes.count( attributeName ) == 0, - GEOS_FMT( "XML Node at '{}' with name={} contains unused attribute '{}'.\n" + GEOS_FMT( "Error in {}: XML Node at '{}' contains unused attribute '{}'.\n" "Valid attributes are:\n{}\nFor more details, please refer to documentation at:\n" "http://geosx-geosx.readthedocs-hosted.com/en/latest/docs/sphinx/userGuide/Index.html", - targetNode.path(), m_dataContext->toString(), attributeName, + getDataContext(), targetNode.path(), attributeName, dumpInputOptions() ), InputError ); } @@ -264,11 +276,10 @@ void Group::registerDataOnMeshRecursive( Group & meshBodies ) Group * Group::createChild( string const & childKey, string const & childName ) { - GEOS_ERROR_IF( !(CatalogInterface::hasKeyName( childKey )), - "KeyName ("<) - * - * @tparam T The type of the Group to add/register. This should be a type that derives from Group. - * @tparam TBASE The type whose type catalog will be used to look up the new sub-group type - * @param[in] name The name of the group to use as a string key. - * @param[in] catalogName The catalog name of the new type. - * @return A pointer to the newly registered Group. - * - * Creates and registers a Group or class derived from Group as a subgroup of this Group. - */ - template< typename T = Group, typename TBASE = Group > - T & registerGroup( string const & name, string const & catalogName ) - { - std::unique_ptr< TBASE > newGroup = TBASE::CatalogInterface::Factory( catalogName, name, this ); - return registerGroup< T >( name, std::move( newGroup ) ); - } - /** * @brief Removes a child group from this group. * @param name the name of the child group to remove from this group. @@ -769,6 +751,26 @@ class Group */ void postRestartInitializationRecursive(); + /** + * @return The validated name for a given Group from its xml value: If the node has a "name" + * attribute, it is validated after the `groupName` rtType regex, and its value is + * returned. Else if the Group name is not "Required", the node tag name is used. + * @param targetNode The XML node whose name is to be processed. It throws if not of element type. + * @param targetNodePos The position of the target node within the XML document. + * @param parentNodeName The name of the parent node, used for error reporting. + * @param parentNodePos The position of the parent node, used for error reporting. + * @param siblingNames A set containing the names of sibling nodes (to verify that there are no + * duplicates). The function will populate this set if the attribute name is + * used and if no error is found. + * @throws InputError if the node type is not an xml element or if there are duplicate names + * among xml siblings. + */ + static string processInputName( xmlWrapper::xmlNode const & targetNode, + xmlWrapper::xmlNodePos const & targetNodePos, + string_view parentNodeName, + xmlWrapper::xmlNodePos const & parentNodePos, + std::set< string > & siblingNames ); + /** * @brief Recursively read values using ProcessInputFile() from the input * file and put them into the wrapped values for this group. @@ -783,11 +785,11 @@ class Group * but allow to reuse an existing xmlNodePos. * @param[in] xmlDocument the XML document that contains the targetNode. * @param[in] targetNode the XML node that to extract input values from. - * @param[in] nodePos the target node position, typically obtained with xmlDocument::getNodePosition(). + * @param[in] targetNodePos the target node position, typically obtained with xmlDocument::getNodePosition(). */ void processInputFileRecursive( xmlWrapper::xmlDocument & xmlDocument, xmlWrapper::xmlNode & targetNode, - xmlWrapper::xmlNodePos const & nodePos ); + xmlWrapper::xmlNodePos const & targetNodePos ); /** * @brief Recursively call postInputInitialization() to apply post processing after diff --git a/src/coreComponents/dataRepository/GroupContext.cpp b/src/coreComponents/dataRepository/GroupContext.cpp index d21771cfbdb..c60880bd768 100644 --- a/src/coreComponents/dataRepository/GroupContext.cpp +++ b/src/coreComponents/dataRepository/GroupContext.cpp @@ -25,7 +25,7 @@ namespace dataRepository { -GroupContext::GroupContext( Group & group, string const & objectName ): +GroupContext::GroupContext( Group & group, string_view objectName ): DataContext( objectName ), m_group( group ) {} diff --git a/src/coreComponents/dataRepository/GroupContext.hpp b/src/coreComponents/dataRepository/GroupContext.hpp index 476540f98eb..166543b883f 100644 --- a/src/coreComponents/dataRepository/GroupContext.hpp +++ b/src/coreComponents/dataRepository/GroupContext.hpp @@ -57,7 +57,7 @@ class GroupContext : public DataContext * @param group The reference to the Group related to this GroupContext. * @param objectName Target object name. */ - GroupContext( Group & group, string const & objectName ); + GroupContext( Group & group, string_view objectName ); /// The reference to the Group related to this GroupContext. Group & m_group; diff --git a/src/coreComponents/dataRepository/ObjectCatalog.hpp b/src/coreComponents/dataRepository/ObjectCatalog.hpp index dbf089033a6..2fd0762c97b 100644 --- a/src/coreComponents/dataRepository/ObjectCatalog.hpp +++ b/src/coreComponents/dataRepository/ObjectCatalog.hpp @@ -30,6 +30,7 @@ #include "common/logger/Logger.hpp" #include "common/format/StringUtilities.hpp" #include "LvArray/src/system.hpp" +#include "DataContext.hpp" #include #include @@ -171,39 +172,48 @@ class CatalogInterface /** * @brief Static method to create a new object that derives from BASETYPE * @param[in] objectTypeName the key to the catalog entry that is able to create the correct type. + * @param context The data context of the Group for which we attempt to create a sub-group. * @param args these are the arguments to the constructor of the target type * @return passes a unique_ptr to the newly allocated class. - * - * @note The simulation is killed if the builder is not found. + * @note Generate a fatal error: + * - if the object type to create is not found in this Catalog, + * - if the builder is not found. */ //START_SPHINX_2 static std::unique_ptr< BASETYPE > factory( std::string const & objectTypeName, + DataContext const & context, ARGS... args ) { - // We stop the simulation if the product is not found - if( !hasKeyName( objectTypeName ) ) - { - std::list< typename CatalogType::key_type > keys = getKeys(); - string const tmp = stringutilities::join( keys.cbegin(), keys.cend(), ",\n" ); - - string errorMsg = "Could not find keyword \"" + objectTypeName + "\" in this context. "; - errorMsg += "Please be sure that all your keywords are properly spelled or that input file parameters have not changed.\n"; - errorMsg += "All available keys are: [\n" + tmp + "\n]"; - GEOS_ERROR( errorMsg ); - } + // We stop the simulation if the type to create is not found + GEOS_ERROR_IF( !hasKeyName( objectTypeName ), unknownTypeError( objectTypeName, context, getKeys() ) ); // We also stop the simulation if the builder is not here. CatalogInterface< BASETYPE, ARGS... > const * builder = getCatalog().at( objectTypeName ).get(); - if( builder == nullptr ) - { - const string errorMsg = "\"" + objectTypeName + "\" could be found. But the builder is invalid.\n"; - GEOS_ERROR( errorMsg ); - } + GEOS_ERROR_IF( builder == nullptr, + GEOS_FMT( "Type \"{}\" is valid in {}, but the builder is invalid.", + objectTypeName, context ) ); return builder->allocate( args ... ); } //STOP_SPHINX + /** + * @return Generates a formatted error message for an unknown type for a catalog. + * @param objectTypeName The name of the object type that is invalid. + * @param context The data context of the Group for which the erroneous type creation was attempted. + * @param allowedKeys A container of allowed keys, which will be listed in the error message. + * @tparam KEYS_CONTAINER_T A container type holding the allowed keys. + */ + template< typename KEYS_CONTAINER_T > + static string unknownTypeError( string const & objectTypeName, DataContext const & context, + KEYS_CONTAINER_T const & allowedKeys ) + { + return GEOS_FMT( "The tag \"{}\" is invalid within {}. Please verify the keywords spelling and that " + "input file parameters have not changed.\n" + "All available tags are: {}\n", + objectTypeName, context, stringutilities::join( allowedKeys, ", " ) ); + } + /** * @brief Downcast base type reference to derived type * @tparam TYPE type to cast to @@ -303,7 +313,7 @@ class CatalogEntry final : public CatalogInterface< BASETYPE, ARGS... > */ CatalogEntry & operator=( CatalogEntry && source ) { - CatalogInterface< BASETYPE, ARGS... >::operator=( std::move(source)); + CatalogInterface< BASETYPE, ARGS... >::operator=( std::move( source )); } /** @@ -312,7 +322,7 @@ class CatalogEntry final : public CatalogInterface< BASETYPE, ARGS... > * @return a unique_ptr to the newly allocated class. */ //START_SPHINX_4 - virtual std::unique_ptr< BASETYPE > allocate( ARGS... args ) const override + virtual std::unique_ptr< BASETYPE > allocate( ARGS ... args ) const override { #if OBJECTCATALOGVERBOSE > 0 GEOS_LOG( "Creating type " << LvArray::system::demangle( typeid(TYPE).name()) @@ -589,7 +599,7 @@ class CatalogEntry< BASETYPE, TYPE > final : public CatalogInterface< BASETYPE > */ CatalogEntry & operator=( CatalogEntry && source ) { - CatalogInterface< BASETYPE >::operator=( std::move(source)); + CatalogInterface< BASETYPE >::operator=( std::move( source )); } /** diff --git a/src/coreComponents/dataRepository/xmlWrapper.cpp b/src/coreComponents/dataRepository/xmlWrapper.cpp index e0ada302c49..c84439c874f 100644 --- a/src/coreComponents/dataRepository/xmlWrapper.cpp +++ b/src/coreComponents/dataRepository/xmlWrapper.cpp @@ -32,18 +32,29 @@ using namespace dataRepository; namespace xmlWrapper { -void validateString( string const & value, Regex const & regex ) +void validateString( string const & value, const Regex & regex ) { - std::smatch m; - bool inputValidated = std::regex_match( value, m, std::regex( regex.m_regexStr ) ); - if( !inputValidated ) + std::regex regexInstance{ regex.m_regexStr }; + // if validation fails, let's try to find the error start & end to underline it. If it fails, underline the whole message. + if( !std::regex_match( value, regexInstance )) { - ptrdiff_t errorId = ( m.size()>0 && m.position( 0 )==0 ) ? m.length() : 0; + std::smatch m; + std::regex_search( value, m, regexInstance ); + ptrdiff_t const errorStart = ( m.size()>0 && m.position( 0 )==0 ) ? m.length() : 0; + + size_t errorLength; + for( errorLength = 1; errorStart + errorLength < value.length(); ++errorLength ) + { + if( std::regex_match( value.substr( errorStart + errorLength, 1 ), regexInstance ) ) + break; + } + + string const underline = string( errorStart, ' ' ) + string( errorLength, '^' ); GEOS_THROW( GEOS_FMT( "Input string validation failed at:\n" " \"{}\"\n" - " {:>{}}\n" + " {}\n" " Expected format: {}", - value, '^', errorId+1, regex.m_formatDescription ), + value, underline, regex.m_formatDescription ), InputError ); } } diff --git a/src/coreComponents/discretizationMethods/NumericalMethodsManager.cpp b/src/coreComponents/discretizationMethods/NumericalMethodsManager.cpp index 6e4fb3c5cf5..f9dbbb06bc3 100644 --- a/src/coreComponents/discretizationMethods/NumericalMethodsManager.cpp +++ b/src/coreComponents/discretizationMethods/NumericalMethodsManager.cpp @@ -41,6 +41,7 @@ NumericalMethodsManager::~NumericalMethodsManager() Group * NumericalMethodsManager::createChild( string const & GEOS_UNUSED_PARAM( childKey ), string const & GEOS_UNUSED_PARAM( childName ) ) { + // Unused as all children are created within the constructor return nullptr; } diff --git a/src/coreComponents/events/EventBase.cpp b/src/coreComponents/events/EventBase.cpp index 6212dc89d61..52b3c386c2d 100644 --- a/src/coreComponents/events/EventBase.cpp +++ b/src/coreComponents/events/EventBase.cpp @@ -123,7 +123,8 @@ EventBase::CatalogInterface::CatalogType & EventBase::getCatalog() Group * EventBase::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< EventBase > event = EventBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< EventBase > event = + EventBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup< EventBase >( childName, std::move( event ) ); } diff --git a/src/coreComponents/events/EventManager.cpp b/src/coreComponents/events/EventManager.cpp index 69c3a46dd95..3f709bce94b 100644 --- a/src/coreComponents/events/EventManager.cpp +++ b/src/coreComponents/events/EventManager.cpp @@ -94,7 +94,8 @@ EventManager::~EventManager() Group * EventManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< EventBase > event = EventBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< EventBase > event = + EventBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup< EventBase >( childName, std::move( event ) ); } diff --git a/src/coreComponents/events/tasks/TasksManager.cpp b/src/coreComponents/events/tasks/TasksManager.cpp index 9f95aefb2ac..7dfec005c62 100644 --- a/src/coreComponents/events/tasks/TasksManager.cpp +++ b/src/coreComponents/events/tasks/TasksManager.cpp @@ -39,7 +39,8 @@ TasksManager::~TasksManager() Group * TasksManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< TaskBase > task = TaskBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< TaskBase > task = + TaskBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup< TaskBase >( childName, std::move( task ) ); } diff --git a/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp b/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp index 155596304d6..94bbe9ea443 100644 --- a/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp +++ b/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp @@ -52,7 +52,8 @@ FieldSpecificationManager & FieldSpecificationManager::getInstance() Group * FieldSpecificationManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< FieldSpecificationBase > bc = FieldSpecificationBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< FieldSpecificationBase > bc = + FieldSpecificationBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup( childName, std::move( bc ) ); } diff --git a/src/coreComponents/fileIO/Outputs/OutputManager.cpp b/src/coreComponents/fileIO/Outputs/OutputManager.cpp index 42ea71af77a..7fcafb8aa98 100644 --- a/src/coreComponents/fileIO/Outputs/OutputManager.cpp +++ b/src/coreComponents/fileIO/Outputs/OutputManager.cpp @@ -40,7 +40,8 @@ OutputManager::~OutputManager() Group * OutputManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< OutputBase > output = OutputBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< OutputBase > output = + OutputBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup< OutputBase >( childName, std::move( output ) ); } diff --git a/src/coreComponents/finiteElement/FiniteElementDiscretizationManager.cpp b/src/coreComponents/finiteElement/FiniteElementDiscretizationManager.cpp index 8bd43f0cdaf..9bc4e70c5ce 100644 --- a/src/coreComponents/finiteElement/FiniteElementDiscretizationManager.cpp +++ b/src/coreComponents/finiteElement/FiniteElementDiscretizationManager.cpp @@ -41,7 +41,8 @@ Group * FiniteElementDiscretizationManager::createChild( string const & childKey { // These objects should probably not be registered on managed group... GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< Group > fem = Group::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< Group > fem = + Group::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup( childName, std::move( fem ) ); } diff --git a/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp b/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp index 500d1da5ecf..b4329194c67 100644 --- a/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp +++ b/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp @@ -50,7 +50,8 @@ Group * FiniteVolumeManager::createChild( string const & childKey, string const } else { - std::unique_ptr< FluxApproximationBase > approx = FluxApproximationBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< FluxApproximationBase > approx = FluxApproximationBase::CatalogInterface::factory( childKey, getDataContext(), + childName, this ); return &this->registerGroup< FluxApproximationBase >( childName, std::move( approx )); } } diff --git a/src/coreComponents/functions/FunctionManager.cpp b/src/coreComponents/functions/FunctionManager.cpp index 0cd5a92ba9c..2f8425beb99 100644 --- a/src/coreComponents/functions/FunctionManager.cpp +++ b/src/coreComponents/functions/FunctionManager.cpp @@ -54,7 +54,8 @@ Group * FunctionManager::createChild( string const & functionCatalogKey, string const & functionName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), functionCatalogKey, functionName ) ); - std::unique_ptr< FunctionBase > function = FunctionBase::CatalogInterface::factory( functionCatalogKey, functionName, this ); + std::unique_ptr< FunctionBase > function = + FunctionBase::CatalogInterface::factory( functionCatalogKey, getDataContext(), functionName, this ); return &this->registerGroup< FunctionBase >( functionName, std::move( function ) ); } diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp index 652c1d3193d..890f550c752 100644 --- a/src/coreComponents/mainInterface/ProblemManager.cpp +++ b/src/coreComponents/mainInterface/ProblemManager.cpp @@ -159,7 +159,10 @@ ProblemManager::~ProblemManager() Group * ProblemManager::createChild( string const & GEOS_UNUSED_PARAM( childKey ), string const & GEOS_UNUSED_PARAM( childName ) ) -{ return nullptr; } +{ + // Unused as all children are created within the constructor + return nullptr; +} void ProblemManager::problemSetup() @@ -455,52 +458,44 @@ void ProblemManager::parseXMLDocument( xmlWrapper::xmlDocument & xmlDocument ) meshManager.generateMeshLevels( domain ); Group & meshBodies = domain.getMeshBodies(); - // Parse element regions - xmlWrapper::xmlNode elementRegionsNode = xmlProblemNode.child( MeshLevel::groupStructKeys::elemManagerString() ); - - for( xmlWrapper::xmlNode regionNode : elementRegionsNode.children() ) + auto parseRegions = [&]( string_view regionManagerKey, bool const hasParticles ) { - string const regionName = regionNode.attribute( "name" ).value(); - try + xmlWrapper::xmlNode regionsNode = xmlProblemNode.child( regionManagerKey.data() ); + xmlWrapper::xmlNodePos regionsNodePos = xmlDocument.getNodePosition( regionsNode ); + std::set< string > regionNames; + + for( xmlWrapper::xmlNode regionNode : regionsNode.children() ) { - string const - regionMeshBodyName = ElementRegionBase::verifyMeshBodyName( meshBodies, - regionNode.attribute( "meshBody" ).value() ); + auto const regionNodePos = xmlDocument.getNodePosition( regionNode ); + string const regionName = Group::processInputName( regionNode, regionNodePos, + regionsNode.name(), regionsNodePos, regionNames ); + try + { + string const regionMeshBodyName = + ElementRegionBase::verifyMeshBodyName( meshBodies, + regionNode.attribute( "meshBody" ).value() ); - MeshBody & meshBody = domain.getMeshBody( regionMeshBodyName ); - meshBody.forMeshLevels( [&]( MeshLevel & meshLevel ) + MeshBody & meshBody = domain.getMeshBody( regionMeshBodyName ); + meshBody.setHasParticles( hasParticles ); + meshBody.forMeshLevels( [&]( MeshLevel & meshLevel ) + { + ObjectManagerBase & elementManager = hasParticles ? + static_cast< ObjectManagerBase & >( meshLevel.getParticleManager() ): + static_cast< ObjectManagerBase & >( meshLevel.getElemManager() ); + Group * newRegion = elementManager.createChild( regionNode.name(), regionName ); + newRegion->processInputFileRecursive( xmlDocument, regionNode ); + } ); + } + catch( InputError const & e ) { - ElementRegionManager & elementManager = meshLevel.getElemManager(); - Group * newRegion = elementManager.createChild( regionNode.name(), regionName ); - newRegion->processInputFileRecursive( xmlDocument, regionNode ); - } ); - } - catch( InputError const & e ) - { - string const nodePosString = xmlDocument.getNodePosition( regionNode ).toString(); - throw InputError( e, "Error while parsing region " + regionName + " (" + nodePosString + "):\n" ); + throw InputError( e, GEOS_FMT( "Error while parsing region {} ({}):\n", + regionName, regionNodePos.toString() ) ); + } } - } + }; - // Parse particle regions - xmlWrapper::xmlNode particleRegionsNode = xmlProblemNode.child( MeshLevel::groupStructKeys::particleManagerString() ); - - for( xmlWrapper::xmlNode regionNode : particleRegionsNode.children() ) - { - string const regionName = regionNode.attribute( "name" ).value(); - string const - regionMeshBodyName = ElementRegionBase::verifyMeshBodyName( meshBodies, - regionNode.attribute( "meshBody" ).value() ); - - MeshBody & meshBody = domain.getMeshBody( regionMeshBodyName ); - meshBody.setHasParticles( true ); - meshBody.forMeshLevels( [&]( MeshLevel & meshLevel ) - { - ParticleManager & particleManager = meshLevel.getParticleManager(); - Group * newRegion = particleManager.createChild( regionNode.name(), regionName ); - newRegion->processInputFileRecursive( xmlDocument, regionNode ); - } ); - } + parseRegions( MeshLevel::groupStructKeys::elemManagerString(), false ); + parseRegions( MeshLevel::groupStructKeys::particleManagerString(), true ); } } diff --git a/src/coreComponents/mesh/CellElementRegionSelector.cpp b/src/coreComponents/mesh/CellElementRegionSelector.cpp index f5bbc006116..f6106a88a2e 100644 --- a/src/coreComponents/mesh/CellElementRegionSelector.cpp +++ b/src/coreComponents/mesh/CellElementRegionSelector.cpp @@ -67,10 +67,10 @@ CellElementRegionSelector::getMatchingCellblocks( CellElementRegion const & regi "Available cellBlock list: {{ {} }}\nAvailable region attribute list: {{ {} }}", region.getWrapperDataContext( ViewKeys::sourceCellBlockNamesString() ), matchPattern, - stringutilities::joinLamda( m_regionAttributesOwners, ", ", - []( auto pair ) { return pair->first; } ), - stringutilities::joinLamda( m_cellBlocksOwners, ", ", - []( auto pair ) { return pair->first; } ) ), + stringutilities::joinLambda( m_regionAttributesOwners, ", ", + []( auto pair ) { return pair->first; } ), + stringutilities::joinLambda( m_cellBlocksOwners, ", ", + []( auto pair ) { return pair->first; } ) ), InputError ); return matchedCellBlocks; } @@ -86,8 +86,8 @@ CellElementRegionSelector::verifyRequestedCellBlocks( CellElementRegion const & GEOS_FMT( "{}: No cellBlock named '{}'.\nAvailable cellBlock list: {{ {} }}", region.getWrapperDataContext( ViewKeys::sourceCellBlockNamesString() ), requestedCellBlockName, - stringutilities::joinLamda( m_cellBlocksOwners, ", ", - []( auto pair ) { return pair->first; } ) ), + stringutilities::joinLambda( m_cellBlocksOwners, ", ", + []( auto pair ) { return pair->first; } ) ), InputError ); } } @@ -162,7 +162,7 @@ void CellElementRegionSelector::checkSelectionConsistency() const multipleRefsErrors.push_back( GEOS_FMT( "The {} '{}' has been referenced in multiple {}:\n{}", qualifierType, qualifier, CellElementRegion::catalogName(), - stringutilities::joinLamda( owningRegions, '\n', getRegionStr ) ) ); + stringutilities::joinLambda( owningRegions, '\n', getRegionStr ) ) ); } } GEOS_THROW_IF( !multipleRefsErrors.empty(), stringutilities::join( multipleRefsErrors, "\n\n" ), InputError ); diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index 18da04ec066..2652103216e 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -69,31 +69,32 @@ void ElementRegionManager::setMaxGlobalIndex() m_maxGlobalIndex = MpiWrapper::max( m_localMaxGlobalIndex, MPI_COMM_GEOS ); } - +auto const & getUserAvailableKeys() +{ + static std::set< string > keys = { + CellElementRegion::catalogName(), + WellElementRegion::catalogName(), + SurfaceElementRegion::catalogName(), + }; + return keys; +} Group * ElementRegionManager::createChild( string const & childKey, string const & childName ) { - GEOS_ERROR_IF( !(CatalogInterface::hasKeyName( childKey )), - "KeyName ("<getGroup( ElementRegionManager::groupKeyStruct::elementRegionsGroup() ); return &elementRegions.registerGroup( childName, - CatalogInterface::factory( childKey, childName, &elementRegions ) ); + CatalogInterface::factory( childKey, getDataContext(), + childName, &elementRegions ) ); } void ElementRegionManager::expandObjectCatalogs() { - ObjectManagerBase::CatalogInterface::CatalogType const & catalog = ObjectManagerBase::getCatalog(); - for( ObjectManagerBase::CatalogInterface::CatalogType::const_iterator iter = catalog.begin(); - iter!=catalog.end(); - ++iter ) + for( string const & key : getUserAvailableKeys() ) { - string const key = iter->first; - if( key.find( "ElementRegion" ) != string::npos ) - { - this->createChild( key, key ); - } + this->createChild( key, key ); } } diff --git a/src/coreComponents/mesh/ExternalDataSourceBase.cpp b/src/coreComponents/mesh/ExternalDataSourceBase.cpp index 157fe3de31e..9c1285764d7 100644 --- a/src/coreComponents/mesh/ExternalDataSourceBase.cpp +++ b/src/coreComponents/mesh/ExternalDataSourceBase.cpp @@ -28,8 +28,9 @@ ExternalDataSourceBase::ExternalDataSourceBase( string const & name, Group * con Group * ExternalDataSourceBase::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< ExternalDataSourceBase > event = ExternalDataSourceBase::CatalogInterface::factory( childKey, childName, this ); - return &this->registerGroup< ExternalDataSourceBase >( childName, std::move( event ) ); + std::unique_ptr< ExternalDataSourceBase > event = + ExternalDataSourceBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); + return &this->registerGroup( childName, std::move( event ) ); } void ExternalDataSourceBase::expandObjectCatalogs() diff --git a/src/coreComponents/mesh/ExternalDataSourceManager.cpp b/src/coreComponents/mesh/ExternalDataSourceManager.cpp index 8a22203af8e..969ab7a043d 100644 --- a/src/coreComponents/mesh/ExternalDataSourceManager.cpp +++ b/src/coreComponents/mesh/ExternalDataSourceManager.cpp @@ -36,7 +36,8 @@ ExternalDataSourceManager::~ExternalDataSourceManager() Group * ExternalDataSourceManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< ExternalDataSourceBase > externalDataSource = ExternalDataSourceBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< ExternalDataSourceBase > externalDataSource = + ExternalDataSourceBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup< ExternalDataSourceBase >( childName, std::move( externalDataSource ) ); } diff --git a/src/coreComponents/mesh/MeshManager.cpp b/src/coreComponents/mesh/MeshManager.cpp index 662b06002b6..36fdfb0bbfa 100644 --- a/src/coreComponents/mesh/MeshManager.cpp +++ b/src/coreComponents/mesh/MeshManager.cpp @@ -43,7 +43,8 @@ MeshManager::~MeshManager() Group * MeshManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< MeshGeneratorBase > mesh = MeshGeneratorBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< MeshGeneratorBase > mesh = + MeshGeneratorBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup< MeshGeneratorBase >( childName, std::move( mesh ) ); } diff --git a/src/coreComponents/mesh/ParticleManager.cpp b/src/coreComponents/mesh/ParticleManager.cpp index 314c23120b4..057533e4f7e 100644 --- a/src/coreComponents/mesh/ParticleManager.cpp +++ b/src/coreComponents/mesh/ParticleManager.cpp @@ -86,14 +86,11 @@ void ParticleManager::setMaxGlobalIndex() Group * ParticleManager::createChild( string const & childKey, string const & childName ) { - GEOS_ERROR_IF( !(CatalogInterface::hasKeyName( childKey )), - "KeyName ("<getGroup( ParticleManager::groupKeyStruct::particleRegionsGroup() ); return &particleRegions.registerGroup( childName, - CatalogInterface::factory( childKey, childName, &particleRegions ) ); - + CatalogInterface::factory( childKey, getDataContext(), + childName, &particleRegions ) ); } void ParticleManager::expandObjectCatalogs() diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index 2dd41c0987a..8ecabcba9ea 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -46,6 +46,7 @@ void CellBlockManager::resize( integer_array const & numElements, Group * CellBlockManager::createChild( string const & GEOS_UNUSED_PARAM( childKey ), string const & GEOS_UNUSED_PARAM( childName ) ) { + // Unused as all children are created within the constructor return nullptr; } diff --git a/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp b/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp index 555bef92167..86eb7f0b849 100644 --- a/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp +++ b/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp @@ -30,8 +30,9 @@ MeshGeneratorBase::MeshGeneratorBase( string const & name, Group * const parent Group * MeshGeneratorBase::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< MeshComponentBase > comp = MeshComponentBase::CatalogInterface::factory( childKey, childName, this ); - return &this->registerGroup< MeshComponentBase >( childName, std::move( comp ) ); + std::unique_ptr< MeshComponentBase > meshComp = + MeshComponentBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); + return &this->registerGroup< MeshComponentBase >( childName, std::move( meshComp ) ); } void MeshGeneratorBase::expandObjectCatalogs() diff --git a/src/coreComponents/mesh/generators/ParticleBlockManager.cpp b/src/coreComponents/mesh/generators/ParticleBlockManager.cpp index d96f11711dc..f219072ab8e 100644 --- a/src/coreComponents/mesh/generators/ParticleBlockManager.cpp +++ b/src/coreComponents/mesh/generators/ParticleBlockManager.cpp @@ -38,6 +38,7 @@ void ParticleBlockManager::resize( integer_array const & numParticles, Group * ParticleBlockManager::createChild( string const & GEOS_UNUSED_PARAM( childKey ), string const & GEOS_UNUSED_PARAM( childName ) ) { + // Unused as all children are created within the constructor return nullptr; } diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp index f64e0e9b7d0..af80ec09c8a 100644 --- a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp +++ b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp @@ -74,20 +74,15 @@ WellGeneratorBase::WellGeneratorBase( string const & name, Group * const parent Group * WellGeneratorBase::createChild( string const & childKey, string const & childName ) { - if( childKey == viewKeyStruct::perforationString() ) - { - ++m_numPerforations; + GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); + const auto childTypes = { viewKeyStruct::perforationString() }; + GEOS_ERROR_IF( childKey != viewKeyStruct::perforationString(), + CatalogInterface::unknownTypeError( childKey, getDataContext(), childTypes ) ); - // keep track of the perforations that have been added - m_perforationList.emplace_back( childName ); - GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - return ®isterGroup< Perforation >( childName ); - } - else - { - GEOS_THROW( "Unrecognized node: " << childKey, InputError ); - } - return nullptr; + ++m_numPerforations; + m_perforationList.emplace_back( childName ); + + return ®isterGroup< Perforation >( childName ); } void WellGeneratorBase::expandObjectCatalogs() diff --git a/src/coreComponents/mesh/simpleGeometricObjects/GeometricObjectManager.cpp b/src/coreComponents/mesh/simpleGeometricObjects/GeometricObjectManager.cpp index c5843dbc747..4209636f4d5 100644 --- a/src/coreComponents/mesh/simpleGeometricObjects/GeometricObjectManager.cpp +++ b/src/coreComponents/mesh/simpleGeometricObjects/GeometricObjectManager.cpp @@ -54,7 +54,8 @@ GeometricObjectManager & GeometricObjectManager::getInstance() Group * GeometricObjectManager::createChild( string const & childKey, string const & childName ) { GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - std::unique_ptr< SimpleGeometricObjectBase > geometriObject = SimpleGeometricObjectBase::CatalogInterface::factory( childKey, childName, this ); + std::unique_ptr< SimpleGeometricObjectBase > geometriObject = + SimpleGeometricObjectBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); return &this->registerGroup< SimpleGeometricObjectBase >( childName, std::move( geometriObject ) ); } diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 669f721b5f4..e0eefbde8c2 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -191,6 +191,7 @@ void PhysicsSolverBase::registerDataOnMesh( Group & meshBodies ) Group * PhysicsSolverBase::createChild( string const & GEOS_UNUSED_PARAM( childKey ), string const & GEOS_UNUSED_PARAM( childName ) ) { + // Unused as all children are created within the constructor return nullptr; } diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverManager.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverManager.cpp index 0267c7df811..9301cf5ee67 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverManager.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverManager.cpp @@ -46,14 +46,10 @@ PhysicsSolverManager::~PhysicsSolverManager() //START_SPHINX_INCLUDE_00 Group * PhysicsSolverManager::createChild( string const & childKey, string const & childName ) { - Group * rval = nullptr; - if( PhysicsSolverBase::CatalogInterface::hasKeyName( childKey ) ) - { - GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); - rval = ®isterGroup( childName, - PhysicsSolverBase::CatalogInterface::factory( childKey, childName, this ) ); - } - return rval; + GEOS_LOG_RANK_0( GEOS_FMT( "{}: adding {} {}", getName(), childKey, childName ) ); + std::unique_ptr< PhysicsSolverBase > solver = + PhysicsSolverBase::CatalogInterface::factory( childKey, getDataContext(), childName, this ); + return ®isterGroup< PhysicsSolverBase >( childName, std::move( solver ) ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 5bb8896898d..c13be912619 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -42,8 +42,6 @@ #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SourceFluxStatistics.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/PhaseVolumeFractionKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseVolumeFractionKernel.hpp" @@ -418,7 +416,7 @@ void CompositionalMultiphaseBase::setConstitutiveNames( ElementSubRegionBase & s string & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); fluidName = getConstitutiveName< MultiFluidBase >( subRegion ); GEOS_THROW_IF( fluidName.empty(), - GEOS_FMT( "{}: Fluid model not found on subregion {}", + GEOS_FMT( "{}: multiphase fluid model not found on subregion {}", getDataContext(), subRegion.getDataContext() ), InputError ); @@ -1333,7 +1331,6 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom if( m_useSimpleAccumulation ) kernelFlags.set( KernelFlags::SimpleAccumulation ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, arrayView1d< string const > const & regionNames ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 318662d00fc..cb9c1fe33d3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -22,6 +22,10 @@ #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/solid/CoupledSolidBase.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp" namespace geos { @@ -98,6 +102,21 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override; + /** + * @brief assembles the accumulation terms for all cells of a spcefici subRegion. + * @tparam SUBREGION_TYPE the subRegion type + * @param dofManager degree-of-freedom manager associated with the linear system + * @param subRegion the subRegion + * @param localMatrix the system matrix + * @param localRhs the system right-hand side vector + * + */ + template< typename SUBREGION_TYPE > + void accumulationAssemblyLaunch( DofManager const & dofManager, + SUBREGION_TYPE const & subRegion, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + virtual void resetStateToBeginningOfStep( DomainPartition & domain ) override; @@ -541,6 +560,58 @@ void CompositionalMultiphaseBase::applyFieldValue( real64 const & time_n, } ); } +template< typename SUBREGION_TYPE > +void CompositionalMultiphaseBase::accumulationAssemblyLaunch( DofManager const & dofManager, + SUBREGION_TYPE const & subRegion, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + constitutive::MultiFluidBase const & fluid = + getConstitutiveModel< constitutive::MultiFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); + constitutive::CoupledSolidBase const & solid = + getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + + string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + + using namespace isothermalCompositionalMultiphaseBaseKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + if( m_useSimpleAccumulation ) + kernelFlags.set( KernelFlags::SimpleAccumulation ); + + if( m_isThermal ) + { + thermalCompositionalMultiphaseBaseKernels:: + AccumulationKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + kernelFlags, + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } + else + { + isothermalCompositionalMultiphaseBaseKernels:: + AccumulationKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + kernelFlags, + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } +} } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index d2dbcb2f3a0..aeb88ef68e3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -54,12 +54,14 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseMobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/AquiferBCKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp" namespace geos { using namespace dataRepository; using namespace constitutive; +using namespace compositionalMultiphaseUtilities; // for ScalingType CompositionalMultiphaseFVM::CompositionalMultiphaseFVM( const string & name, Group * const parent ) @@ -1264,9 +1266,89 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, } -real64 CompositionalMultiphaseFVM::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain ) +void CompositionalMultiphaseFVM::assembleHydrofracFluxTerms( real64 const GEOS_UNUSED_PARAM ( time_n ), + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) { + GEOS_MARK_FUNCTION; + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); + + string const & elemDofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + + using namespace isothermalCompositionalMultiphaseFVMKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_hasCapPressure ) + kernelFlags.set( KernelFlags::CapPressure ); + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + if( m_gravityDensityScheme == GravityDensityScheme::PhasePresence ) + kernelFlags.set( KernelFlags::CheckPhasePresenceInGravity ); + if( m_isThermal ) + { + // should not end up here but just in case + GEOS_ERROR( "Thermal not yet supported in CompositionalMultiphaseFVM::assembleHydrofracFluxTerms" ); + } + if( fluxApprox.upwindingParams().upwindingScheme != UpwindingScheme::PPU ) + { + // a bit tricky to check in advance + GEOS_ERROR( "Only PPU upwinding is supported in CompositionalMultiphaseFVM::assembleHydrofracFluxTerms" ); + } + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & ) + { + fluxApprox.forStencils< CellElementStencilTPFA, FaceElementToCellStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + isothermalCompositionalMultiphaseFVMKernels:: + FluxComputeKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + kernelFlags, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + multiphasePoromechanicsConformingFracturesKernels:: + FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + kernelFlags, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView(), + dR_dAper ); + } ); + } ); +} + +real64 CompositionalMultiphaseFVM::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain ) +{ if( m_targetFlowCFL < 0 ) return CompositionalMultiphaseBase::setNextDt( currentDt, domain ); else diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index 896dc87b7f1..14d66667035 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -163,6 +163,14 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const override; + virtual void + assembleHydrofracFluxTerms( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; @@ -205,15 +213,6 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase static constexpr char const * gravityDensitySchemeString() { return "gravityDensityScheme"; } }; - /** - * @brief Solution scaling type - */ - enum class ScalingType : integer - { - Global, ///< Scale the Newton update with a unique scaling factor - Local ///< Scale the Newton update locally (modifies the Newton direction) - }; - /** * @brief Storage for value and element location, used to determine global max + location */ @@ -255,7 +254,7 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase } m_dbcParams; /// Solution scaling type - ScalingType m_scalingType; + compositionalMultiphaseUtilities::ScalingType m_scalingType; /// scheme for density treatment in gravity GravityDensityScheme m_gravityDensityScheme; @@ -293,10 +292,6 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase }; -ENUM_STRINGS( CompositionalMultiphaseFVM::ScalingType, - "Global", - "Local" ); - } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index 021ce8fd96b..53415a507cd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -541,7 +541,7 @@ bool CompositionalMultiphaseHybridFVM::checkSystemSolution( DomainPartition & do SolutionCheckKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping, m_allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType::Global, + compositionalMultiphaseUtilities::ScalingType::Global, scalingFactor, pressure, compDens, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp index 3af76d4928d..cbac0f13f2d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp @@ -28,6 +28,19 @@ namespace geos namespace compositionalMultiphaseUtilities { +/** + * @brief Solution scaling type, used in CompositionalMultiphaseFVM + */ +enum class ScalingType : integer +{ + Global, ///< Scale the Newton update with a unique scaling factor + Local ///< Scale the Newton update locally (modifies the Newton direction) +}; + +ENUM_STRINGS( ScalingType, + "Global", + "Local" ); + /** * @brief In each block, shift the elements from 0 to numRowsToShift-1 one position ahead and replaces the first element * with the sum of all elements from 0 to numRowsToShift-1 of the input block one-dimensional array of values. diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index c081f20d4ce..e9a55b86c46 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -143,6 +143,7 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies ) [&]( localIndex const, ElementSubRegionBase & subRegion ) { + subRegion.registerField< fields::flow::deltaVolume >( getName() ); subRegion.registerField< fields::flow::gravityCoefficient >( getName() ). setApplyDefaultValue( 0.0 ); subRegion.registerField< fields::flow::netToGross >( getName() ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 4352e059405..b97716d0b00 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -146,6 +146,28 @@ class FlowSolverBase : public PhysicsSolverBase real64 const & timeAtBeginningOfStep, real64 const & dt ); + /** + * @brief assembles the flux terms for all cells for the hydrofracture case + * @param time_n previous time value + * @param dt time step + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix the system matrix + * @param localRhs the system right-hand side vector + * @param dR_dAper + */ + virtual void assembleHydrofracFluxTerms( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + { + GEOS_UNUSED_VAR ( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper ); + GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." ); + } + virtual void initializeFluidState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } virtual void initializeThermalState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 98934988d8d..94f3d46a640 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -93,15 +93,12 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies ) [&]( localIndex const, ElementSubRegionBase & subRegion ) { - subRegion.registerField< deltaVolume >( getName() ); - subRegion.registerField< mobility >( getName() ); subRegion.registerField< dMobility_dPressure >( getName() ); subRegion.registerField< fields::flow::mass >( getName() ); subRegion.registerField< fields::flow::mass_n >( getName() ); - if( m_isThermal ) { subRegion.registerField< dMobility_dTemperature >( getName() ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index db72beaa75a..a4049c28a87 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -205,25 +205,6 @@ class SinglePhaseBase : public FlowSolverBase arrayView1d< real64 > const & localRhs, string const & jumpDofKey ) = 0; - /** - * @brief assembles the flux terms for all cells for the hydrofracture case - * @param time_n previous time value - * @param dt time step - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix the system matrix - * @param localRhs the system right-hand side vector - * @param dR_dAper - */ - virtual void - assembleHydrofracFluxTerms( real64 const time_n, - real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) = 0; - struct viewKeyStruct : FlowSolverBase::viewKeyStruct { static constexpr char const * elemDofFieldString() { return "singlePhaseVariables"; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp index e0b9732cec2..5884b6fc416 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp @@ -311,25 +311,6 @@ void SinglePhaseHybridFVM::assembleEDFMFluxTerms( real64 const GEOS_UNUSED_PARAM localRhs ); } -void SinglePhaseHybridFVM::assembleHydrofracFluxTerms( real64 const time_n, - real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) -{ - GEOS_UNUSED_VAR ( time_n ); - GEOS_UNUSED_VAR ( dt ); - GEOS_UNUSED_VAR ( domain ); - GEOS_UNUSED_VAR ( dofManager ); - GEOS_UNUSED_VAR ( localMatrix ); - GEOS_UNUSED_VAR ( localRhs ); - GEOS_UNUSED_VAR ( dR_dAper ); - - GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." ); -} - void SinglePhaseHybridFVM::applyBoundaryConditions( real64 const time_n, real64 const dt, DomainPartition & domain, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp index 7a804ccf7ee..b2910490508 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp @@ -142,15 +142,6 @@ class SinglePhaseHybridFVM : public SinglePhaseBase arrayView1d< real64 > const & localRhs, string const & jumpDofKey ) override final; - virtual void - assembleHydrofracFluxTerms( real64 const time_n, - real64 const dt, - DomainPartition const & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs, - CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; - virtual void applyAquiferBC( real64 const time, real64 const dt, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp index f92335dd726..3d79301623a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp @@ -107,6 +107,7 @@ struct C1PPUPhaseFlux real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp], real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] ) { + real64 dPotGrad_dTrans = 0.0; real64 dPresGrad_dP[numFluxSupportPoints]{}; real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; @@ -114,7 +115,8 @@ struct C1PPUPhaseFlux PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, checkPhasePresenceInGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, - phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + potGrad, dPotGrad_dTrans, dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); // gravity head diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp index 78271f54416..45882953307 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp @@ -252,6 +252,7 @@ class FluxComputeKernel : public FluxComputeKernelBase real64 compFlux[numComp]{}; real64 dCompFlux_dP[numFluxSupportPoints][numComp]{}; real64 dCompFlux_dC[numFluxSupportPoints][numComp][numComp]{}; + real64 dCompFlux_dTrans[numComp]{}; real64 const trans[numFluxSupportPoints] = { stack.transmissibility[connectionIndex][0], stack.transmissibility[connectionIndex][1] }; @@ -350,7 +351,8 @@ class FluxComputeKernel : public FluxComputeKernelBase dPhaseFlux_dC, compFlux, dCompFlux_dP, - dCompFlux_dC ); + dCompFlux_dC, + dCompFlux_dTrans ); } // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp index 9365a96c830..97f5bdd5a2e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp @@ -1760,7 +1760,7 @@ struct IHUPhaseFlux real64 dummy[numComp]; real64 dDummy_dP[numFluxSupportPoints][numComp]; real64 dDummy_dC[numFluxSupportPoints][numComp][numComp]; - + real64 dDummy_dTrans[numComp]; for( integer jp = 0; jp < numPhase; ++jp ) { @@ -1777,7 +1777,7 @@ struct IHUPhaseFlux phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, k_up_ppu, potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, - dummy, dDummy_dP, dDummy_dC ); + dummy, dDummy_dP, dDummy_dC, dDummy_dTrans ); totFlux += phaseFlux; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp index 468adabb4b1..14fc73f029d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp @@ -68,6 +68,10 @@ struct PPUPhaseFlux * @param phaseFlux phase flux * @param dPhaseFlux_dP derivative of phase flux wrt pressure * @param dPhaseFlux_dC derivative of phase flux wrt comp density + * @param compFlux component flux + * @param dCompFlux_dP derivative of phase flux wrt pressure + * @param dCompFlux_dC derivative of phase flux wrt comp density + * @param dCompFlux_dTrans derivative of phase flux wrt transmissibility */ template< integer numComp, integer numFluxSupportPoints > GEOS_HOST_DEVICE @@ -96,13 +100,15 @@ struct PPUPhaseFlux ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, localIndex & k_up, real64 & potGrad, - real64 ( &phaseFlux ), + real64 & phaseFlux, real64 ( & dPhaseFlux_dP )[numFluxSupportPoints], real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp], real64 ( & compFlux )[numComp], real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp], - real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] ) + real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp], + real64 ( & dCompFlux_dTrans )[numComp] ) { + real64 dPotGrad_dTrans = 0; real64 dPresGrad_dP[numFluxSupportPoints]{}; real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; @@ -111,8 +117,9 @@ struct PPUPhaseFlux seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, - phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, - dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + potGrad, dPotGrad_dTrans, dPresGrad_dP, + dPresGrad_dC, dGravHead_dP, dGravHead_dC ); // *** upwinding *** @@ -125,6 +132,11 @@ struct PPUPhaseFlux real64 const mobility = phaseMob[er_up][esr_up][ei_up][ip]; + // compute phase flux using upwind mobility + phaseFlux = mobility * potGrad; + + real64 const dPhaseFlux_dTrans = mobility * dPotGrad_dTrans; + // pressure gradient depends on all points in the stencil for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { @@ -136,8 +148,6 @@ struct PPUPhaseFlux dPhaseFlux_dC[ke][jc] *= mobility; } } - // compute phase flux using upwind mobility. - phaseFlux = mobility * potGrad; real64 const dMob_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP]; arraySlice1d< real64 const, compflow::USD_PHASE_DC - 2 > dPhaseMobSub = @@ -151,8 +161,8 @@ struct PPUPhaseFlux } //distribute on phaseComponentFlux here - PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, phaseFlux - , dPhaseFlux_dP, dPhaseFlux_dC, compFlux, dCompFlux_dP, dCompFlux_dC ); + PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans, compFlux, dCompFlux_dP, dCompFlux_dC, dCompFlux_dTrans ); } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp index 3d7f1cba38b..f7810873e5e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp @@ -119,6 +119,94 @@ struct PhaseComponentFlux } } } + + /** + * @brief Compute the component flux for a given phase + * @tparam numComp number of components + * @tparam numFluxSupportPoints number of flux support points + * @param ip phase index + * @param k_up uptream index for this phase + * @param seri arraySlice of the stencil-implied element region index + * @param sesri arraySlice of the stencil-implied element subregion index + * @param sei arraySlice of the stencil-implied element index + * @param phaseCompFrac phase component fraction + * @param dPhaseCompFrac derivative of phase component fraction wrt pressure, temperature, component fraction + * @param dCompFrac_dCompDens derivative of component fraction wrt component density + * @param phaseFlux phase flux + * @param dPhaseFlux_dP derivative of phase flux wrt pressure + * @param dPhaseFlux_dC derivative of phase flux wrt comp density + * @param dPhaseFlux_dTrans derivative of phase flux wrt transmissibility + * @param compFlux component flux + * @param dCompFlux_dP derivative of phase flux wrt pressure + * @param dCompFlux_dC derivative of phase flux wrt comp density + * @param dCompFlux_dTrans derivative of phase flux wrt transmissibility + */ + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + compute( localIndex const ip, + localIndex const k_up, + localIndex const ( &seri )[numFluxSupportPoints], + localIndex const ( &sesri )[numFluxSupportPoints], + localIndex const ( &sei )[numFluxSupportPoints], + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, + ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, + ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, + real64 const & phaseFlux, + real64 const ( &dPhaseFlux_dP )[numFluxSupportPoints], + real64 const ( &dPhaseFlux_dC )[numFluxSupportPoints][numComp], + real64 const & dPhaseFlux_dTrans, + real64 ( & compFlux )[numComp], + real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp], + real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp], + real64 ( & dCompFlux_dTrans )[numComp] ) + { + localIndex const er_up = seri[k_up]; + localIndex const esr_up = sesri[k_up]; + localIndex const ei_up = sei[k_up]; + + real64 dProp_dC[numComp]{}; + + // slice some constitutive arrays to avoid too much indexing in component loop + arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE_COMP-3 > phaseCompFracSub = + phaseCompFrac[er_up][esr_up][ei_up][0][ip]; + arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC-3 > dPhaseCompFracSub = + dPhaseCompFrac[er_up][esr_up][ei_up][0][ip]; + + // compute component fluxes and derivatives using upstream cell composition + for( integer ic = 0; ic < numComp; ++ic ) + { + real64 const ycp = phaseCompFracSub[ic]; + compFlux[ic] += phaseFlux * ycp; + + // derivatives stemming from phase flux + for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) + { + dCompFlux_dP[ke][ic] += dPhaseFlux_dP[ke] * ycp; + for( integer jc = 0; jc < numComp; ++jc ) + { + dCompFlux_dC[ke][ic][jc] += dPhaseFlux_dC[ke][jc] * ycp; + } + } + + // additional derivatives stemming from upstream cell phase composition + dCompFlux_dP[k_up][ic] += phaseFlux * dPhaseCompFracSub[ic][Deriv::dP]; + + // convert derivatives of comp fraction w.r.t. comp fractions to derivatives w.r.t. comp densities + applyChainRule( numComp, + dCompFrac_dCompDens[er_up][esr_up][ei_up], + dPhaseCompFracSub[ic], + dProp_dC, + Deriv::dC ); + for( integer jc = 0; jc < numComp; ++jc ) + { + dCompFlux_dC[k_up][ic][jc] += phaseFlux * dProp_dC[jc]; + } + + dCompFlux_dTrans[ic] += dPhaseFlux_dTrans * ycp; + } + } + }; } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp index b783b061712..f9ca248b4dc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp @@ -62,6 +62,7 @@ struct PotGrad ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, real64 & potGrad, + real64 & dPotGrad_dTrans, real64 ( & dPresGrad_dP )[numFluxSupportPoints], real64 ( & dPresGrad_dC )[numFluxSupportPoints][numComp], real64 ( & dGravHead_dP )[numFluxSupportPoints], @@ -85,7 +86,9 @@ struct PotGrad real64 dDensMean_dC[numFluxSupportPoints][numComp]{}; real64 presGrad = 0.0; + real64 dPresGrad_dTrans = 0.0; real64 gravHead = 0.0; + real64 dGravHead_dTrans = 0.0; real64 dCapPressure_dC[numComp]{}; calculateMeanDensity( checkPhasePresenceInGravity, ip, @@ -126,22 +129,25 @@ struct PotGrad } } - presGrad += trans[i] * (pres[er][esr][ei] - capPressure); - dPresGrad_dP[i] += trans[i] * (1 - dCapPressure_dP) - + dTrans_dPres[i] * (pres[er][esr][ei] - capPressure); + real64 const dP = pres[er][esr][ei] - capPressure; + presGrad += trans[i] * dP; + dPresGrad_dTrans += dP; + dPresGrad_dP[i] += trans[i] * (1 - dCapPressure_dP) + dTrans_dPres[i] * dP; for( integer jc = 0; jc < numComp; ++jc ) { dPresGrad_dC[i][jc] += -trans[i] * dCapPressure_dC[jc]; } - real64 const gravD = trans[i] * gravCoef[er][esr][ei]; - real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei]; + real64 const gC = gravCoef[er][esr][ei]; + real64 const gravD = trans[i] * gC; + real64 const dGravD_dTrans = gC; + real64 const dGravD_dP = dTrans_dPres[i] * gC; // the density used in the potential difference is always a mass density // unlike the density used in the phase mobility, which is a mass density // if useMass == 1 and a molar density otherwise gravHead += densMean * gravD; - + dGravHead_dTrans += densMean * dGravD_dTrans; // need to add contributions from both cells the mean density depends on for( integer j = 0; j < numFluxSupportPoints; ++j ) { @@ -155,7 +161,7 @@ struct PotGrad // compute phase potential gradient potGrad = presGrad - gravHead; - + dPotGrad_dTrans = dPresGrad_dTrans - dGravHead_dTrans; } template< integer numComp, integer numFluxSupportPoints > diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp index 5cbf59b932d..89d681bb02c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp @@ -61,7 +61,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer */ SolutionCheckKernel( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView2d< real64 const, compflow::USD_COMP > const compDens, @@ -217,7 +217,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer StackVariables & stack, FUNC && kernelOp = NoOpFunc{} ) const { - bool const localScaling = m_scalingType == CompositionalMultiphaseFVM::ScalingType::Local; + bool const localScaling = m_scalingType == compositionalMultiphaseUtilities::ScalingType::Local; real64 const newPres = m_pressure[ei] + (localScaling ? m_pressureScalingFactor[ei] : m_scalingFactor) * m_localSolution[stack.localRow]; if( newPres < 0 ) @@ -280,7 +280,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer real64 const m_scalingFactor; /// scaling type (global or local) - CompositionalMultiphaseFVM::ScalingType const m_scalingType; + compositionalMultiphaseUtilities::ScalingType const m_scalingType; }; @@ -306,7 +306,7 @@ class SolutionCheckKernelFactory static SolutionCheckKernel::StackVariables createAndLaunch( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView2d< real64 const, compflow::USD_COMP > const compDens, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp index 35c969bb891..b0447b33782 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp @@ -60,7 +60,7 @@ class SolutionCheckKernel : public isothermalCompositionalMultiphaseBaseKernels: */ SolutionCheckKernel( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView1d< real64 const > const temperature, @@ -103,7 +103,7 @@ class SolutionCheckKernel : public isothermalCompositionalMultiphaseBaseKernels: { Base::computeSolutionCheck( ei, stack, [&] () { - bool const localScaling = m_scalingType == CompositionalMultiphaseFVM::ScalingType::Local; + bool const localScaling = m_scalingType == compositionalMultiphaseUtilities::ScalingType::Local; // compute the change in temperature real64 const newTemp = m_temperature[ei] + (localScaling ? m_temperatureScalingFactor[ei] : m_scalingFactor * m_localSolution[stack.localRow + m_temperatureOffset]); if( newTemp < minTemperature ) @@ -149,7 +149,7 @@ class SolutionCheckKernelFactory static SolutionCheckKernel::StackVariables createAndLaunch( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView1d< real64 const > const temperature, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index b891467863d..e33e4ec0b17 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -1535,7 +1535,7 @@ CompositionalMultiphaseWell::checkSystemSolution( DomainPartition & domain, compositionalMultiphaseWellKernels:: SolutionCheckKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping, - CompositionalMultiphaseFVM::ScalingType::Global, + compositionalMultiphaseUtilities::ScalingType::Global, scalingFactor, pressure, compDens, @@ -1573,7 +1573,7 @@ CompositionalMultiphaseWell::checkSystemSolution( DomainPartition & domain, { //integer const m_allowCompDensChopping(true); integer const m_allowNegativePressure( false ); - CompositionalMultiphaseFVM::ScalingType const m_scalingType( CompositionalMultiphaseFVM::ScalingType::Global ); + compositionalMultiphaseUtilities::ScalingType const m_scalingType( compositionalMultiphaseUtilities::ScalingType::Global ); arrayView1d< real64 const > const pressure = subRegion.getField< fields::well::pressure >(); arrayView1d< real64 const > const temperature = diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp index e7496035ec6..1b7ed01bb4e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp @@ -67,17 +67,10 @@ WellSolverBase::WellSolverBase( string const & name, Group *WellSolverBase::createChild( string const & childKey, string const & childName ) { - Group *rval = nullptr; - - if( childKey == keys::wellControls ) - { - rval = ®isterGroup< WellControls >( childName ); - } - else - { - PhysicsSolverBase::createChild( childKey, childName ); - } - return rval; + const auto childTypes = { keys::wellControls }; + GEOS_ERROR_IF( childKey != keys::wellControls, + PhysicsSolverBase::CatalogInterface::unknownTypeError( childKey, getDataContext(), childTypes ) ); + return ®isterGroup< WellControls >( childName ); } void WellSolverBase::expandObjectCatalogs() diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index c2e9078ba47..5722613372d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -810,7 +810,7 @@ class SolutionCheckKernelFactory template< typename POLICY > static isothermalCompositionalMultiphaseBaseKernels::SolutionCheckKernel::StackVariables createAndLaunch( integer const allowCompDensChopping, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView2d< real64 const, compflow::USD_COMP > const compDens, diff --git a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt index 89a34cbc276..2332116db38 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt @@ -9,6 +9,7 @@ set( physicsSolvers_headers multiphysics/HydrofractureSolver.hpp multiphysics/HydrofractureSolverKernels.hpp multiphysics/MultiphasePoromechanics.hpp + multiphysics/MultiphasePoromechanicsConformingFractures.hpp multiphysics/PhaseFieldFractureSolver.hpp multiphysics/PoromechanicsInitialization.hpp multiphysics/PoromechanicsFields.hpp @@ -16,6 +17,7 @@ set( physicsSolvers_headers multiphysics/LogLevelsInfo.hpp multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp + multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp multiphysics/poromechanicsKernels/PoromechanicsBase.hpp multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp @@ -46,6 +48,7 @@ set( physicsSolvers_sources multiphysics/FlowProppantTransportSolver.cpp multiphysics/HydrofractureSolver.cpp multiphysics/MultiphasePoromechanics.cpp + multiphysics/MultiphasePoromechanicsConformingFractures.cpp multiphysics/PhaseFieldFractureSolver.cpp multiphysics/PoromechanicsInitialization.cpp multiphysics/SinglePhasePoromechanics.cpp diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp index 6b01f872f5c..0d11624cfc1 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp @@ -32,7 +32,7 @@ #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp" -//#include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" +#include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" //#include "physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp" namespace geos @@ -361,7 +361,7 @@ void MultiphasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::updateBulkDensity } template class MultiphasePoromechanics<>; -//template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsLagrangeContact >; +template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsLagrangeContact >; //template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsEmbeddedFractures >; template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> >; diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp new file mode 100644 index 00000000000..8ef737ad998 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp @@ -0,0 +1,833 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultiphasePoromechanicsConformingFractures.cpp + */ + +#include "MultiphasePoromechanicsConformingFractures.hpp" + +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" +#include "finiteVolume/FluxApproximationBase.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; +using namespace fields; + +template< typename FLOW_SOLVER > +MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::MultiphasePoromechanicsConformingFractures( const string & name, + Group * const parent ) + : Base( name, parent ) +{ + // TODO: MGR recipe + // LinearSolverParameters & params = this->m_linearSolverParameters.get(); + // params.mgr.strategy = LinearSolverParameters::MGR::StrategyType::multiphasePoromechanicsConformingFractures; + // params.mgr.separateComponents = false; + // params.mgr.displacementFieldName = solidMechanics::totalDisplacement::key(); + // params.dofsPerNode = 3; +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::postInputInitialization() +{ + Base::postInputInitialization(); + + if( this->flowSolver()->isThermal() ) + { + GEOS_ERROR( "Thermal flow is not yet supported for multiphase poromechanics conforming fractures" ); + } +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const +{ + /// We need to add 2 coupling terms: + // 1. Poromechanical coupling in the bulk + Base::setupCoupling( domain, dofManager ); + + // 2. Traction - pressure coupling in the fracture + dofManager.addCoupling( m_flowDofKey, + fields::contact::traction::key(), + DofManager::Connector::Elem ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity ) +{ + GEOS_MARK_FUNCTION; + + GEOS_UNUSED_VAR( setSparsity ); + + /// 1. Add all coupling terms handled directly by the DofManager + dofManager.setDomain( domain ); + this->setupDofs( domain, dofManager ); + dofManager.reorderByRank(); + + /// 2. Add coupling terms not added by the DofManager. + localIndex const numLocalRows = dofManager.numLocalDofs(); + + SparsityPattern< globalIndex > patternOriginal; + dofManager.setSparsityPattern( patternOriginal ); + + // Get the original row lengths (diagonal blocks only) + array1d< localIndex > rowLengths( patternOriginal.numRows() ); + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); + } + + // Add the number of nonzeros induced by coupling + addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); + + // Create a new pattern with enough capacity for coupled matrix + SparsityPattern< globalIndex > pattern; + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), + patternOriginal.numColumns(), + rowLengths.data() ); + + // Copy the original nonzeros + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + } + + // Add the nonzeros from coupling + addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); + + localMatrix.setName( this->getName() + "/matrix" ); + localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); + + rhs.setName( this->getName() + "/rhs" ); + rhs.create( numLocalRows, MPI_COMM_GEOS ); + + solution.setName( this->getName() + "/solution" ); + solution.create( numLocalRows, MPI_COMM_GEOS ); + + setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); + + // if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + // { + // createPreconditioner( domain ); + // } +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::assembleSystem( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + this->solidMechanicsSolver()->synchronizeFractureState( domain ); + + assembleElementBasedContributions( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs ); + + // Assemble fluxes 3D/2D and get dFluidResidualDAperture + this->flowSolver()->assembleHydrofracFluxTerms( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs, + getDerivativeFluxResidual_dNormalJump() ); + + // This step must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. + assembleCouplingTerms( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::assembleElementBasedContributions( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_UNUSED_VAR( time_n, dt ); + + /// 3. assemble Force Residual w.r.t. pressure and Flow mass residual w.r.t. displacement + + Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); + + // Flow accumulation for fractures + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs ); + } ); + } ); + + this->solidMechanicsSolver()->assembleContact( domain, dofManager, localMatrix, localRhs ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::assembleCouplingTerms( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_UNUSED_VAR( time_n, dt ); + // These 2 steps need to occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & regionNames ) + { + /// 3. assemble Force Residual w.r.t. pressure and Flow mass residual w.r.t. displacement + assembleForceResidualDerivativeWrtPressure( mesh, regionNames, dofManager, localMatrix, localRhs ); + assembleFluidMassResidualDerivativeWrtDisplacement( mesh, regionNames, dofManager, localMatrix, localRhs ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +setUpDflux_dApertureMatrix( DomainPartition & domain, + DofManager const & GEOS_UNUSED_PARAM( dofManager ), + CRSMatrix< real64, globalIndex > & localMatrix ) +{ + integer const numComp = this->flowSolver()->numFluidComponents(); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & regionNames ) + { + std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = this->getRefDerivativeFluxResidual_dAperture(); + + { + // calculate number of fracture elements + localIndex numRows = 0; + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, FaceElementSubRegion const & subRegion ) + { + numRows += subRegion.size(); + } ); + // number of columns (derivatives) = number of fracture elements + localIndex numCol = numRows; + // number of rows (equations) = number of fracture elements * number of components + numRows *= numComp; + + derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numCol ); + derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" ); + + derivativeFluxResidual_dAperture->reserveNonZeros( localMatrix.numNonZeros() ); + localIndex maxRowSize = -1; + for( localIndex row = 0; row < localMatrix.numRows(); ++row ) + { + localIndex const rowSize = localMatrix.numNonZeros( row ); + maxRowSize = maxRowSize > rowSize ? maxRowSize : rowSize; + } + // TODO This is way too much. The With the full system rowSize is not a good estimate for this. + for( localIndex row = 0; row < numRows * numComp; ++row ) + { + derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize ); + } + } + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) + { + for( integer ic = 0; ic < numComp; ic++ ) + { + derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0] * numComp + ic, sei[iconn][k1], 0.0 ); + } + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +addTransmissibilityCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, // meshBodyName, + MeshLevel const & mesh, + arrayView1d< string const > const & ) // regionNames + { + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const flowDofKey = dofManager.getKey( m_flowDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & stabilizationMethod = fvManager.getFluxApproximation( this->solidMechanicsSolver()->getStabilizationName() ); + + stabilizationMethod.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn=0; iconn( sesri[iconn][0] ); + + ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst(); + + arrayView1d< globalIndex const > const faceElementDofNumber = + elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + for( localIndex k0=0; k0= 0 && rowNumber < rowLengths.size() ) + { + for( localIndex k1=0; k1 +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +addTransmissibilityCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const +{ + GEOS_MARK_FUNCTION; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_flowDofKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + // Get the finite volume method used to compute the stabilization + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + GEOS_ERROR_IF( !fractureSubRegion.hasWrapper( flow::pressure::key() ), + this->getDataContext() << ": The fracture subregion must contain pressure field." ); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + + arrayView1d< globalIndex const > const & + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + + // A fracture connector has to be an edge shared by two faces + if( numFluxElems == 2 ) + { + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + // First index: face element. Second index: node + for( localIndex kf = 0; kf < 2; ++kf ) + { + // Set row DOF index + // Note that the 1-kf index is intentional, as this is coupling the pressure of one face cell + // to the nodes of the adjacent cell + localIndex const rowIndex = flowDofNumber[sei[iconn][1-kf]] - rankOffset; + + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + + // Get fracture, face and region/subregion/element indices (for elements on both sides) + localIndex const fractureIndex = sei[iconn][kf]; + + // Get the number of nodes + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] ); + + // Loop over the two sides of each fracture element + for( localIndex kf1 = 0; kf1 < 2; ++kf1 ) + { + localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1]; + + // Save the list of DOF associated with nodes + for( localIndex a=0; a( i ); + pattern.insertNonZero( rowIndex, colIndex ); + } + } + } + } + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +assembleForceResidualDerivativeWrtPressure( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + EdgeManager const & edgeManager = mesh.getEdgeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); + arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); + arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); + arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + arrayView1d< real64 const > faceAreas = faceManager.faceArea(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_flowDofKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + // Get the coordinates for all nodes + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + arrayView1d< globalIndex const > const & + flowDofNumber = subRegion.getReference< globalIndex_array >( flowDofKey ); + arrayView1d< real64 const > const & pressure = subRegion.getReference< array1d< real64 > >( flow::pressure::key() ); + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + + forAll< serialPolicy >( subRegion.size(), [=]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + + real64 Nbar[3]; + Nbar[ 0 ] = faceNormal[elemsToFaces[kfe][0]][0] - faceNormal[elemsToFaces[kfe][1]][0]; + Nbar[ 1 ] = faceNormal[elemsToFaces[kfe][0]][1] - faceNormal[elemsToFaces[kfe][1]][1]; + Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; + LvArray::tensorOps::normalize< 3 >( Nbar ); + globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types + real64 nodeRHS[3 * m_maxFaceNodes]; + stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); + globalIndex colDOF[1]; + colDOF[0] = flowDofNumber[kfe]; // pressure is always first + + for( localIndex kf=0; kf<2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe][kf]; + + // Compute local area contribution for each node + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + for( localIndex a=0; a( globalNodalForce, Nbar, nodalForceMag ); + + for( localIndex i=0; i<3; ++i ) + { + rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf ); + + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf ); + } + } + + for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRow< parallelHostAtomic >( localRow, + colDOF, + &dRdP[idof], + 1 ); + RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] ); + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) +{ + GEOS_MARK_FUNCTION; + + integer const numComp = this->flowSolver()->numFluidComponents(); + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + EdgeManager const & edgeManager = mesh.getEdgeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); + arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); + arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); + arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + arrayView1d< real64 const > faceAreas = faceManager.faceArea(); + + CRSMatrixView< real64 const, localIndex const > const & + dFluxResidual_dNormalJump = getDerivativeFluxResidual_dNormalJump().toViewConst(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_flowDofKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + // Get the coordinates for all nodes + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >(); + + arrayView1d< globalIndex const > const & flowDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + arrayView1d< real64 const > const & area = subRegion.getElementArea().toViewConst(); + + arrayView1d< integer const > const & fractureState = subRegion.getField< fields::contact::fractureState >(); + + forAll< serialPolicy >( subRegion.size(), [&]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + globalIndex nodeDOF[2*3*m_maxFaceNodes]; + + real64 Nbar[3]; + Nbar[ 0 ] = faceNormal[kf0][0] - faceNormal[kf1][0]; + Nbar[ 1 ] = faceNormal[kf0][1] - faceNormal[kf1][1]; + Nbar[ 2 ] = faceNormal[kf0][2] - faceNormal[kf1][2]; + LvArray::tensorOps::normalize< 3 >( Nbar ); + + stackArray2d< real64, 2*3*m_maxFaceNodes * MultiFluidBase::MAX_NUM_COMPONENTS > dRdU( MultiFluidBase::MAX_NUM_COMPONENTS, 2*3*m_maxFaceNodes ); + + bool const isFractureOpen = ( fractureState[kfe] == fields::contact::FractureState::Open ); + + // Accumulation derivative + if( isFractureOpen ) + { + for( localIndex kf=0; kf<2; ++kf ) + { + // Compute local area contribution for each node + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + + // TODO: move to something like this plus a static method. + // localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[kfe][kf] ); + // stackArray1d nodalArea( numNodesPerFace ); + + for( localIndex a=0; a( i ); + real64 const dAper_dU = -pow( -1, kf ) * Nbar[i]; + for( integer ic = 0; ic < numComp; ic++ ) + { + dRdU[ ic ][ kf*3*numNodesPerFace + 3*a + i ] = dVolume_dAperture * dAper_dU * compDens[kfe][ic]; // assuming poro=1 + } + } + } + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( flowDofNumber[kfe] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + integer const numRows = numComp; + for( integer i = 0; i < numRows; ++i ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow+i, + nodeDOF, + dRdU[i], + 2 * 3 * numNodesPerFace ); + } + } + } + + // flux derivative + bool skipAssembly = true; + for( integer ic = 0; ic < numComp; ic++ ) + { + localIndex const numColumns = dFluxResidual_dNormalJump.numNonZeros( kfe * numComp + ic ); + arraySlice1d< localIndex const > const & columns = dFluxResidual_dNormalJump.getColumns( kfe * numComp + ic ); + // TODO uncomment and fix build + arraySlice1d< real64 const > const & values = dFluxResidual_dNormalJump.getEntries( kfe * numComp + ic ); + + skipAssembly &= !isFractureOpen; + + for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) + { + real64 const dR_dAper = values[kfe1]; + localIndex const kfe2 = columns[kfe1]; + + bool const isOpen = ( fractureState[kfe2] == fields::contact::FractureState::Open ); + skipAssembly &= !isOpen; + + for( localIndex kf = 0; kf < 2; ++kf ) + { + //TODO: We should avoid allocating LvArrays inside kernel + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe2][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + + for( localIndex a=0; a( i ); + real64 const dAper_dU = -pow( -1, kf ) * Nbar[i] * ( nodalArea[a] / area[kfe2] ); + dRdU[ ic ][ kf*3*numNodesPerFace + 3*a + i ] = dR_dAper * dAper_dU; + } + } + } + } + + if( !skipAssembly ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( flowDofNumber[kfe] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + integer const numRows = numComp; + for( integer i = 0; i < numRows; ++i ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow+i, + nodeDOF, + dRdU[i], + 2 * 3 * numNodesPerFace ); + } + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::updateState( DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + + // call base poromechanics update + Base::updateState( domain ); + // need to call solid mechanics update separately to compute face displacement jump + this->solidMechanicsSolver()->updateState( domain ); + + // remove the contribution of the hydraulic aperture from the stencil weights + this->flowSolver()->prepareStencilWeights( domain ); + + updateHydraulicApertureAndFracturePermeability( domain ); + + // update the stencil weights using the updated hydraulic aperture + this->flowSolver()->updateStencilWeights( domain ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ) +{ + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion & subRegion ) + { + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + arrayView1d< real64 const > const area = subRegion.getElementArea(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView2d< real64 const > const fractureTraction = subRegion.getField< fields::contact::traction >(); + arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< fields::flow::aperture0 >(); + + arrayView1d< real64 > const aperture = subRegion.getElementAperture(); + arrayView1d< real64 > const hydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); + arrayView1d< real64 > const deltaVolume = subRegion.getField< flow::deltaVolume >(); + + string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( porousSolidName ); + + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); + + constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicAperture ) + { + using HydraulicApertureType = TYPEOFREF( castedHydraulicAperture ); + typename HydraulicApertureType::KernelWrapper hydraulicApertureWrapper = castedHydraulicAperture.createKernelWrapper(); + + constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + { + typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); + + poromechanicsFracturesKernels::StateUpdateKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + porousMaterialWrapper, + hydraulicApertureWrapper, + dispJump, + pressure, + area, + volume, + deltaVolume, + aperture, + oldHydraulicAperture, + hydraulicAperture, + fractureTraction ); + + } ); + } ); + } ); + } ); +} + +template class MultiphasePoromechanicsConformingFractures<>; +//template class MultiphasePoromechanicsConformingFractures< MultiphaseReservoirAndWells<> >; + +namespace +{ +//typedef MultiphasePoromechanicsConformingFractures< MultiphaseReservoirAndWells<> > +// MultiphasePoromechanicsConformingFractures; +//REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphasePoromechanicsConformingFractures, string const &, Group * const ) +typedef MultiphasePoromechanicsConformingFractures<> MultiphasePoromechanicsConformingFractures; +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphasePoromechanicsConformingFractures, string const &, Group * const ) +} + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp new file mode 100644 index 00000000000..cd8d7c0073b --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp @@ -0,0 +1,203 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultiphasePoromechanicsConformingFractures.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ + +#include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp" +#include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" + +namespace geos +{ + +template< typename FLOW_SOLVER = CompositionalMultiphaseBase > +class MultiphasePoromechanicsConformingFractures : public MultiphasePoromechanics< FLOW_SOLVER, SolidMechanicsLagrangeContact > +{ +public: + + using Base = MultiphasePoromechanics< FLOW_SOLVER, SolidMechanicsLagrangeContact >; + using Base::m_solvers; + using Base::m_dofManager; + using Base::m_localMatrix; + using Base::m_rhs; + using Base::m_solution; + + /// String used to form the solverName used to register solvers in CoupledSolver + static string coupledSolverAttributePrefix() { return "poromechanicsConformingFractures"; } + + /** + * @brief main constructor for MultiphasePoromechanicsConformingFractures objects + * @param name the name of this instantiation of MultiphasePoromechanicsConformingFractures in the repository + * @param parent the parent group of this instantiation of MultiphasePoromechanicsConformingFractures + */ + MultiphasePoromechanicsConformingFractures( const string & name, + dataRepository::Group * const parent ); + + /// Destructor for the class + ~MultiphasePoromechanicsConformingFractures() override {} + + /** + * @brief name of the node manager in the object catalog + * @return string that contains the catalog name to generate a new SinglePhasePoromechanicsConformingFractures object through the object + * catalog. + */ + static string catalogName() + { + if constexpr ( std::is_same_v< FLOW_SOLVER, CompositionalMultiphaseBase > ) + { + return "MultiphasePoromechanicsConformingFractures"; + } + else + { + return FLOW_SOLVER::catalogName() + "PoromechanicsConformingFractures"; + } + } + + /** + * @copydoc PhysicsSolverBase::getCatalogName() + */ + string getCatalogName() const override { return catalogName(); } + + /** + * @defgroup Solver Interface Functions + * + * These functions provide the primary interface that is required for derived classes + */ + /**@{*/ + + virtual void setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const override; + + virtual void setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity = true ) override; + + virtual void assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; + + virtual void updateState( DomainPartition & domain ) override final; + + /**@}*/ + +protected: + + virtual void postInputInitialization() override; + +private: + + struct viewKeyStruct : public Base::viewKeyStruct + {}; + + static const localIndex m_maxFaceNodes=11; // Maximum number of nodes on a contact face + + void assembleElementBasedContributions( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + virtual void assembleCouplingTerms( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; + + void assembleForceResidualDerivativeWrtPressure( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + /** + * @Brief add the nnz induced by the flux-aperture coupling + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param rowLenghts the nnz in each row + */ + void addTransmissibilityCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const; + + /** + * @Brief add the sparsity pattern induced by the flux-aperture coupling + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param pattern the sparsity pattern + */ + void addTransmissibilityCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const; + + /** + * @brief Set up the Dflux_dApertureMatrix object + * + * @param domain + * @param dofManager + * @param localMatrix + */ + void setUpDflux_dApertureMatrix( DomainPartition & domain, + DofManager const & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix ); + + /** + * @brief + * + * @param domain + */ + void updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ); + + + std::unique_ptr< CRSMatrix< real64, localIndex > > & getRefDerivativeFluxResidual_dAperture() + { + return m_derivativeFluxResidual_dAperture; + } + + CRSMatrixView< real64, localIndex const > getDerivativeFluxResidual_dNormalJump() + { + return m_derivativeFluxResidual_dAperture->toViewConstSizes(); + } + + CRSMatrixView< real64 const, localIndex const > getDerivativeFluxResidual_dNormalJump() const + { + return m_derivativeFluxResidual_dAperture->toViewConst(); + } + + std::unique_ptr< CRSMatrix< real64, localIndex > > m_derivativeFluxResidual_dAperture; + + string const m_flowDofKey = CompositionalMultiphaseBase::viewKeyStruct::elemDofFieldString(); + +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index 02fa02be59b..11ec7caa5d4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -24,6 +24,7 @@ #include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsStatistics.hpp" #include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp" @@ -136,6 +137,7 @@ execute( real64 const time_n, namespace { typedef PoromechanicsInitialization< MultiphasePoromechanics<> > MultiphasePoromechanicsInitialization; +typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFractures<> > MultiphasePoromechanicsConformingFracturesInitialization; typedef PoromechanicsInitialization< MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > > MultiphaseReservoirPoromechanicsInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanics<> > SinglePhasePoromechanicsInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures<> > SinglePhasePoromechanicsConformingFracturesInitialization; @@ -144,6 +146,7 @@ typedef PoromechanicsInitialization< SinglePhasePoromechanicsEmbeddedFractures > typedef PoromechanicsInitialization< SinglePhasePoromechanics< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsInitialization; typedef PoromechanicsInitialization< HydrofractureSolver< SinglePhasePoromechanics<> > > HydrofractureInitialization; REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, MultiphaseReservoirPoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 62a772e0278..9db368c00b7 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -21,10 +21,10 @@ #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_ -#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/multiphysics/CoupledSolver.hpp" -#include "physicsSolvers/multiphysics/PoromechanicsFields.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/multiphysics/PoromechanicsFields.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "constitutive/solid/PorousSolid.hpp" #include "constitutive/contact/HydraulicApertureBase.hpp" diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp index 0e657441f11..1965dc5a457 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp @@ -21,10 +21,7 @@ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ #include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" -#include "physicsSolvers/multiphysics/CoupledSolver.hpp" #include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" -#include "dataRepository/Group.hpp" namespace geos { diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index 6c6e2be5bc4..f25329cfe14 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -95,17 +95,10 @@ void SinglePhasePoromechanicsEmbeddedFractures::initializePostInitialConditionsP updateState( this->getGroupByPath< DomainPartition >( "/Problem/domain" ) ); } -void SinglePhasePoromechanicsEmbeddedFractures::setupDofs( DomainPartition const & domain, - DofManager & dofManager ) const +void SinglePhasePoromechanicsEmbeddedFractures::setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const { - GEOS_MARK_FUNCTION; - solidMechanicsSolver()->setupDofs( domain, dofManager ); - flowSolver()->setupDofs( domain, dofManager ); - - // Add coupling between displacement and cell pressures - dofManager.addCoupling( fields::solidMechanics::totalDisplacement::key(), - SinglePhaseBase::viewKeyStruct::elemDofFieldString(), - DofManager::Connector::Elem ); + Base::setupCoupling( domain, dofManager ); map< std::pair< string, string >, array1d< string > > meshTargets; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp index 95c3e7b7f82..f90d4ffe0be 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -58,8 +58,8 @@ class SinglePhasePoromechanicsEmbeddedFractures : public SinglePhasePoromechanic bool const setSparsity = true ) override; virtual void - setupDofs( DomainPartition const & domain, - DofManager & dofManager ) const override; + setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const override; virtual void assembleSystem( real64 const time, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp new file mode 100644 index 00000000000..33b07aaaaa1 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp @@ -0,0 +1,410 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultiphasePoromechanicsConformingFractures.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP + +#include "physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp" +//#include "physicsSolvers/fluidFlow/FluxKernelsHelper.hpp" + +namespace geos +{ + +namespace multiphasePoromechanicsConformingFracturesKernels +{ + +template< integer NUM_EQN, integer NUM_DOF > +class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper > +{ +public: + + /** + * @brief The type for element-based data. Consists entirely of ArrayView's. + * + * Can be converted from ElementRegionManager::ElementViewConstAccessor + * by calling .toView() or .toViewConst() on an accessor instance + */ + template< typename VIEWTYPE > + using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + + using AbstractBase = isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernelBase; + using DofNumberAccessor = AbstractBase::DofNumberAccessor; + using CompFlowAccessors = AbstractBase::CompFlowAccessors; + using MultiFluidAccessors = AbstractBase::MultiFluidAccessors; + using CapPressureAccessors = AbstractBase::CapPressureAccessors; + using PermeabilityAccessors = AbstractBase::PermeabilityAccessors; + using FracturePermeabilityAccessors = StencilMaterialAccessors< constitutive::PermeabilityBase, + fields::permeability::dPerm_dDispJump >; + + using AbstractBase::m_dt; + using AbstractBase::m_dofNumber; + using AbstractBase::m_gravCoef; + using AbstractBase::m_pres; + + using Base = isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper >; + using Base::numDof; + using Base::numEqn; + using Base::maxNumElems; + using Base::maxNumConns; + using Base::maxStencilSize; + using Base::m_stencilWrapper; + using Base::m_seri; + using Base::m_sesri; + using Base::m_sei; + using Base::numFluxSupportPoints; + using Base::numComp; + using Base::m_numPhases; + using Base::m_kernelFlags; + + using Base::m_permeability; + using Base::m_dPerm_dPres; + using Base::m_phaseMob; + using Base::m_dPhaseMob; + using Base::m_phaseVolFrac; + using Base::m_dPhaseVolFrac; + using Base::m_phaseCompFrac; + using Base::m_dPhaseCompFrac; + using Base::m_dCompFrac_dCompDens; + using Base::m_phaseMassDens; + using Base::m_dPhaseMassDens; + using Base::m_phaseCapPressure; + using Base::m_dPhaseCapPressure_dPhaseVolFrac; + + FluxComputeKernel( integer const numPhases, + globalIndex const rankOffset, + SurfaceElementStencilWrapper const & stencilWrapper, + DofNumberAccessor const & dofNumberAccessor, + CompFlowAccessors const & compFlowAccessors, + MultiFluidAccessors const & multiFluidAccessors, + CapPressureAccessors const & capPressureAccessors, + PermeabilityAccessors const & permeabilityAccessors, + FracturePermeabilityAccessors const & fracturePermeabilityAccessors, + real64 const dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + : Base( numPhases, + rankOffset, + stencilWrapper, + dofNumberAccessor, + compFlowAccessors, + multiFluidAccessors, + capPressureAccessors, + permeabilityAccessors, + dt, + localMatrix, + localRhs, + kernelFlags ), + m_dR_dAper( dR_dAper ), + m_dPerm_dDispJump( fracturePermeabilityAccessors.get( fields::permeability::dPerm_dDispJump {} ) ) + {} + + + /** + * @struct StackVariables + * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack + */ + struct StackVariables : public Base::StackVariables + { +public: + + /** + * @brief Constructor for the stack variables + * @param[in] size size of the stencil for this connection + * @param[in] numElems number of elements for this connection + */ + GEOS_HOST_DEVICE + StackVariables( localIndex const size, localIndex numElems ) + : Base::StackVariables( size, numElems ), + localColIndices( numElems ), + // TODO + dFlux_dAperture( numComp, numElems, size ) + {} + + stackArray1d< localIndex, maxNumElems > localColIndices; + + // TODO + stackArray3d< real64, numEqn * maxNumElems * maxStencilSize > dFlux_dAperture; + + /// Derivatives of transmissibility with respect to the dispJump + real64 dTrans_dDispJump[maxNumConns][2][3]{}; + }; + + /** + * @brief Performs the setup phase for the kernel. + * @param[in] iconn the connection index + * @param[in] stack the stack variables + */ + GEOS_HOST_DEVICE + void setup( localIndex const iconn, + StackVariables & stack ) const + { + // set degrees of freedom indices for this face + for( integer i = 0; i < stack.stencilSize; ++i ) + { + globalIndex const offset = m_dofNumber[m_seri( iconn, i )][m_sesri( iconn, i )][m_sei( iconn, i )]; + for( integer jdof = 0; jdof < numDof; ++jdof ) + { + stack.dofColIndices[i * numDof + jdof] = offset + jdof; + } + stack.localColIndices[ i ] = m_sei( iconn, i ); + } + } + + /** + * @brief Compute the local flux contributions to the residual and Jacobian + * @tparam FUNC the type of the function that can be used to customize the computation of the flux + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + * @param[in] NoOpFunc the function used to customize the computation of the flux + */ + template< typename FUNC = NoOpFunc > + GEOS_HOST_DEVICE + void computeFlux( localIndex const iconn, + StackVariables & stack, + FUNC && compFluxKernelOp = NoOpFunc{} ) const + { + m_stencilWrapper.computeWeights( iconn, + m_permeability, + m_dPerm_dPres, + m_dPerm_dDispJump, + stack.transmissibility, + stack.dTrans_dPres, + stack.dTrans_dDispJump ); + + localIndex k[numFluxSupportPoints]; + localIndex connectionIndex = 0; + for( k[0] = 0; k[0] < stack.numConnectedElems; ++k[0] ) + { + for( k[1] = k[0] + 1; k[1] < stack.numConnectedElems; ++k[1] ) + { + /// cell indices + localIndex const seri[numFluxSupportPoints] = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )}; + localIndex const sesri[numFluxSupportPoints] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )}; + localIndex const sei[numFluxSupportPoints] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )}; + + // clear working arrays + real64 compFlux[numComp]{}; + real64 dCompFlux_dP[numFluxSupportPoints][numComp]{}; + real64 dCompFlux_dC[numFluxSupportPoints][numComp][numComp]{}; + real64 dCompFlux_dTrans[numComp]{}; + + real64 const trans[numFluxSupportPoints] = { stack.transmissibility[connectionIndex][0], + stack.transmissibility[connectionIndex][1] }; + + real64 const dTrans_dPres[numFluxSupportPoints] = { stack.dTrans_dPres[connectionIndex][0], + stack.dTrans_dPres[connectionIndex][1] }; + + //***** calculation of flux ***** + // loop over phases, compute and upwind phase flux and sum contributions to each component's flux + for( integer ip = 0; ip < m_numPhases; ++ip ) + { + // create local work arrays + real64 potGrad = 0.0; + real64 phaseFlux = 0.0; + real64 dPhaseFlux_dP[numFluxSupportPoints]{}; + real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{}; + + localIndex k_up = -1; + + + isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute< numComp, numFluxSupportPoints > + ( m_numPhases, + ip, + m_kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure ), + m_kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CheckPhasePresenceInGravity ), + seri, sesri, sei, + trans, + dTrans_dPres, + m_pres, + m_gravCoef, + m_phaseMob, m_dPhaseMob, + m_phaseVolFrac, + m_dPhaseVolFrac, + m_phaseCompFrac, m_dPhaseCompFrac, + m_dCompFrac_dCompDens, + m_phaseMassDens, m_dPhaseMassDens, + m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, + k_up, + potGrad, + phaseFlux, + dPhaseFlux_dP, + dPhaseFlux_dC, + compFlux, + dCompFlux_dP, + dCompFlux_dC, + dCompFlux_dTrans ); + + // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives + // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase + compFluxKernelOp( ip, k, seri, sesri, sei, connectionIndex, + k_up, seri[k_up], sesri[k_up], sei[k_up], potGrad, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + + } // loop over phases + + /// populate local flux vector and derivatives + for( integer ic = 0; ic < numComp; ++ic ) + { + integer const eqIndex0 = k[0] * numEqn + ic; + integer const eqIndex1 = k[1] * numEqn + ic; + + stack.localFlux[eqIndex0] += m_dt * compFlux[ic]; + stack.localFlux[eqIndex1] -= m_dt * compFlux[ic]; + + for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) + { + localIndex const localDofIndexPres = k[ke] * numDof; + stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dCompFlux_dP[ke][ic]; + stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dCompFlux_dP[ke][ic]; + + for( integer jc = 0; jc < numComp; ++jc ) + { + localIndex const localDofIndexComp = localDofIndexPres + jc + 1; + stack.localFluxJacobian[eqIndex0][localDofIndexComp] += m_dt * dCompFlux_dC[ke][ic][jc]; + stack.localFluxJacobian[eqIndex1][localDofIndexComp] -= m_dt * dCompFlux_dC[ke][ic][jc]; + } + } + } + + for( integer ic = 0; ic < numComp; ++ic ) + { + real64 dFlux_dAper[2] = {0.0, 0.0}; + dFlux_dAper[0] = m_dt * dCompFlux_dTrans[ic] * stack.dTrans_dDispJump[connectionIndex][0][0]; + dFlux_dAper[1] = -m_dt * dCompFlux_dTrans[ic] * stack.dTrans_dDispJump[connectionIndex][1][0]; + + stack.dFlux_dAperture[ic][k[0]][k[0]] += dFlux_dAper[0]; + stack.dFlux_dAperture[ic][k[0]][k[1]] += dFlux_dAper[1]; + stack.dFlux_dAperture[ic][k[1]][k[0]] -= dFlux_dAper[0]; + stack.dFlux_dAperture[ic][k[1]][k[1]] -= dFlux_dAper[1]; + } +// + connectionIndex++; + } + } + } + + /** + * @brief Performs the complete phase for the kernel. + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + */ + template< typename FUNC = NoOpFunc > + GEOS_HOST_DEVICE + void complete( localIndex const iconn, + StackVariables & stack, + FUNC && kernelOp = NoOpFunc{} ) const + { + // Call Base::complete to assemble the mass balance equations + // In the lambda, fill the dR_dAper matrix + Base::complete( iconn, stack, [&] ( integer const i, + localIndex const localRow ) + { + // TODO + localIndex const ei = LvArray::integerConversion< localIndex >( m_sei( iconn, i ) ); + for( integer ic = 0; ic < numComp; ++ic ) + { + localIndex const row = ei * numComp + ic; + + m_dR_dAper.addToRowBinarySearch< parallelDeviceAtomic >( row, + stack.localColIndices.data(), + stack.dFlux_dAperture[ic][i].dataIfContiguous(), + stack.stencilSize ); + } + + // call the lambda to assemble additional terms, such as thermal terms + kernelOp( i, localRow ); + } ); + } + +private: + + CRSMatrixView< real64, localIndex const > m_dR_dAper; + + ElementViewConst< arrayView4d< real64 const > > const m_dPerm_dDispJump; +}; + + + +/** + * @class FluxComputeKernelFactory + */ +class FluxComputeKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @tparam STENCILWRAPPER the type of the stencil wrapper + * @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) + * @param[in] elemManager reference to the element region manager + * @param[in] stencilWrapper reference to the stencil wrapper + * @param[in] dt time step size + * @param[inout] localMatrix the local CRS matrix + * @param[inout] localRhs the local right-hand side vector + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComps, + integer const numPhases, + globalIndex const rankOffset, + string const & dofKey, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, + string const & solverName, + ElementRegionManager const & elemManager, + SurfaceElementStencilWrapper const & stencilWrapper, + real64 const dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + { + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC ) + { + integer constexpr NUM_COMP = NC(); + integer constexpr NUM_DOF = NC() + 1; + + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); + dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + + using kernelType = FluxComputeKernel< NUM_COMP, NUM_DOF >; + typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); + typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); + typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName ); + typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); + typename kernelType::FracturePermeabilityAccessors fracPermAccessors( elemManager, solverName ); + + kernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor, + compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors, fracPermAccessors, + dt, localMatrix, localRhs, kernelFlags, dR_dAper ); + + kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); + } ); + } +}; + +} // namespace multiphasePoromechanicsConformingFracturesKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp index 302c8113b3d..00350d1f905 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp @@ -22,6 +22,7 @@ #include "physicsSolvers/multiphysics/PoromechanicsFields.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsBase.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" namespace geos { diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index e130b1e56d4..27cb9f836a9 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -3245,7 +3245,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m - + @@ -3268,7 +3268,7 @@ Local- Add jump stabilization on interior of macro elements--> - + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 3abf5b863f0..64384a62a01 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -486,7 +486,7 @@ - + @@ -572,7 +572,7 @@ - + @@ -609,7 +609,7 @@ - + @@ -660,7 +660,7 @@ - + @@ -701,7 +701,7 @@ - + @@ -734,7 +734,7 @@ - + @@ -745,7 +745,7 @@ - + @@ -758,7 +758,7 @@ - + @@ -771,7 +771,7 @@ - + @@ -787,7 +787,7 @@ - + @@ -821,7 +821,7 @@ - + @@ -884,7 +884,7 @@ - + @@ -915,7 +915,7 @@ - + @@ -928,7 +928,7 @@ - + @@ -941,7 +941,7 @@ - + @@ -954,7 +954,7 @@ - + @@ -967,7 +967,7 @@ - + @@ -982,7 +982,7 @@ - + @@ -993,7 +993,7 @@ - + @@ -1006,7 +1006,7 @@ - + @@ -1017,7 +1017,7 @@ - + @@ -1028,7 +1028,7 @@ - + @@ -1039,7 +1039,7 @@ - + @@ -1050,7 +1050,7 @@ - + @@ -1063,7 +1063,7 @@ - + @@ -1074,7 +1074,7 @@ - + @@ -1085,7 +1085,7 @@ - + @@ -1098,7 +1098,7 @@ - + @@ -1113,7 +1113,7 @@ - + @@ -1128,7 +1128,7 @@ - + @@ -1141,7 +1141,7 @@ - + @@ -1156,7 +1156,7 @@ - + @@ -1167,7 +1167,7 @@ - + @@ -1180,7 +1180,7 @@ - + @@ -1193,7 +1193,7 @@ - + @@ -1208,7 +1208,7 @@ - + @@ -1224,7 +1224,7 @@ - + @@ -1239,7 +1239,7 @@ - + @@ -1256,7 +1256,7 @@ - + @@ -1273,7 +1273,7 @@ - + @@ -1290,7 +1290,7 @@ - + @@ -1305,7 +1305,7 @@ - + @@ -1318,7 +1318,7 @@ - + @@ -1357,7 +1357,7 @@ - + @@ -1386,7 +1386,7 @@ - + @@ -1479,7 +1479,7 @@ - + @@ -3103,7 +3103,7 @@ - + @@ -3131,7 +3131,7 @@ - + @@ -3150,11 +3150,11 @@ - + - + @@ -3164,7 +3164,7 @@ - + @@ -3174,11 +3174,11 @@ - + - + @@ -3188,7 +3188,7 @@ - + @@ -3198,7 +3198,7 @@ - + @@ -3208,7 +3208,7 @@ - + @@ -3232,7 +3232,7 @@ - + @@ -3250,7 +3250,7 @@ - + @@ -3262,7 +3262,7 @@ - + @@ -3274,7 +3274,7 @@ - + @@ -3282,11 +3282,11 @@ - + - + @@ -3309,7 +3309,7 @@ - + @@ -3335,7 +3335,7 @@ - + @@ -3356,7 +3356,7 @@ - + @@ -3386,7 +3386,7 @@ - + @@ -3400,7 +3400,7 @@ - + @@ -3427,7 +3427,7 @@ - + @@ -3466,7 +3466,7 @@ - + diff --git a/src/coreComponents/unitTests/dataRepositoryTests/testObjectCatalog.cpp b/src/coreComponents/unitTests/dataRepositoryTests/testObjectCatalog.cpp index 586e4fb3bc3..737981ea002 100644 --- a/src/coreComponents/unitTests/dataRepositoryTests/testObjectCatalog.cpp +++ b/src/coreComponents/unitTests/dataRepositoryTests/testObjectCatalog.cpp @@ -17,6 +17,7 @@ // Source includes #include "dataRepository/ObjectCatalog.hpp" +#include "dataRepository/DataContext.hpp" #include "common/logger/Logger.hpp" #include "mainInterface/initialization.hpp" @@ -67,7 +68,6 @@ class Derived1 : public Base } static string catalogName() { return "derived1"; } string getCatalogName() { return catalogName(); } - }; REGISTER_CATALOG_ENTRY( Base, Derived1, int &, double const & ) //STOP_SPHINX @@ -99,14 +99,17 @@ TEST( testObjectCatalog, testRegistration ) GEOS_LOG( "EXECUTING MAIN" ); int junk = 1; double junk2 = 3.14; + dataRepository::DataFileContext const context = dataRepository::DataFileContext( "Base Test Class", __FILE__, __LINE__ ); // allocate a new Derived1 object std::unique_ptr< Base > - derived1 = Base::CatalogInterface::factory( "derived1", junk, junk2 ); + derived1 = Base::CatalogInterface::factory( "derived1", context, + junk, junk2 ); // allocate a new Derived2 object std::unique_ptr< Base > - derived2 = Base::CatalogInterface::factory( "derived2", junk, junk2 ); + derived2 = Base::CatalogInterface::factory( "derived2", context, + junk, junk2 ); EXPECT_STREQ( derived1->getCatalogName().c_str(), Derived1::catalogName().c_str() );