diff --git a/src/coreComponents/fieldSpecification/AquiferBoundaryCondition.cpp b/src/coreComponents/fieldSpecification/AquiferBoundaryCondition.cpp index 4a9c99a7b2f..391171c0872 100644 --- a/src/coreComponents/fieldSpecification/AquiferBoundaryCondition.cpp +++ b/src/coreComponents/fieldSpecification/AquiferBoundaryCondition.cpp @@ -114,6 +114,8 @@ AquiferBoundaryCondition::AquiferBoundaryCondition( string const & name, Group * void AquiferBoundaryCondition::postInputInitialization() { + FieldSpecificationBase::postInputInitialization(); + GEOS_THROW_IF_LE_MSG( m_permeability, 0.0, "The aquifer permeability cannot be equal to zero or negative", InputError, getDataContext() ); diff --git a/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp b/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp index 6e0f1181e6e..21e12864c45 100644 --- a/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp +++ b/src/coreComponents/fieldSpecification/EquilibriumInitialCondition.cpp @@ -96,6 +96,7 @@ EquilibriumInitialCondition::EquilibriumInitialCondition( string const & name, G void EquilibriumInitialCondition::postInputInitialization() { + FieldSpecificationBase::postInputInitialization(); FunctionManager const & functionManager = FunctionManager::getInstance(); diff --git a/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp b/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp index 94e78e7e06a..210defcd825 100644 --- a/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp +++ b/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp @@ -15,6 +15,8 @@ #include "FieldSpecificationBase.hpp" +#include "common/format/StringUtilities.hpp" +#include "common/logger/Logger.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" namespace geos @@ -39,6 +41,11 @@ FieldSpecificationBase::FieldSpecificationBase( string const & name, Group * par setInputFlag( InputFlags::OPTIONAL ). setDescription( "Path to the target field" ); + registerWrapper( viewKeyStruct::regionNamesString(), &m_regionNames ). + setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Names of the regions where the field specification is applied." ); + registerWrapper( viewKeyStruct::fieldNameString(), &m_fieldName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setInputFlag( InputFlags::OPTIONAL ). @@ -61,6 +68,13 @@ FieldSpecificationBase::FieldSpecificationBase( string const & name, Group * par setInputFlag( InputFlags::OPTIONAL ). setDescription( "Name of function that specifies variation of the boundary condition." ); + registerWrapper( viewKeyStruct::functionNamesString(), &m_functionNames ). + setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). + setInputFlag( InputFlags::OPTIONAL ). + setSizedFromParent( 0 ). + setDescription( "Names of per-component functions that specifies variation of the boundary condition.\n" + "Either left empty or sized exactly like 'scales'.\n" ); + registerWrapper( viewKeyStruct::bcApplicationTableNameString(), &m_bcApplicationFunctionName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setInputFlag( InputFlags::OPTIONAL ). @@ -71,6 +85,11 @@ FieldSpecificationBase::FieldSpecificationBase( string const & name, Group * par setInputFlag( InputFlags::OPTIONAL ). setDescription( "Apply a scaling factor for the value of the boundary condition." ); + registerWrapper( viewKeyStruct::scalesString(), &m_scales ). + setInputFlag( InputFlags::OPTIONAL ). + setSizedFromParent( 0 ). + setDescription( "Apply scaling factors for the values of every component of the boundary condition.\n" ); + registerWrapper( viewKeyStruct::initialConditionString(), &m_initialCondition ). setApplyDefaultValue( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -110,22 +129,66 @@ FieldSpecificationBase::getCatalog() } +void FieldSpecificationBase::postInputInitialization() +{ + GEOS_THROW_IF( !m_functionNames.empty() && + m_functionNames.size() != static_cast< string_array::size_type >( m_scales.size() ), + GEOS_FMT ( "Size mismatch: '{}' has {} entries but '{}' has {}. " + "Either leave '{}' empty or size it exactly like '{}'", + viewKeyStruct::functionNamesString(), m_functionNames.size(), + viewKeyStruct::scalesString(), m_scales.size(), + viewKeyStruct::functionNamesString(), viewKeyStruct::scalesString() ), + InputError, + getDataContext() ); + + if( usesNonScalarValues() ) + { + GEOS_THROW_IF( m_component != -1, + GEOS_FMT ( "'{}' must not be set when '{}' is set.", + viewKeyStruct::componentString(), + viewKeyStruct::scalesString() ), + InputError, + getDataContext() ); + + GEOS_THROW_IF( !m_functionName.empty(), + GEOS_FMT ( "'{}' must not be set when '{}' is set." + "Use '{}' to provide one function per component instead", + viewKeyStruct::functionNameString(), + viewKeyStruct::scalesString(), + viewKeyStruct::functionNamesString() ), + InputError, + getDataContext() ); + } + + if( !m_regionNames.empty() ) + { + GEOS_THROW_IF( !m_objectPath.empty(), + GEOS_FMT ( "'{}' must not be set when '{}' is set.", + viewKeyStruct::objectPathString(), + viewKeyStruct::regionNamesString() ), + InputError, + getDataContext() ); + } +} void FieldSpecificationBase::setMeshObjectPath( Group const & meshBodies ) { + string const path = m_regionNames.empty() + ? m_objectPath + : "ElementRegions/{" + stringutilities::join( m_regionNames, ' ' ) + "}"; try { - m_meshObjectPaths = std::make_unique< MeshObjectPath >( m_objectPath, meshBodies ); + m_meshObjectPaths = std::make_unique< MeshObjectPath >( path, meshBodies ); } catch( std::exception const & e ) { ErrorLogger::global().modifyCurrentExceptionMessage() .addToMsg( getWrapperDataContext( viewKeyStruct::objectPathString() ).toString() + - " is a wrong objectPath: " + m_objectPath + "\n" ) + " is a wrong objectPath: " + path + "\n" ) .addContextInfo( getWrapperDataContext( viewKeyStruct::objectPathString() ).getContextInfo() .setPriority( 2 ) ); throw InputError( e, getWrapperDataContext( viewKeyStruct::objectPathString() ).toString() + - " is a wrong objectPath: " + m_objectPath + "\n" ); + " is a wrong objectPath: " + path + "\n" ); } } diff --git a/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp b/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp index 1905a4e5c36..c0206a7b456 100644 --- a/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp +++ b/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp @@ -355,6 +355,54 @@ class FieldSpecificationBase : public dataRepository::Group arrayView1d< real64 > const & rhsContribution, LAMBDA && lambda ) const; + /** + * @brief Compute the contributions that will be added/enforced to the right-hand side, + * and collect the corresponding dof numbers + * @tparam FIELD_OP A wrapper struct to define how the boundary condition operates on the variables. + * Either \ref OpEqual or \ref OpAdd. + * @tparam POLICY Execution policy to use when iterating over target set. + * @tparam LAMBDA The type of lambda function passed into the parameter list. + * @param[in] component The field component to apply the boundary condition to. + * @param[in] scale The scale factor to apply to the boundary condition value. + * @param[in] functionName The name of the function used to evaluate the boundary condition value. + * @param[in] targetSet The set of indices which the boundary condition will be applied. + * @param[in] time The time at which any time dependent functions are to be evaluated as part of the + * application of the boundary condition. + * @param[in] dt time step size which is applied as a factor to bc values + * @param[in] dataGroup The Group that contains the field to apply the boundary condition to. + * @param[in] dofMap The map from the local index of the primary field to the global degree of + * freedom number. + * @param[in] dofRankOffset Offset of dof indices on current rank. + * @param[inout] matrix Local part of the system matrix. + * @param[inout] dof array storing the degrees of freedom of the rhsContribution, to know where + * in the rhs they will be added/enforced + * @param[inout] rhsContribution array storing the values that will be added/enforced to the right-hand side + * @param[in] lambda A lambda function which defines how the value that is passed into the functions + * provided by the FIELD_OP templated type. + * + * This overload behaves like the legacy computeRhsContribution, but takes the component, scale + * and functionName as explicit arguments rather than reading them from the corresponding members. + * This is the variant called when using non-scalar valued field specifications/boundary conditions, + * where each component is applied in turn. + * + * Note that this function only computes the rhs contributions, but does not apply them to the right-hand side. + * The application of these rhs contributions is done in applyBoundaryConditionToSystem. + */ + template< typename FIELD_OP, typename POLICY, typename LAMBDA > + void + computeRhsContribution( integer const component, + real64 const scale, + string const & functionName, + SortedArrayView< localIndex const > const & targetSet, + real64 const time, + real64 const dt, + dataRepository::Group const & dataGroup, + arrayView1d< globalIndex const > const & dofMap, + globalIndex const dofRankOffset, + CRSMatrixView< real64, globalIndex const > const & matrix, + arrayView1d< globalIndex > const & dof, + arrayView1d< real64 > const & rhsContribution, + LAMBDA && lambda ) const; /** * @brief Function to zero matrix rows to apply boundary conditions @@ -371,6 +419,7 @@ class FieldSpecificationBase : public dataRepository::Group arrayView1d< globalIndex const > const & dofMap, CRSMatrixView< real64, globalIndex const > const & matrix ) const; + /** * @brief View keys */ @@ -382,6 +431,8 @@ class FieldSpecificationBase : public dataRepository::Group constexpr static char const * constitutivePathString() { return "constitutivePath"; } /// @return The key for objectPath constexpr static char const * objectPathString() { return "objectPath"; } + /// @return The key for regionNames + constexpr static char const * regionNamesString() { return "regionNames"; } /// @return The key for fieldName constexpr static char const * fieldNameString() { return "fieldName"; } /// @return The key for dataType @@ -394,8 +445,12 @@ class FieldSpecificationBase : public dataRepository::Group constexpr static char const * bcApplicationTableNameString() { return "bcApplicationTableName"; } /// @return The key for scale constexpr static char const * scaleString() { return "scale"; } + /// @return The key for scales + constexpr static char const * scalesString() { return "scales"; } /// @return The key for functionName constexpr static char const * functionNameString() { return "functionName"; } + /// @return The key for functionNames + constexpr static char const * functionNamesString() { return "functionNames"; } /// @return The key for initialCondition constexpr static char const * initialConditionString() { return "initialCondition"; } /// @return The key for beginTime @@ -415,6 +470,15 @@ class FieldSpecificationBase : public dataRepository::Group return m_functionName; } + /** + * Accessor + * @return const reference to m_functionNames + */ + string_array const & getFunctionNames() const + { + return m_functionNames; + } + /** * Accessor * @return const reference to m_objectPath @@ -424,6 +488,15 @@ class FieldSpecificationBase : public dataRepository::Group return m_objectPath; } + /** + * Accessor + * @return const reference to m_regionNames + */ + string_array const & getRegionNames() const + { + return m_regionNames; + } + /** * Accessor * @return const reference to m_fieldName @@ -497,6 +570,15 @@ class FieldSpecificationBase : public dataRepository::Group return m_scale; } + /** + * Accessor + * @return const m_scales + */ + arrayView1d< real64 const > getScales() const + { + return m_scales.toViewConst(); + } + /** * Mutator * @param[in] fieldName The name of the field @@ -515,6 +597,15 @@ class FieldSpecificationBase : public dataRepository::Group m_objectPath = objectPath; } + /** + * Mutator + * @param[in] regionNames The region names + */ + void setRegionNames( string_array const & regionNames ) + { + m_regionNames = regionNames; + } + /** * Mutator * @param[in] scale Scaling factor @@ -524,6 +615,27 @@ class FieldSpecificationBase : public dataRepository::Group m_scale = scale; } + /** + * Mutator + * @brief Set the per-component scale factors + * @param[in] scales The tensor-valued scale + */ + void setScales( array1d< real64 > const & scales ) + { + m_scales = scales; + } + + /** + * Mutator + * @brief Set the per-component function names + * @param[in] functionNames The per-component function names. Must have the same + * size as @p m_scales or be empty. + */ + void setFunctionNames( string_array const & functionNames ) + { + m_functionNames = functionNames; + } + /** * Mutator * @param[in] isInitialCondition Logical value to indicate if it is an initial condition @@ -559,9 +671,38 @@ class FieldSpecificationBase : public dataRepository::Group return *(m_meshObjectPaths.get()); } + /** + * @brief Query whether this field specification uses non-scalar scales + */ + bool usesNonScalarValues() const + { + return !m_scales.empty() || !m_functionNames.empty(); + } + + template< typename LAMBDA > + void forEachComponent( LAMBDA && lambda ) const + { + if( usesNonScalarValues() ) + { + localIndex const numComponents = m_scales.size(); + for( localIndex comp = 0; comp < numComponents; ++comp ) + { + string const & functionName = m_functionName.empty() ? string{} + : m_functionNames[ comp ]; + lambda( comp, m_scales[ comp ], functionName ); + } + } + else + { + lambda( getComponent(), m_scale, m_functionName ); + } + } + protected: + virtual void postInputInitialization() override; + private: @@ -574,6 +715,9 @@ class FieldSpecificationBase : public dataRepository::Group std::unique_ptr< MeshObjectPath > m_meshObjectPaths; + /// the region names where the field specification is applied + string_array m_regionNames; + /// the name of the field the boundary condition is applied to or a key string to use for /// determining whether or not to apply the boundary condition. string m_fieldName; @@ -591,9 +735,16 @@ class FieldSpecificationBase : public dataRepository::Group /// The name of the function used to generate values for application. string m_functionName; + /// Per-component function names used when @p m_scales.size() > 1 + /// Either empty or sized exactly like @p m_scales + string_array m_functionNames; + /// The scale factor to use on the value of the boundary condition. real64 m_scale; + /// Per-component scale factors + array1d< real64 > m_scales; + /// Time after which the bc is allowed to be applied real64 m_beginTime; @@ -614,25 +765,28 @@ void FieldSpecificationBase::applyFieldValueKernel( ArrayView< T, N, USD > const real64 const time, Group & dataGroup ) const { - integer const component = getComponent(); FunctionManager & functionManager = FunctionManager::getInstance(); - if( m_functionName.empty() ) + forEachComponent( [&]( integer const component, + real64 const scale, + string const & functionName ) { - real64 const value = m_scale; - forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + if( functionName.empty() ) { - localIndex const a = targetSet[ i ]; - FIELD_OP::SpecifyFieldValue( field, a, component, value ); - } ); - } - else - { + real64 const value = scale; + forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + localIndex const a = targetSet[ i ]; + FIELD_OP::SpecifyFieldValue( field, a, component, value ); + } ); + return; + } + FunctionBase const & function = [&]() -> FunctionBase const & { try { - return functionManager.getGroup< FunctionBase >( m_functionName ); + return functionManager.getGroup< FunctionBase >( functionName ); } catch( std::exception const & e ) { @@ -648,7 +802,7 @@ void FieldSpecificationBase::applyFieldValueKernel( ArrayView< T, N, USD > const if( function.isFunctionOfTime()==2 ) { - real64 const value = m_scale * function.evaluate( &time ); + real64 const value = scale * function.evaluate( &time ); forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { localIndex const a = targetSet[ i ]; @@ -660,14 +814,14 @@ void FieldSpecificationBase::applyFieldValueKernel( ArrayView< T, N, USD > const real64_array result( static_cast< localIndex >( targetSet.size() ) ); function.evaluate( dataGroup, time, targetSet, result ); arrayView1d< real64 const > const & resultView = result.toViewConst(); - real64 const scale = m_scale; + real64 const localScale = scale; forAll< POLICY >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { localIndex const a = targetSet[ i ]; - FIELD_OP::SpecifyFieldValue( field, a, component, scale * resultView[i] ); + FIELD_OP::SpecifyFieldValue( field, a, component, localScale * resultView[i] ); } ); } - } + } ); } @@ -703,13 +857,38 @@ void FieldSpecificationBase::applyBoundaryConditionToSystemKernel( SortedArrayVi arrayView1d< real64 > const & rhs, ArrayView< T const, NDIM, USD > const & fieldView ) const { - integer const component = getComponent(); - this->applyBoundaryConditionToSystem< FIELD_OP, POLICY >( targetSet, time, dataGroup, dofMap, dofRankOffset, matrix, rhs, - [fieldView, component] GEOS_HOST_DEVICE ( localIndex const a ) + forEachComponent( [&]( integer const component, + real64 const scale, + string const & functionName ) { - real64 value = 0.0; - FieldSpecificationEqual::readFieldValue( fieldView, a, component, value ); - return value; + integer const comp = ( component >= 0 ) ? component : 0; + + array1d< globalIndex > dofArray( targetSet.size() ); + arrayView1d< globalIndex > const & dof = dofArray.toView(); + + array1d< real64 > rhsContributionArray( targetSet.size() ); + arrayView1d< real64 > const & rhsContribution = rhsContributionArray.toView(); + + computeRhsContribution< FIELD_OP, POLICY >( comp, + scale, + functionName, + targetSet, + time, + 1.0, + dataGroup, + dofMap, + dofRankOffset, + matrix, + dof, + rhsContribution, + [fieldView, component] GEOS_HOST_DEVICE ( localIndex const a ) + { + real64 value = 0.0; + FieldSpecificationEqual::readFieldValue( fieldView, a, component, value ); + return value; + } ); + + FIELD_OP::template prescribeRhsValues< POLICY >( rhs, dof, dofRankOffset, rhsContribution ); } ); } @@ -779,30 +958,43 @@ FieldSpecificationBase:: arrayView1d< real64 > const & rhs, LAMBDA && lambda ) const { - array1d< globalIndex > dofArray( targetSet.size() ); - arrayView1d< globalIndex > const & dof = dofArray.toView(); - - array1d< real64 > rhsContributionArray( targetSet.size() ); - arrayView1d< real64 > const & rhsContribution = rhsContributionArray.toView(); - - computeRhsContribution< FIELD_OP, POLICY, LAMBDA >( targetSet, - time, - dt, - dataGroup, - dofMap, - dofRankOffset, - matrix, - dof, - rhsContribution, - std::forward< LAMBDA >( lambda ) ); - - FIELD_OP::template prescribeRhsValues< POLICY >( rhs, dof, dofRankOffset, rhsContribution ); + forEachComponent( [&]( integer const component, + real64 const scale, + string const & functionName ) + { + integer const comp = ( component >= 0 ) ? component : 0; + + array1d< globalIndex > dofArray( targetSet.size() ); + arrayView1d< globalIndex > const & dof = dofArray.toView(); + + array1d< real64 > rhsContributionArray( targetSet.size() ); + arrayView1d< real64 > const & rhsContribution = rhsContributionArray.toView(); + + computeRhsContribution< FIELD_OP, POLICY >( comp, + scale, + functionName, + targetSet, + time, + dt, + dataGroup, + dofMap, + dofRankOffset, + matrix, + dof, + rhsContribution, + lambda ); + + FIELD_OP::template prescribeRhsValues< POLICY >( rhs, dof, dofRankOffset, rhsContribution ); + } ); } template< typename FIELD_OP, typename POLICY, typename LAMBDA > void FieldSpecificationBase:: - computeRhsContribution( SortedArrayView< localIndex const > const & targetSet, + computeRhsContribution( integer const component, + real64 const scale, + string const & functionName, + SortedArrayView< localIndex const > const & targetSet, real64 const time, real64 const dt, dataRepository::Group const & dataGroup, @@ -813,8 +1005,6 @@ FieldSpecificationBase:: arrayView1d< real64 > const & rhsContribution, LAMBDA && lambda ) const { - integer const component = ( getComponent() >=0 ) ? getComponent() : 0; - string const & functionName = getReference< string >( viewKeyStruct::functionNameString() ); FunctionManager & functionManager = FunctionManager::getInstance(); // Compute the value of the rhs terms, and collect the dof numbers @@ -822,7 +1012,7 @@ FieldSpecificationBase:: if( functionName.empty() || functionManager.getGroup< FunctionBase >( functionName ).isFunctionOfTime() == 2 ) { - real64 value = m_scale * dt; + real64 value = scale * dt; if( !functionName.empty() ) { FunctionBase const & function = functionManager.getGroup< FunctionBase >( functionName ); @@ -849,7 +1039,7 @@ FieldSpecificationBase:: real64_array resultsArray( targetSet.size() ); function.evaluate( dataGroup, time, targetSet, resultsArray ); arrayView1d< real64 const > const & results = resultsArray.toViewConst(); - real64 const value = m_scale * dt; + real64 const value = scale * dt; forAll< POLICY >( targetSet.size(), [targetSet, dof, dofMap, dofRankOffset, component, matrix, rhsContribution, results, value, lambda] GEOS_HOST_DEVICE ( @@ -867,6 +1057,36 @@ FieldSpecificationBase:: } } +template< typename FIELD_OP, typename POLICY, typename LAMBDA > +void +FieldSpecificationBase:: + computeRhsContribution( SortedArrayView< localIndex const > const & targetSet, + real64 const time, + real64 const dt, + dataRepository::Group const & dataGroup, + arrayView1d< globalIndex const > const & dofMap, + globalIndex const dofRankOffset, + CRSMatrixView< real64, globalIndex const > const & matrix, + arrayView1d< globalIndex > const & dof, + arrayView1d< real64 > const & rhsContribution, + LAMBDA && lambda ) const +{ + computeRhsContribution< FIELD_OP, POLICY, LAMBDA >( ( getComponent() >= 0 ) ? getComponent() : 0, + m_scale, + m_functionName, + targetSet, + time, + dt, + dataGroup, + dofMap, + dofRankOffset, + matrix, + dof, + rhsContribution, + std::forward< LAMBDA >( lambda ) ); + +} + template< typename POLICY > void FieldSpecificationBase::zeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet, @@ -874,19 +1094,25 @@ void FieldSpecificationBase::zeroSystemRowsForBoundaryCondition( SortedArrayView CRSMatrixView< real64, globalIndex const > const & matrix ) const { - integer const component = ( getComponent() >=0 ) ? getComponent() : 0; - forAll< POLICY >( targetSet.size(), [targetSet, dofMap, matrix, component] GEOS_HOST_DEVICE ( localIndex const i ) + forEachComponent( [&]( integer const component, + real64 const scale, + string const & functionName ) { - localIndex const a = targetSet[ i ]; - globalIndex const dof = dofMap[ a ] + component; + GEOS_UNUSED_VAR( scale, functionName ); + integer const comp = ( component >= 0 ) ? component : 0; + forAll< POLICY >( targetSet.size(), [targetSet, dofMap, matrix, component] GEOS_HOST_DEVICE ( localIndex const i ) + { + localIndex const a = targetSet[ i ]; + globalIndex const dof = dofMap[ a ] + component; - arraySlice1d< real64 > const entries = matrix.getEntries( dof ); - localIndex const numEntries = matrix.numNonZeros( dof ); + arraySlice1d< real64 > const entries = matrix.getEntries( dof ); + localIndex const numEntries = matrix.numNonZeros( dof ); - for( localIndex j = 0; j < numEntries; ++j ) - { - entries[ j ] = 0; - } + for( localIndex j = 0; j < numEntries; ++j ) + { + entries[ j ] = 0; + } + } ); } ); } diff --git a/src/coreComponents/fieldSpecification/PerfectlyMatchedLayer.cpp b/src/coreComponents/fieldSpecification/PerfectlyMatchedLayer.cpp index 7e71480536c..ae56811c914 100644 --- a/src/coreComponents/fieldSpecification/PerfectlyMatchedLayer.cpp +++ b/src/coreComponents/fieldSpecification/PerfectlyMatchedLayer.cpp @@ -78,6 +78,8 @@ PerfectlyMatchedLayer::PerfectlyMatchedLayer( string const & name, Group * const void PerfectlyMatchedLayer::postInputInitialization() { + FieldSpecificationBase::postInputInitialization(); + GEOS_THROW_IF( (m_xMax[0]( getDirection() ) < 1e-20,