diff --git a/inputFiles/constitutiveDriver/relperm/relperm.ats b/inputFiles/constitutiveDriver/relperm/relperm.ats
new file mode 100644
index 00000000000..30c06d25a23
--- /dev/null
+++ b/inputFiles/constitutiveDriver/relperm/relperm.ats
@@ -0,0 +1,20 @@
+from geos.ats.test_builder import TestDeck, RestartcheckParameters, generate_geos_tests
+
+restartcheck_params = {'atol': 1e-08, 'rtol': 4e-07}
+
+def create_relperm_test(name, description=None):
+ if description is None: description = name
+ return TestDeck(
+ name=name,
+ description=description,
+ partitions=[(1, 1, 1)],
+ restart_step=0,
+ check_step=0,
+ restartcheck_params=RestartcheckParameters(**restartcheck_params))
+
+decks = [
+ create_relperm_test("testRelperm_Table_2_phase", "Relperm driver for 2 phase tables"),
+ create_relperm_test("testRelperm_Table_3_phase", "Relperm driver for 3 phase tables")
+]
+
+generate_geos_tests(decks)
diff --git a/inputFiles/constitutiveDriver/relperm/testRelperm_Table_2_phase.xml b/inputFiles/constitutiveDriver/relperm/testRelperm_Table_2_phase.xml
new file mode 100644
index 00000000000..94574aafe28
--- /dev/null
+++ b/inputFiles/constitutiveDriver/relperm/testRelperm_Table_2_phase.xml
@@ -0,0 +1,118 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/constitutiveDriver/relperm/testRelperm_Table_3_phase.xml b/inputFiles/constitutiveDriver/relperm/testRelperm_Table_3_phase.xml
new file mode 100644
index 00000000000..715e5e58421
--- /dev/null
+++ b/inputFiles/constitutiveDriver/relperm/testRelperm_Table_3_phase.xml
@@ -0,0 +1,41 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/constitutive/relativePermeability/RelativePermeabilitySelector.hpp b/src/coreComponents/constitutive/relativePermeability/RelativePermeabilitySelector.hpp
index ee20a4c4f6e..fbba10578b8 100644
--- a/src/coreComponents/constitutive/relativePermeability/RelativePermeabilitySelector.hpp
+++ b/src/coreComponents/constitutive/relativePermeability/RelativePermeabilitySelector.hpp
@@ -51,7 +51,7 @@ void constitutiveUpdatePassThru( RelativePermeabilityBase const & relPerm,
BrooksCoreyBakerRelativePermeability,
BrooksCoreyStone2RelativePermeability,
TableRelativePermeability,
- constitutive::TableRelativePermeabilityHysteresis,
+ TableRelativePermeabilityHysteresis,
VanGenuchtenBakerRelativePermeability,
VanGenuchtenStone2RelativePermeability >::execute( relPerm, std::forward< LAMBDA >( lambda ) );
}
@@ -64,7 +64,7 @@ void constitutiveUpdatePassThru( RelativePermeabilityBase & relPerm,
BrooksCoreyBakerRelativePermeability,
BrooksCoreyStone2RelativePermeability,
TableRelativePermeability,
- constitutive::TableRelativePermeabilityHysteresis,
+ TableRelativePermeabilityHysteresis,
VanGenuchtenBakerRelativePermeability,
VanGenuchtenStone2RelativePermeability >::execute( relPerm, std::forward< LAMBDA >( lambda ) );
}
diff --git a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.cpp b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.cpp
index ca17f97223f..4291eebf48e 100644
--- a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.cpp
+++ b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.cpp
@@ -13,204 +13,175 @@
* ------------------------------------------------------------------------------------------------------------
*/
-#include "common/MpiWrapper.hpp"
-#include "functions/FunctionManager.hpp"
-#include "functions/TableFunction.hpp"
+#include "RelpermDriver.hpp"
+
#include "constitutive/ConstitutiveManager.hpp"
#include "constitutiveDrivers/LogLevelsInfo.hpp"
#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp"
#include "constitutive/relativePermeability/RelativePermeabilitySelector.hpp"
-#include "RelpermDriver.hpp"
-
namespace geos
{
using namespace dataRepository;
using namespace constitutive;
-RelpermDriver::RelpermDriver( const geos::string & name,
- geos::dataRepository::Group * const parent )
- :
- TaskBase( name, parent )
+RelpermDriver::RelpermDriver( const string & name,
+ Group * const parent )
+ : ConstitutiveDriver( name, parent )
{
registerWrapper( viewKeyStruct::relpermNameString(), &m_relpermName ).
setRTTypeName( rtTypes::CustomTypes::groupNameRef ).
setInputFlag( InputFlags::REQUIRED ).
- setDescription( "Relperm model to test" );
+ setDescription( "Relative permeability model to test" );
- registerWrapper( viewKeyStruct::numStepsString(), &m_numSteps ).
- setInputFlag( InputFlags::REQUIRED ).
- setDescription( "Number of saturation steps to take" );
-
- registerWrapper( viewKeyStruct::outputString(), &m_outputFile ).
+ registerWrapper( viewKeyStruct::historicalSaturationsString(), &m_historicalSaturations ).
setInputFlag( InputFlags::OPTIONAL ).
- setApplyDefaultValue( "none" ).
- setDescription( "Output file" );
-
- registerWrapper( viewKeyStruct::baselineString(), &m_baselineFile ).
- setInputFlag( InputFlags::OPTIONAL ).
- setApplyDefaultValue( "none" ).
- setDescription( "Baseline file" );
-
- addLogLevel< logInfo::LogOutput >();
+ setDescription( "Historical saturations for each phase." );
}
-
-void RelpermDriver::outputResults()
+void RelpermDriver::postInputInitialization()
{
- // TODO: improve file path output to grab command line -o directory
- // for the moment, we just use the specified m_outputFile directly
+ ConstitutiveDriver::postInputInitialization();
- FILE * fp = fopen( m_outputFile.c_str(), "w" );
+ RelativePermeabilityBase const & baseRelperm = getRelperm();
- fprintf( fp, "# column 1 = time\n" );
- fprintf( fp, "# columns %d-%d = phase vol fractions\n", 2, 1 + m_numPhases );
- fprintf( fp, "# columns %d-%d = phase relperm\n", 2 + m_numPhases, 1 + 2 * m_numPhases );
+ integer const numPhases = baseRelperm.numFluidPhases();
- if( ( m_numPhases == 2 && m_table.size( 1 ) > 5 ) || m_table.size( 1 ) > 7 )
- {
- fprintf( fp, "# columns %d-%d = phase relperm (hyst)\n", 1 + 2 * m_numPhases, 1 + 3 * m_numPhases );
- }
+ // Must be 2-phase or 3-phase
+ GEOS_ERROR_IF( numPhases < 2 || 3 < numPhases,
+ "Number of phases for relative permeability model must be 2 or 3",
+ getWrapperDataContext( viewKeyStruct::relpermNameString() ) );
-
- for( integer n = 0; n < m_table.size( 0 ); ++n )
+ // Historical saturations must be the same number as the phases
+ if( !m_historicalSaturations.empty())
{
- for( integer col = 0; col < m_table.size( 1 ); ++col )
- {
- fprintf( fp, "%.4e ", m_table( n, col ) );
- }
- fprintf( fp, "\n" );
+ GEOS_ERROR_IF( m_historicalSaturations.size() != numPhases,
+ "Number of historical saturations must be the same as the number of phases",
+ getWrapperDataContext( viewKeyStruct::historicalSaturationsString() ) );
}
- fclose( fp );
-
-
-}
+ string_array columnNames;
+ getColumnNames( columnNames );
+ integer const numCols = static_cast< integer >(columnNames.size());
-void RelpermDriver::postInputInitialization()
-{
- ConstitutiveManager
- & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" );
- RelativePermeabilityBase
- & baseRelperm = constitutiveManager.getGroup< RelativePermeabilityBase >( m_relpermName );
-
- m_numPhases = baseRelperm.numFluidPhases();
-
+ allocateTable( numCols, numPhases );
}
-
-bool RelpermDriver::execute( const geos::real64 GEOS_UNUSED_PARAM( time_n ),
- const geos::real64 GEOS_UNUSED_PARAM( dt ),
- const geos::integer GEOS_UNUSED_PARAM( cycleNumber ),
- const geos::integer GEOS_UNUSED_PARAM( eventCounter ),
- const geos::real64 GEOS_UNUSED_PARAM( eventProgress ),
- geos::DomainPartition &
- GEOS_UNUSED_PARAM( domain ) )
+bool RelpermDriver::execute()
{
- // this code only makes sense in serial
-
- GEOS_THROW_IF( MpiWrapper::commRank() > 0, "RelpermDriver should only be run in serial", geos::RuntimeError );
+ RelativePermeabilityBase & baseRelperm = getRelperm();
-
- ConstitutiveManager
- & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" );
- RelativePermeabilityBase
- & baseRelperm = constitutiveManager.getGroup< RelativePermeabilityBase >( m_relpermName );
+ integer const numPhases = baseRelperm.numFluidPhases();
GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, "Launching Relperm Driver" );
- GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Relperm .................. " << m_relpermName );
+ GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Relperm ................ " << m_relpermName );
GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Type ................... " << baseRelperm.getCatalogName() );
- GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " No. of Phases .......... " << m_numPhases );
+ GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " No. of Phases .......... " << numPhases );
GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Steps .................. " << m_numSteps );
- GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Output ................. " << m_outputFile );
- GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Baseline ............... " << m_baselineFile );
+ if( !m_outputFile.empty())
+ {
+ GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Output ................. " << m_outputFile );
+ }
+
+ initializeTable();
// create a dummy discretization with one quadrature point for
// storing constitutive data
-
conduit::Node node;
dataRepository::Group rootGroup( "root", node );
dataRepository::Group discretization( "discretization", &rootGroup );
- discretization.resize( 1 ); // one element
+ // Allocate as many elements as the number of rows
+ integer const numRows = m_table.size( 0 );
+ discretization.resize( numRows ); // numRows elements
baseRelperm.allocateConstitutiveData( discretization, 1 ); // one quadrature point
constitutiveUpdatePassThru( baseRelperm, [&]( auto & selectedRelpermModel )
{
using RELPERM_TYPE = TYPEOFREF( selectedRelpermModel );
- resizeTables< RELPERM_TYPE >();
runTest< RELPERM_TYPE >( selectedRelpermModel, m_table );
} );
// move table back to host for output
m_table.move( LvArray::MemorySpace::host );
- if( m_outputFile != "none" )
+ return false;
+}
+
+void RelpermDriver::getColumnNames( string_array & columnNames ) const
+{
+ RelativePermeabilityBase const & baseRelperm = getRelperm();
+ bool const has_hysteresis = (dynamic_cast< constitutive::TableRelativePermeabilityHysteresis const * >(&baseRelperm) != nullptr);
+
+ integer const numPhases = baseRelperm.numFluidPhases();
+ string_array const & phaseNames = baseRelperm.phaseNames();
+
+ columnNames.emplace_back( "index" );
+ for( integer ip = 0; ip < numPhases; ip++ )
{
- outputResults();
+ columnNames.emplace_back( GEOS_FMT( "saturation,{}", phaseNames[ip] ));
}
-
- if( m_baselineFile != "none" )
+ if( has_hysteresis )
{
- compareWithBaseline();
+ for( integer ip = 0; ip < numPhases; ip++ )
+ {
+ columnNames.emplace_back( GEOS_FMT( "historical saturation,{}", phaseNames[ip] ));
+ }
+ }
+ for( integer ip = 0; ip < numPhases; ip++ )
+ {
+ columnNames.emplace_back( GEOS_FMT( "relperm,{}", phaseNames[ip] ));
}
-
- return false;
}
+void RelpermDriver::allocateTable( integer numColumns, integer numPhases )
+{
+ // For 3-phase we have m_numSteps+1 points for each of the other two phases
+ integer const numRows = (numPhases == 3) ? (m_numSteps+1)*(m_numSteps+1) : (m_numSteps+1);
+ m_table.resize( numRows, numColumns );
+ for( integer index = 0; index < numRows; ++index )
+ {
+ m_table( index, TIME ) = index;
+ }
+}
-template< typename RELPERM_TYPE >
-void RelpermDriver::resizeTables()
+void RelpermDriver::initializeTable()
{
- ConstitutiveManager
- & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" );
- RelativePermeabilityBase
- & baseRelperm = constitutiveManager.getGroup< RelativePermeabilityBase >( m_relpermName );
+ RelativePermeabilityBase & baseRelperm = getRelperm();
using PT = RelativePermeabilityBase::PhaseType;
integer const ipWater = baseRelperm.getPhaseOrder()[PT::WATER];
integer const ipOil = baseRelperm.getPhaseOrder()[PT::OIL];
integer const ipGas = baseRelperm.getPhaseOrder()[PT::GAS];
- real64 minSw = 0., minSnw = 0.;
- if( baseRelperm.numFluidPhases() > 2 )
- {
- minSw = baseRelperm.getWettingPhaseMinVolumeFraction();
- minSnw = baseRelperm.getNonWettingMinVolumeFraction();
- }
- else
- {
- if( ipWater < 0 )// a.k.a o/g
- {
- minSw = 0;
- minSnw = baseRelperm.getNonWettingMinVolumeFraction();
- }
- else if( ipGas < 0 || ipOil < 0 )// a.k.a w/o or w/g
- {
- minSnw = 0;
- minSw = baseRelperm.getWettingPhaseMinVolumeFraction();
- }
- }
+ integer const numPhases = baseRelperm.numFluidPhases();
+
+ auto const [ipWetting, ipNonWetting] = baseRelperm.wettingAndNonWettingPhaseIndices();
+ real64 const min_wetting_saturation = baseRelperm.getPhaseMinVolumeFraction()[ipWetting];
+ real64 const min_non_wetting_saturation = baseRelperm.getPhaseMinVolumeFraction()[ipNonWetting];
- real64 const dSw = ( 1 - minSw - minSnw ) / m_numSteps;
- // set input columns
+ real64 const dSw = ( 1.0 - min_wetting_saturation - min_non_wetting_saturation ) / m_numSteps;
+
+ // Offset for saturations in table
+ constexpr integer SATURATION = 1;
- resizeTable< RELPERM_TYPE >();
// 3-phase branch
- if( m_numPhases > 2 )
+ if( numPhases == 3 )
{
+ real64 swat = 0.0;
+ real64 sgas = 0.0;
for( integer ni = 0; ni < m_numSteps + 1; ++ni )
{
+ swat = min_wetting_saturation + ni*dSw;
for( integer nj = 0; nj < m_numSteps + 1; ++nj )
{
+ sgas = min_non_wetting_saturation + nj*dSw;
integer index = ni * ( m_numSteps + 1 ) + nj;
- m_table( index, TIME ) = minSw + index * dSw;
- m_table( index, ipWater + 1 ) = minSw + nj * dSw;
- m_table( index, ipGas + 1 ) = minSnw + ni * dSw;
- m_table( index, ipOil + 1 ) =
- 1. - m_table( index, ipWater + 1 ) - m_table( index, ipOil + 1 );
+ m_table( index, ipWater + SATURATION ) = swat;
+ m_table( index, ipGas + SATURATION ) = sgas;
+ m_table( index, ipOil + SATURATION ) = 1.0 - swat - sgas;
}
}
}
@@ -218,118 +189,23 @@ void RelpermDriver::resizeTables()
{
for( integer ni = 0; ni < m_numSteps + 1; ++ni )
{
- integer index = ni;
- m_table( index, TIME ) = minSw + index * dSw;
- if( ipWater < 0 )
- {
- m_table( index, ipGas + 1 ) = minSnw + ni * dSw;
- m_table( index, ipOil + 1 ) = 1. - m_table( index, ipGas + 1 );
- }
- else if( ipGas < 0 )
- {
- m_table( index, ipWater + 1 ) = minSw + ni * dSw;
- m_table( index, ipOil + 1 ) = 1. - m_table( index, ipWater + 1 );
- }
- else if( ipOil < 0 )
- {
- m_table( index, ipWater + 1 ) = minSw + ni * dSw;
- m_table( index, ipGas + 1 ) = 1. - m_table( index, ipWater + 1 );
- }
+ real64 const s_nw = min_non_wetting_saturation + ni * dSw;
+ m_table( ni, ipNonWetting + SATURATION ) = s_nw;
+ m_table( ni, ipWetting + SATURATION ) = 1.0 - s_nw;
}
-
}
-
-
}
-
-template< typename RELPERM_TYPE >
-std::enable_if_t< std::is_same< TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void >
-RelpermDriver::resizeTable()
+RelativePermeabilityBase & RelpermDriver::getRelperm()
{
- if( m_numPhases > 2 )
- {
- m_table.resize( ( m_numSteps + 1 ) * ( m_numSteps + 1 ), 1 + 3 * m_numPhases );
- }
- else
- {
- m_table.resize( m_numSteps + 1, 1 + 3 * m_numPhases );
- }
-
+ return getConstitutiveManager().getGroup< RelativePermeabilityBase >( m_relpermName );
}
-template< typename RELPERM_TYPE >
-std::enable_if_t< !std::is_same< TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void >
-RelpermDriver::resizeTable()
+RelativePermeabilityBase const & RelpermDriver::getRelperm() const
{
- if( m_numPhases > 2 )
- {
- m_table.resize( ( m_numSteps + 1 ) * ( m_numSteps + 1 ), 1 + 2 * m_numPhases );
- }
- else
- {
- m_table.resize( m_numSteps + 1, 1 + 2 * m_numPhases );
- }
-}
-
-
-//TODO refactor - duplication
-void RelpermDriver::compareWithBaseline()
-{
- // open baseline file
-
- std::ifstream file( m_baselineFile.c_str() );
- GEOS_THROW_IF( !file.is_open(),
- GEOS_FMT( "Can't seem to open the baseline file {}", m_baselineFile ),
- InputError );
-
- // discard file header
-
- string line;
- for( integer row = 0; row < 7; ++row )
- {
- getline( file, line );
- }
-
- // read data block. we assume the file size is consistent with m_table,
- // but check for a premature end-of-file. we then compare results value by value.
- // we ignore the newton iteration and residual columns, as those may be platform
- // specific.
-
- real64 value;
- //table is redim to fit the layout of relperm so the second dimension is numGaussPt
- // and always of size 1
- for( integer row = 0; row < m_table.size( 0 ); ++row )
- {
- for( integer col = 0; col < m_table.size( 1 ); ++col )
- {
- GEOS_THROW_IF( file.eof(), "Baseline file appears shorter than internal results", geos::RuntimeError );
- file >> value;
-
- real64 const error = fabs( m_table[row][col] - value ) / ( fabs( value ) + 1 );
- GEOS_THROW_IF( error > m_baselineTol,
- GEOS_FMT( "Results do not match baseline at data row {} (row {} with header) and column {}",
- row + 1,
- row + m_numColumns,
- col + 1 ),
- geos::RuntimeError );
- }
- }
-
- // check we actually reached the end of the baseline file
-
- file >> value;
- GEOS_THROW_IF( !file.eof(), "Baseline file appears longer than internal results", geos::RuntimeError );
-
- // success
-
- GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Comparison ............. Internal results consistent with baseline." );
-
- file.close();
+ return getConstitutiveManager().getGroup< RelativePermeabilityBase >( m_relpermName );
}
-
-
REGISTER_CATALOG_ENTRY( TaskBase,
RelpermDriver,
string const &, dataRepository::Group * const )
diff --git a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.hpp b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.hpp
index 426c66bf130..89e4f8b5cce 100644
--- a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.hpp
+++ b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriver.hpp
@@ -13,111 +13,73 @@
* ------------------------------------------------------------------------------------------------------------
*/
-#ifndef GEOS_RELPERMDRIVER_HPP_
-#define GEOS_RELPERMDRIVER_HPP_
+#ifndef GEOS_CONSTITUTIVEDRIVERS_RELATIVEPERMEABILITY_RELPERMDRIVER_HPP
+#define GEOS_CONSTITUTIVEDRIVERS_RELATIVEPERMEABILITY_RELPERMDRIVER_HPP
-#include "events/tasks/TaskBase.hpp"
-#include "constitutive/relativePermeability/TableRelativePermeabilityHysteresis.hpp"
+#include "constitutiveDrivers/ConstitutiveDriver.hpp"
namespace geos
{
-class RelpermDriver : public TaskBase
+namespace constitutive
{
+class RelativePermeabilityBase;
+}
+/**
+ * @class RelpermDriver
+ *
+ * Class to allow for testing Relative permeability models without the
+ * complexity of setting up a full simulation.
+ */
+class RelpermDriver : public ConstitutiveDriver
+{
public:
- RelpermDriver( const string & name,
- Group * const parent );
+ RelpermDriver( const string & name, Group * const parent );
- static string catalogName()
- { return "RelpermDriver"; }
+ static string catalogName() { return "RelpermDriver"; }
void postInputInitialization() override;
- virtual bool execute( real64 const GEOS_UNUSED_PARAM( time_n ),
- real64 const GEOS_UNUSED_PARAM( dt ),
- integer const GEOS_UNUSED_PARAM( cycleNumber ),
- integer const GEOS_UNUSED_PARAM( eventCounter ),
- real64 const GEOS_UNUSED_PARAM( eventProgress ),
- DomainPartition &
- GEOS_UNUSED_PARAM( domain ) ) override;
+ bool execute() override;
+
+ void getColumnNames( string_array & columnNames ) const override;
/**
* @brief Run test using loading protocol in table
- * @param i Relperm constitutive model
+ * @param relperm Relperm constitutive model
* @param table Table with input/output time history
*/
template< typename RELPERM_TYPE >
- std::enable_if_t< std::is_same< constitutive::TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void >
- runTest( RELPERM_TYPE & relperm,
- const arrayView2d< real64, 1 > & table );
-
- template< typename RELPERM_TYPE >
- std::enable_if_t< !std::is_same< constitutive::TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void >
- runTest( RELPERM_TYPE & relperm,
- const arrayView2d< real64, 1 > & table );
+ void runTest( RELPERM_TYPE & relperm, const arrayView2d< real64, 1 > & table );
+private:
/**
- * @brief Ouput table to file for easy plotting
+ * @brief Get the relative permeability model from the catalog
*/
- void outputResults();
+ constitutive::RelativePermeabilityBase & getRelperm();
+ constitutive::RelativePermeabilityBase const & getRelperm() const;
+
+ void allocateTable( integer numColumns, integer numPhases );
/**
- * @brief Read in a baseline table from file and compare with computed one (for unit testing purposes)
+ * @brief Initialises the table by filling in primary variables
*/
- void compareWithBaseline();
-
-private:
-
- template< typename RELPERM_TYPE >
- void resizeTables();
-
- template< typename RELPERM_TYPE >
- std::enable_if_t< std::is_same< constitutive::TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void >
- resizeTable();
-
- template< typename RELPERM_TYPE >
- std::enable_if_t< !std::is_same< constitutive::TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void >
- resizeTable();
+ void initializeTable();
/**
* @struct viewKeyStruct holds char strings and viewKeys for fast lookup
*/
struct viewKeyStruct
{
- constexpr static char const * relpermNameString()
- { return "relperm"; }
-
- constexpr static char const * numStepsString()
- { return "steps"; }
-
- constexpr static char const * outputString()
- { return "output"; }
-
- constexpr static char const * baselineString()
- { return "baseline"; }
+ constexpr static char const * relpermNameString() { return "relperm"; }
+ constexpr static char const * historicalSaturationsString() { return "historicalSaturations"; }
};
- integer m_numSteps; ///< Number of load steps
- integer m_numColumns; ///< Number of columns in data table (depends on number of fluid phases)
- integer m_numPhases; ///< Number of fluid phases
-
string m_relpermName; ///< relPermType identifier
- string m_outputFile; ///< Output file (optional, no output if not specified)
-
- array2d< real64 > m_table; ///< Table storing time-history of input/output
-
- Path m_baselineFile; ///< Baseline file (optional, for unit testing of solid models)
-
- enum columnKeys
- {
- TIME
- }; ///< Enumeration of "input" column keys for readability
-
- static constexpr real64 m_baselineTol = 1e-3; ///< Comparison tolerance for baseline results
+ array1d< real64 > m_historicalSaturations;
};
-
}
-#endif //GEOS_RELPERMDRIVER_HPP_
+#endif //GEOS_CONSTITUTIVEDRIVERS_RELATIVEPERMEABILITY_RELPERMDRIVER_HPP
diff --git a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverRunTest.hpp b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverRunTest.hpp
index eabe5616eee..e430e8fab3f 100644
--- a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverRunTest.hpp
+++ b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverRunTest.hpp
@@ -13,239 +13,119 @@
* ------------------------------------------------------------------------------------------------------------
*/
-#ifndef GEOS_RELPERMDRIVERRUNTEST_HPP_
-#define GEOS_RELPERMDRIVERRUNTEST_HPP_
+#ifndef GEOS_CONSTITUTIVEDRIVERS_RELATIVEPERMEABILITY_RELPERMDRIVERRUNTEST_HPP
+#define GEOS_CONSTITUTIVEDRIVERS_RELATIVEPERMEABILITY_RELPERMDRIVERRUNTEST_HPP
#include "constitutiveDrivers/relativePermeability/RelpermDriver.hpp"
+#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp"
#include "constitutive/relativePermeability/RelativePermeabilityFields.hpp"
-#include "constitutive/relativePermeability/Layouts.hpp"
#include "constitutive/relativePermeability/KilloughHysteresis.hpp"
-
+#include "common/DataLayouts.hpp"
namespace geos
{
-//specific to Hysteresis
+// Hysteresis traits
+template< typename RELPERM_TYPE >
+struct HasHysteresis : std::false_type
+{};
+
template< typename RELPERM_TYPE >
-std::enable_if_t< std::is_same< constitutive::TableRelativePermeabilityHysteresis,
- RELPERM_TYPE >::value, void >
-RelpermDriver::runTest( RELPERM_TYPE & relperm,
- const arrayView2d< real64 > & table )
+void
+RelpermDriver::runTest( RELPERM_TYPE & relperm, const arrayView2d< real64 > & table )
{
- // get number of phases and components
+ // Get the number of phases
integer const numPhases = relperm.numFluidPhases();
- // create kernel wrapper
+ // Create the kernel wrapper
+ typename RELPERM_TYPE::KernelWrapper const kernelWrapper = relperm.createKernelWrapper();
- typename constitutive::TableRelativePermeabilityHysteresis::KernelWrapper const kernelWrapper = relperm.createKernelWrapper();
+ // Offset for saturations in table
+ constexpr integer SATURATION = 1;
- // set saturation to user specified feed
- // it is more convenient to provide input in molar, so perform molar to mass conversion here
+ // Offset for relative permeability data
+ integer RELPERM = SATURATION + numPhases;
- array2d< real64, compflow::LAYOUT_PHASE > saturationValues;
- if( numPhases > 2 )
- {
- saturationValues.resize( ( m_numSteps + 1 ) * ( m_numSteps + 1 ), numPhases );
- }
- else
- {
- saturationValues.resize( m_numSteps + 1, numPhases );
- }
- using PT = typename RELPERM_TYPE::PhaseType;
- integer const ipWater = relperm.getPhaseOrder()[PT::WATER];
- integer const ipOil = relperm.getPhaseOrder()[PT::OIL];
- integer const ipGas = relperm.getPhaseOrder()[PT::GAS];
- localIndex offset = std::max( std::max( ipOil, ipWater ), std::max( ipOil, ipGas ) ) + 1;
-
- integer ipWetting = -1, ipNonWetting = -1;
- std::tie( ipWetting, ipNonWetting ) = relperm.wettingAndNonWettingPhaseIndices();
+ // Number of "cells"
+ integer const numRows = m_table.size( 0 );
- for( integer n = 0; n < table.size( 0 ); ++n )
+ // If we have hysteresis, we need to populate the historical saturations
+ if constexpr (HasHysteresis< RELPERM_TYPE >::value)
{
- if( m_numPhases > 2 )
- {
- saturationValues[n][ipWater] = table( n, ipWater + 1 );
- saturationValues[n][ipOil] = table( n, ipOil + 1 );
- saturationValues[n][ipGas] = table( n, ipGas + 1 );
- }
- else//two-phase
- {
- if( ipWater < 0 )
- {
- saturationValues[n][ipOil] = table( n, ipOil + 1 );
- saturationValues[n][ipGas] = table( n, ipGas + 1 );
- }
- else if( ipGas < 0 )
- {
- saturationValues[n][ipWater] = table( n, ipWater + 1 );
- saturationValues[n][ipOil] = table( n, ipOil + 1 );
- }
- else if( ipOil < 0 )
- {
- saturationValues[n][ipWater] = table( n, ipWater + 1 );
- saturationValues[n][ipGas] = table( n, ipGas + 1 );
- }
- }
- }
-
-
- arrayView2d< real64 const, compflow::USD_PHASE > const saturation = saturationValues.toViewConst();
+ // Shift the relative permeability offset
+ RELPERM += numPhases;
- auto const & phaseHasHysteresis = relperm.template getReference< array1d< integer > >( constitutive::TableRelativePermeabilityHysteresis::viewKeyStruct::phaseHasHysteresisString());
+ // Offset for the historical saturations
+ integer const OFFSET = SATURATION + numPhases;
- arrayView2d< real64, compflow::USD_PHASE > phaseMaxHistoricalVolFraction = relperm.template getField< fields::relperm::phaseMaxHistoricalVolFraction >().reference();
- arrayView2d< real64, compflow::USD_PHASE > phaseMinHistoricalVolFraction = relperm.template getField< fields::relperm::phaseMinHistoricalVolFraction >().reference();
+ using CurveType = constitutive::KilloughHysteresis::HysteresisCurve;
+ using KeyStruct = typename RELPERM_TYPE::viewKeyStruct;
-// arrayView1d< real64 > const drainagePhaseMinVolFraction = relperm.template getReference< array1d< real64 > >(
-// constitutive::TableRelativePermeabilityHysteresis::viewKeyStruct::drainagePhaseMinVolumeFractionString());
-// arrayView1d< real64 > const drainagePhaseMaxVolFraction = relperm.template getReference< array1d< real64 > >(
-// constitutive::TableRelativePermeabilityHysteresis::viewKeyStruct::drainagePhaseMaxVolumeFractionString());
- constitutive::KilloughHysteresis::HysteresisCurve const wettingCurve = relperm.template getReference< constitutive::KilloughHysteresis::HysteresisCurve >(
- constitutive::TableRelativePermeabilityHysteresis::viewKeyStruct::wettingCurveString());
+ auto [ipWetting, ipNonWetting] = relperm.wettingAndNonWettingPhaseIndices();
- constitutive::KilloughHysteresis::HysteresisCurve const nonWettingCurve = relperm.template getReference< constitutive::KilloughHysteresis::HysteresisCurve >(
- constitutive::TableRelativePermeabilityHysteresis::viewKeyStruct::nonWettingCurveString());
- //setting for drainage
- {
- if( phaseHasHysteresis[ipNonWetting] )
- {
- phaseMaxHistoricalVolFraction[0][ipNonWetting] = nonWettingCurve.m_extremumPhaseVolFraction;
- GEOS_LOG( GEOS_FMT( "New max non-wetting phase historical phase volume fraction: {}", phaseMaxHistoricalVolFraction[0][ipNonWetting] ) );
- }
- if( phaseHasHysteresis[ipWetting] )
- {
- phaseMinHistoricalVolFraction[0][ipWetting] = wettingCurve.m_extremumPhaseVolFraction;
- GEOS_LOG( GEOS_FMT( "New min wetting phase historical phase volume fraction: {}", phaseMinHistoricalVolFraction[0][ipWetting] ) );
- }
- }
+ stackArray1d< real64, constitutive::RelativePermeabilityBase::MAX_NUM_PHASES > historicalSaturations( numPhases );
+ historicalSaturations.zero();
-
-
- forAll< parallelDevicePolicy<> >( saturation.size( 0 ),
- [numPhases, kernelWrapper, saturation, table,
- offset] GEOS_HOST_DEVICE ( integer const n )
- {
- // nw phase set max to snw_max to get the imbibition bounding curve
- kernelWrapper.update( 0, 0, saturation[n] );
- for( integer p = 0; p < numPhases; ++p )
+ if( m_historicalSaturations.empty())
{
- table( n, offset + 1 + p ) = kernelWrapper.relperm()( 0, 0, p );
+ auto const & phaseHasHysteresis = relperm.template getReference< array1d< integer > >( KeyStruct::phaseHasHysteresisString());
+ if( phaseHasHysteresis[ipNonWetting] )
+ {
+ CurveType const nonWettingCurve = relperm.template getReference< CurveType >( KeyStruct::nonWettingCurveString() );
+ historicalSaturations[ipNonWetting] = nonWettingCurve.m_extremumPhaseVolFraction;
+ }
+ if( phaseHasHysteresis[ipWetting] )
+ {
+ CurveType const wettingCurve = relperm.template getReference< CurveType >( KeyStruct::wettingCurveString());
+ historicalSaturations[ipWetting] = wettingCurve.m_extremumPhaseVolFraction;
+ }
}
- } );
-
- //loop in charge of hysteresis values
- offset += numPhases;
-
-//setting for imbibition
- {
- if( phaseHasHysteresis[ipNonWetting] )
+ else
{
-
- phaseMaxHistoricalVolFraction[0][ipNonWetting] = nonWettingCurve.m_criticalDrainagePhaseVolFraction;
- GEOS_LOG( GEOS_FMT( "New max non-wetting phase historical phase volume fraction: {}", phaseMaxHistoricalVolFraction[0][ipNonWetting] ) );
+ for( integer p = 0; p < numPhases; ++p )
+ {
+ historicalSaturations[p] = m_historicalSaturations[p];
+ }
}
- if( phaseHasHysteresis[ipWetting] )
+ auto historicalView = historicalSaturations.toView();
+
+ arrayView2d< real64, compflow::USD_PHASE > phaseMaxHistoricalVolFraction = relperm.template getField< fields::relperm::phaseMaxHistoricalVolFraction >().reference();
+ arrayView2d< real64, compflow::USD_PHASE > phaseMinHistoricalVolFraction = relperm.template getField< fields::relperm::phaseMinHistoricalVolFraction >().reference();
+ forAll< parallelDevicePolicy<> >( numRows,
+ [numPhases, table,
+ ipWetting, ipNonWetting,
+ OFFSET,
+ phaseMaxHistoricalVolFraction,
+ phaseMinHistoricalVolFraction,
+ historicalView] GEOS_HOST_DEVICE ( integer const n )
{
- phaseMinHistoricalVolFraction[0][ipWetting] = wettingCurve.m_criticalDrainagePhaseVolFraction;
- GEOS_LOG( GEOS_FMT( "New min wetting phase historical phase volume fraction: {}", phaseMinHistoricalVolFraction[0][ipWetting] ) );
- }
+ phaseMaxHistoricalVolFraction( n, ipNonWetting ) = historicalView[ipNonWetting];
+ phaseMinHistoricalVolFraction( n, ipWetting ) = historicalView[ipWetting];
+ for( integer p = 0; p < numPhases; ++p )
+ {
+ table( n, p + OFFSET ) = historicalView[p];
+ }
+ } );
}
-
-
- forAll< parallelDevicePolicy<> >( saturation.size( 0 ),
- [numPhases, kernelWrapper, saturation, table,
- offset] GEOS_HOST_DEVICE ( integer const n )
+ forAll< parallelDevicePolicy<> >( numRows,
+ [numPhases, kernelWrapper, table,
+ RELPERM] GEOS_HOST_DEVICE ( integer const n )
{
- // nw phase set max to snw_max to get the imbibition bounding curve
+ StackArray< real64, 2, constitutive::RelativePermeabilityBase::MAX_NUM_PHASES, compflow::LAYOUT_PHASE > saturation( 1, numPhases );
- kernelWrapper.update( 0, 0, saturation[n] );
for( integer p = 0; p < numPhases; ++p )
{
- table( n, offset + 1 + p ) = kernelWrapper.relperm()( 0, 0, p );
+ saturation[0][p] = table( n, SATURATION + p );
}
- } );
-
-
-}
-
-template< typename RELPERM_TYPE >
-std::enable_if_t< !std::is_same< constitutive::TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void >
-RelpermDriver::runTest( RELPERM_TYPE & relperm,
- const arrayView2d< real64 > & table )
-{
- // get number of phases and components
-
- integer const numPhases = relperm.numFluidPhases();
-
- // create kernel wrapper
-
- typename RELPERM_TYPE::KernelWrapper const kernelWrapper = relperm.createKernelWrapper();
-
- // set saturation to user specified feed
- // it is more convenient to provide input in molar, so perform molar to mass conversion here
-
- array2d< real64, compflow::LAYOUT_PHASE > saturationValues;
- if( numPhases > 2 )
- {
- saturationValues.resize(( m_numSteps + 1 ) * ( m_numSteps + 1 ), numPhases );
- }
- else
- {
- saturationValues.resize( m_numSteps + 1, numPhases );
- }
- using PT = typename RELPERM_TYPE::PhaseType;
- integer const ipWater = relperm.getPhaseOrder()[PT::WATER];
- integer const ipOil = relperm.getPhaseOrder()[PT::OIL];
- integer const ipGas = relperm.getPhaseOrder()[PT::GAS];
- const localIndex offset = std::max( std::max( ipOil, ipWater ), std::max( ipOil, ipGas ) ) + 1;
-
- for( integer n = 0; n < table.size( 0 ); ++n )
- {
-
-
- if( m_numPhases > 2 )
- {
- saturationValues[n][ipWater] = table( n, ipWater + 1 );
- saturationValues[n][ipOil] = table( n, ipOil + 1 );
- saturationValues[n][ipGas] = table( n, ipGas + 1 );
- }
- else//two-phase
- {
- if( ipWater < 0 )
- {
- saturationValues[n][ipOil] = table( n, ipOil + 1 );
- saturationValues[n][ipGas] = table( n, ipGas + 1 );
- }
- else if( ipGas < 0 )
- {
- saturationValues[n][ipWater] = table( n, ipWater + 1 );
- saturationValues[n][ipOil] = table( n, ipOil + 1 );
- }
- }
-
- }
-
- arrayView2d< real64 const, compflow::USD_PHASE > const saturation = saturationValues.toViewConst();
-
- // perform relperm update using table (Swet,Snonwet) and save resulting total density, etc.
- // note: column indexing should be kept consistent with output file header below.
-
- forAll< parallelDevicePolicy<> >( saturation.size( 0 ),
- [numPhases, kernelWrapper, saturation, table,
- offset] GEOS_HOST_DEVICE ( integer const n )
- {
- kernelWrapper.update( 0, 0, saturation[n] );
+ kernelWrapper.update( n, 0, saturation[0] );
for( integer p = 0; p < numPhases; ++p )
{
- table( n, offset + 1 + p ) = kernelWrapper.relperm()( 0, 0, p );
+ table( n, RELPERM + p ) = kernelWrapper.relperm()( n, 0, p );
}
} );
-
}
-
}
-
-#endif //GEOS_RELPERMDRIVERRUNTEST_HPP_
+#endif //GEOS_CONSTITUTIVEDRIVERS_RELATIVEPERMEABILITY_RELPERMDRIVERRUNTEST_HPP
diff --git a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverTableRelativeHysteresisRunTest.cpp b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverTableRelativeHysteresisRunTest.cpp
index 4ed6420b673..d50bd1bdf19 100644
--- a/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverTableRelativeHysteresisRunTest.cpp
+++ b/src/coreComponents/constitutiveDrivers/relativePermeability/RelpermDriverTableRelativeHysteresisRunTest.cpp
@@ -16,8 +16,12 @@
#include "RelpermDriverRunTest.hpp"
#include "constitutive/relativePermeability/TableRelativePermeabilityHysteresis.hpp"
-
namespace geos
{
+
+template<>
+struct HasHysteresis< constitutive::TableRelativePermeabilityHysteresis > : std::true_type
+{};
+
template void RelpermDriver::runTest< geos::constitutive::TableRelativePermeabilityHysteresis >( geos::constitutive::TableRelativePermeabilityHysteresis &, arrayView2d< real64 > const & );
}
diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index 9a9bb51d843..3449d34d4b7 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -6484,18 +6484,20 @@ Information output from lower logLevels is added with the desired log level
-
-
+
+
-
-
+
+
+
+
-
+
diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other
index c7434b7b99c..15258cad63f 100644
--- a/src/coreComponents/schema/schema.xsd.other
+++ b/src/coreComponents/schema/schema.xsd.other
@@ -1554,7 +1554,10 @@ A field can represent a physical variable. (pressure, temperature, global compos
-
+
+
+
+