diff --git a/src/coreComponents/common/FieldSpecificationOps.hpp b/src/coreComponents/common/FieldSpecificationOps.hpp index d663aae38a7..0a26cbea82b 100644 --- a/src/coreComponents/common/FieldSpecificationOps.hpp +++ b/src/coreComponents/common/FieldSpecificationOps.hpp @@ -67,6 +67,26 @@ struct OpAdd } }; +/** + * @brief OpMultiply Operator that multiplies by a value + */ +struct OpMultiply +{ + /** + * @brief Pointwise update of a value + * @tparam T type of the left-hand side + * @tparam U type of the right-hand side + * @param[in] lhs value to update + * @param[in] rhs input value + */ + template< typename T, typename U > + GEOS_HOST_DEVICE static inline + void apply( T & lhs, U const & rhs ) + { + lhs *= static_cast< T >( rhs ); + } +}; + /** * @brief FieldSpecificationOp */ @@ -681,6 +701,70 @@ struct FieldSpecificationAdd : public FieldSpecificationOp< OpAdd > }; + +/** + * @struct FieldSpecificationMultiply + * this struct a collection of static functions which adhere to an assumed interface for multiplying + * a value for a field. + */ +struct FieldSpecificationMultiply : public FieldSpecificationOp< OpMultiply > +{ + /// Alias for FieldSpecificationOp< OpMultiply > + using base_type = FieldSpecificationOp< OpMultiply >; + using base_type::SpecifyFieldValue; + + /** + * @brief Function to apply a value to a vector field for a single dof. + * @param[in] dof The degree of freedom that is to be modified. + * @param[in] dofRankOffset offset of dof indices on current rank + * @param[in] matrix A ParalleMatrix object: the system matrix. + * @param[out] rhs The rhs contribution to be modified + * @param[in] bcValue The value to multiply to rhs + * @param[in] fieldValue unused. + * + */ + GEOS_HOST_DEVICE + static inline void + SpecifyFieldValue( globalIndex const dof, + globalIndex const dofRankOffset, + CRSMatrixView< real64, globalIndex const > const & matrix, + real64 & rhs, + real64 const bcValue, + real64 const fieldValue ) + { + GEOS_UNUSED_VAR( dof ); + GEOS_UNUSED_VAR( dofRankOffset ); + GEOS_UNUSED_VAR( matrix ); + GEOS_UNUSED_VAR( fieldValue ); + rhs *= bcValue; + } + + /** + * @brief Function to add some values of a vector. + * @tparam POLICY the execution policy to use when setting values + * @param rhs the target right-hand side vector + * @param dof a list of global DOF indices to be set + * @param dofRankOffset offset of dof indices on current rank + * @param values a list of values corresponding to \p dof that will be added to \p rhs. + */ + template< typename POLICY > + static inline void prescribeRhsValues( arrayView1d< real64 > const & rhs, + arrayView1d< globalIndex const > const & dof, + globalIndex const dofRankOffset, + arrayView1d< real64 const > const & values ) + { + GEOS_ASSERT_EQ( dof.size(), values.size() ); + forAll< POLICY >( dof.size(), [rhs, dof, dofRankOffset, values] GEOS_HOST_DEVICE ( localIndex const a ) + { + globalIndex const localRow = dof[ a ] - dofRankOffset; + if( localRow >= 0 && localRow < rhs.size() ) + { rhs[ localRow ] *= values[ a ]; } + } ); + } + +}; + + } //namespace geos #endif //GEOS_COMMON_FIELDSPECIFICATIONOPS_HPP diff --git a/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp b/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp index 693ecc2f4aa..38d912afb47 100644 --- a/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp +++ b/src/coreComponents/fieldSpecification/FieldSpecificationBase.cpp @@ -83,6 +83,11 @@ FieldSpecificationBase::FieldSpecificationBase( string const & name, Group * par setInputFlag( InputFlags::OPTIONAL ). setDescription( "Time at which the boundary condition will stop being applied." ); + registerWrapper( viewKeyStruct::isScalingString(), &m_isScaling ). + setApplyDefaultValue( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Boundary condition is a multiplicative scaling of the values already present." ); + enableLogLevelInput(); } diff --git a/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp b/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp index cdce71bc664..9472dce97f8 100644 --- a/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp +++ b/src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp @@ -392,6 +392,8 @@ class FieldSpecificationBase : public dataRepository::Group constexpr static char const * beginTimeString() { return "beginTime"; } /// @return The key for endTime constexpr static char const * endTimeString() { return "endTime"; } + /// @return The key for isScaling + constexpr static char const * isScalingString() { return "isScaling"; } }; /** @@ -476,6 +478,15 @@ class FieldSpecificationBase : public dataRepository::Group return m_initialCondition; } + /** + * Accessor + * @return const m_isScaling + */ + int isScaling() const + { + return m_isScaling; + } + /** * Accessor * @return const m_scale @@ -591,6 +602,9 @@ class FieldSpecificationBase : public dataRepository::Group /// The name of a function used to turn on and off the boundary condition. string m_bcApplicationFunctionName; + /// Whether or not the boundary condition is a multiplicative scaling of what is already present (from mesh) + int m_isScaling; + }; diff --git a/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp b/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp index 155596304d6..93d2384060e 100644 --- a/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp +++ b/src/coreComponents/fieldSpecification/FieldSpecificationManager.cpp @@ -216,7 +216,7 @@ void FieldSpecificationManager::applyInitialConditions( MeshLevel & mesh ) const { this->forSubGroups< FieldSpecificationBase >( [&] ( FieldSpecificationBase const & fs ) { - if( fs.initialCondition() ) + if( fs.initialCondition() && !fs.isScaling() ) { fs.apply< dataRepository::Group >( mesh, [&]( FieldSpecificationBase const & bc, @@ -225,10 +225,28 @@ void FieldSpecificationManager::applyInitialConditions( MeshLevel & mesh ) const Group & targetGroup, string const fieldName ) { - bc.applyFieldValue< FieldSpecificationEqual >( targetSet, 0.0, targetGroup, fieldName ); + bc.applyFieldValue< FieldSpecificationEqual >( targetSet, 0.0, targetGroup, fieldName ); } ); } } ); } +void FieldSpecificationManager::applyScalingInitialConditions( MeshLevel & mesh ) const +{ + this->forSubGroups< FieldSpecificationBase >( [&] ( FieldSpecificationBase const & fs ) + { + if( fs.initialCondition() && fs.isScaling() ) + { + fs.apply< dataRepository::Group >( mesh, + [&]( FieldSpecificationBase const & bc, + string const &, + SortedArrayView< localIndex const > const & targetSet, + Group & targetGroup, + string const fieldName ) + { + bc.applyFieldValue< FieldSpecificationMultiply >( targetSet, 0.0, targetGroup, fieldName ); + } ); + } + } ); +} } /* namespace geos */ diff --git a/src/coreComponents/fieldSpecification/FieldSpecificationManager.hpp b/src/coreComponents/fieldSpecification/FieldSpecificationManager.hpp index be7f8ffec7b..870f8112408 100644 --- a/src/coreComponents/fieldSpecification/FieldSpecificationManager.hpp +++ b/src/coreComponents/fieldSpecification/FieldSpecificationManager.hpp @@ -178,6 +178,12 @@ class FieldSpecificationManager : public dataRepository::Group */ void applyInitialConditions( MeshLevel & mesh ) const; + /** + * @brief function to apply initial conditions which involve scaling a mesh field + * @param mesh the MeshLevel object + */ + void applyScalingInitialConditions( MeshLevel & mesh ) const; + /** * @brief function to validate the application of boundary conditions * @param mesh the MeshLevel object @@ -251,6 +257,7 @@ FieldSpecificationManager:: Group & targetGroup, string const & targetField ) { + fs.applyFieldValue< FieldSpecificationEqual, POLICY >( targetSet, time, targetGroup, targetField ); lambda( fs, targetSet ); } ); diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp index 652c1d3193d..57e34a34356 100644 --- a/src/coreComponents/mainInterface/ProblemManager.cpp +++ b/src/coreComponents/mainInterface/ProblemManager.cpp @@ -1159,6 +1159,7 @@ void ProblemManager::applyInitialConditions() if( !meshLevel.isShallowCopy() ) // to avoid messages printed three times { m_fieldSpecificationManager->validateBoundaryConditions( meshLevel ); + m_fieldSpecificationManager->applyScalingInitialConditions( meshLevel );// with current implementation, scaling only works with shallow copy mesh levels } m_fieldSpecificationManager->applyInitialConditions( meshLevel ); } );