From cbbd0245e60ac8c4c1fd063090243275db7d646b Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 17 Feb 2026 17:13:17 +0100 Subject: [PATCH 1/8] re-implement velocities --- src/coreComponents/common/DataLayouts.hpp | 9 +++ .../finiteVolume/CellElementStencilTPFA.cpp | 64 +++++++++++++++++++ .../finiteVolume/CellElementStencilTPFA.hpp | 12 ++++ .../fluidFlow/CompositionalMultiphaseBase.cpp | 34 ++++++++++ .../fluidFlow/CompositionalMultiphaseBase.hpp | 5 ++ .../CompositionalMultiphaseBaseFields.hpp | 10 ++- .../fluidFlow/CompositionalMultiphaseFVM.cpp | 5 ++ .../compositional/FluxComputeKernel.hpp | 34 ++++++++++ .../compositional/FluxComputeKernelBase.cpp | 1 + .../compositional/FluxComputeKernelBase.hpp | 12 +++- 10 files changed, 183 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/common/DataLayouts.hpp b/src/coreComponents/common/DataLayouts.hpp index cd09af18ef0..fc051ccbaa8 100644 --- a/src/coreComponents/common/DataLayouts.hpp +++ b/src/coreComponents/common/DataLayouts.hpp @@ -272,6 +272,9 @@ using LAYOUT_OBL_OPERATOR_VALUES = RAJA::PERM_JI; /// OBL operator derivatives derivative array layout using LAYOUT_OBL_OPERATOR_DERIVATIVES = RAJA::PERM_JKI; +/// phase velocity used in dispersion layout +using LAYOUT_PHASE_VELOCITY = RAJA::PERM_JKI; + #else /// Component global density/fraction array layout @@ -298,6 +301,9 @@ using LAYOUT_OBL_OPERATOR_VALUES = RAJA::PERM_IJ; /// OBL operator derivatives derivative array layout using LAYOUT_OBL_OPERATOR_DERIVATIVES = RAJA::PERM_IJK; +/// phase velocity used in dispersion layout +using LAYOUT_PHASE_VELOCITY = RAJA::PERM_IJK; + #endif /// Component global density/fraction unit stride dimension @@ -324,6 +330,9 @@ static constexpr int USD_OBL_VAL = LvArray::typeManipulation::getStrideOneDimens /// OBL operator derivatives unit stride dimension static constexpr int USD_OBL_DER = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_OBL_OPERATOR_DERIVATIVES{} ); +/// Phase velocity used in dispersion +static constexpr int USD_PHASE_VELOCITY = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_PHASE_VELOCITY{} ); + } // namespace compflow } // namespace geos diff --git a/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp b/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp index d95969e7a48..a1e9331dfcf 100644 --- a/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp +++ b/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp @@ -120,4 +120,68 @@ CellElementStencilTPFAWrapper:: m_geometricStabilizationCoef( geometricStabilizationCoef ) {} +void +CellElementStencilTPFAWrapper::computeVelocity( + localIndex iconn, localIndex ip, + const real64 (&phaseFlux), + arraySlice1d< real64 const > const (&cellCartDim)[2], + localIndex const (&ghostRank)[2], + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ) +{ + + real64 surface[2]; + + for( localIndex i = 0; i < 2; i++ ) + { + localIndex const er = m_elementRegionIndices[iconn][i]; + localIndex const esr = m_elementSubRegionIndices[iconn][i]; + localIndex const ei = m_elementIndices[iconn][i]; + + if( ghostRank[i] < 0 ) + { + //halfWeight is d(Cell,Face)/area, we want area + real64 const halfWeight = m_weights[iconn][i]; + real64 c2fVec[3]; + LvArray::tensorOps::copy< 3 >( c2fVec, m_cellToFaceVec[iconn][i] ); + real64 const c2fDistance = LvArray::tensorOps::normalize< 3 >( c2fVec ); + + surface[i] = c2fDistance * halfWeight; + + real64 faceNormal[3], invDist[3], velocityNorm[3], phaseVel[3]; + LvArray::tensorOps::copy< 3 >( faceNormal, m_faceNormal[iconn] ); + + LvArray::tensorOps::scale< 3 >( faceNormal, phaseFlux / surface[i] ); + //change sign + if( LvArray::tensorOps::AiBi< 3 >( m_cellToFaceVec[iconn][i], faceNormal ) < 0.0 ) + { + LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); + } + + LvArray::tensorOps::hadamardProduct< 3 >( velocityNorm, faceNormal, m_cellToFaceVec[iconn][i] ); + for( int dir = 0; dir < 3; ++dir ) + { + invDist[dir] = (LvArray::math::abs( cellCartDim[i][dir] ) > LvArray::NumericLimits< real64 >::epsilon) ? + 1. / cellCartDim[i][dir] : LvArray::NumericLimits< real64 >::epsilon; + } + LvArray::tensorOps::hadamardProduct< 3 >( phaseVel, velocityNorm, invDist ); + LvArray::tensorOps::add< 3 >( phaseVelocity[er][esr][ei][ip], phaseVel ); + } + } +} + +void +CellElementStencilTPFAWrapper::initVelocity( localIndex iconn, + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ) +{ + for( localIndex i = 0; i < 2; i++ ) + { + localIndex const er = m_elementRegionIndices[iconn][i]; + localIndex const esr = m_elementSubRegionIndices[iconn][i]; + + phaseVelocity[er][esr].zero(); + + } +} + + } /* namespace geos */ diff --git a/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp b/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp index 88fa3bb2b27..094aec0b110 100644 --- a/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp +++ b/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp @@ -25,6 +25,7 @@ namespace geos { + /** * Provides access to the cellElement stencil that may be called from a kernel function. */ @@ -32,6 +33,17 @@ class CellElementStencilTPFAWrapper : public StencilWrapperBase< TwoPointStencil { public: + void + computeVelocity( + localIndex iconn, localIndex ip, + const real64 (&phaseFlux), + arraySlice1d< real64 const > const (&cellCartDim)[2], + localIndex const (&ghostRank)[2], + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ); + void + initVelocity( localIndex iconn, + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ); + /// Coefficient view accessory type template< typename VIEWTYPE > using CoefficientAccessor = ElementRegionManager::ElementViewConst< VIEWTYPE >; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index a4a9039eb2c..ca810354ff7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -72,6 +72,7 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, m_hasCapPressure( false ), m_hasDiffusion( false ), m_hasDispersion( false ), + m_hasVelocityComputed( 0 ), m_minScalingFactor( 0.01 ), m_allowCompDensChopping( 1 ), m_useTotalMassEquation( 1 ), @@ -167,6 +168,12 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, setApplyDefaultValue( 1 ). setDescription( "Flag indicating whether simple accumulation form is used" ); + this->registerWrapper( viewKeyStruct::hasVelocityComputedString(), &m_hasVelocityComputed ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "Flag indicating whether reconstructed velocity is computed" ); + this->registerWrapper( viewKeyStruct::minCompDensString(), &m_minCompDens ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -300,6 +307,7 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) // If at least one region has a dispersion model, consider it enabled for all m_hasDispersion |= !getConstitutiveName< DispersionBase >( subRegion ).empty(); + m_hasVelocityComputed += integer( m_hasDispersion );//activate velocity if dispersion is on GEOS_ERROR_IF( m_hasDispersion, "Dispersion is not supported yet, please remove it from this XML file" ); } ); @@ -441,6 +449,18 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) setDimLabels( 1, fluid.componentNames() ). reference().resizeDimension< 1 >( m_numComponents ); + if( m_hasVelocityComputed ) + { + + array1d< std::string > directions( 3 ); + directions[0] = "x"; directions[1] = "y"; directions[2] = "z"; + subRegion.registerField< flow::phaseVelocity >( getName()). + setDimLabels( 1, fluid.phaseNames() ). + setDimLabels( 2, directions ). + reference().resizeDimension< 1, 2 >( m_numPhases, directions.size() ); + } + + } ); FaceManager & faceManager = mesh.getFaceManager(); @@ -1053,6 +1073,7 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, arrayView1d< real64 const > const temperature = subRegion.template getField< flow::temperature >(); diffusionMaterial.initializeTemperatureState( temperature ); } + if( m_hasDispersion ) { string const & dispersionName = subRegion.template getReference< string >( viewKeyStruct::dispersionNamesString() ); @@ -1060,8 +1081,15 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, GEOS_UNUSED_VAR( dispersionMaterial ); // TODO: compute the phase velocities here //dispersionMaterial.saveConvergedVelocitySate( phaseVelovity ); + // arrayView3d< real64 const > const phaseVelocity = subRegion.template getField< fields::flow::phaseVelocity >(); + // if( m_useMass ) + // dispersionMaterial.initializeVelocityState( phaseVelocity, fluid.phaseMassDensity()); + // else + // dispersionMaterial.initializeVelocityState( phaseVelocity, fluid.phaseDensity()); + } + } ); } @@ -2882,7 +2910,13 @@ void CompositionalMultiphaseBase::implicitStepComplete( real64 const & time, GEOS_UNUSED_VAR( dispersionMaterial ); // TODO: compute the total velocity here //dispersionMaterial.saveConvergedVelocitySate( totalVelovity ); + // if( m_useMass ) + // this->saveConvergedVelocityState( velocity, fluidMaterial.phaseMassDensity() ); + // else + // this->saveConvergedVelocityState( velocity, fluidMaterial.phaseDensity() ); } + + } ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 227055bf3ed..54f1896a15d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -289,6 +289,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase static constexpr char const * useMassFlagString() { return "useMass"; } static constexpr char const * formulationTypeString() { return "formulationType"; } + static constexpr char const * hasVelocityComputedString() { return "hasVelocityComputed"; } // time stepping controls @@ -301,6 +302,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase static constexpr char const * targetFlowCFLString() { return "targetFlowCFL"; } + // nonlinear solver parameters static constexpr char const * maxCompFracChangeString() { return "maxCompFractionChange"; } @@ -495,6 +497,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase /// flag to determine whether or not to apply dispersion bool m_hasDispersion; + /// flag to determine whether or not to apply assemble velocities + int m_hasVelocityComputed; + /// maximum (absolute) change in a component fraction in a Newton iteration real64 m_maxCompFracChange; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp index 31120bcbd40..4099f68cde9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp @@ -39,7 +39,7 @@ using array3dLayoutPhase_dC = array3d< real64, compflow::LAYOUT_PHASE_DC >; using array2dLayoutComp = array2d< real64, compflow::LAYOUT_COMP >; using array3dLayoutComp_dC = array3d< real64, compflow::LAYOUT_COMP_DC >; using array3dLayoutPhaseComp = array3d< real64, compflow::LAYOUT_PHASE_COMP >; - +using array3dLayoutPhase = array3d< real64, compflow::LAYOUT_PHASE_VELOCITY >; DECLARE_FIELD( globalCompDensity, "globalCompDensity", array2dLayoutComp, @@ -88,7 +88,13 @@ DECLARE_FIELD( globalCompFraction_n, // NOPLOT, // NO_WRITE, // "Global component fraction updates at the previous sequential iteration" ); - +DECLARE_FIELD( phaseVelocity, + "cellCenterPhaseVelocity", + array3dLayoutPhase, + 1, + LEVEL_0, + WRITE_AND_READ, + "Molar/Mass weighted phase velocities reconstructed at cell center" ); DECLARE_FIELD( bcGlobalCompFraction, "bcGlobalCompFraction", array2dLayoutComp, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 40a8b4ddd62..a477b1b57fc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -257,6 +257,9 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, kernelFlags.set( KernelFlags::TotalMassEquation ); if( m_gravityDensityScheme == GravityDensityScheme::PhasePresence ) kernelFlags.set( KernelFlags::CheckPhasePresenceInGravity ); + if( m_hasVelocityComputed ) + kernelFlags.set( KernelFlags::VelocityCompute ); + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, @@ -417,6 +420,8 @@ void CompositionalMultiphaseFVM::assembleStabilizedFluxTerms( real64 const dt, kernelFlags.set( KernelFlags::CapPressure ); if( m_useTotalMassEquation ) kernelFlags.set( KernelFlags::TotalMassEquation ); + if( m_hasVelocityComputed ) + kernelFlags.set( KernelFlags::VelocityCompute ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp index 4d36066b97c..bc030637d61 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp @@ -44,6 +44,8 @@ namespace geos namespace isothermalCompositionalMultiphaseFVMKernels { +class CellElementStencilTPFAWrapper; + /** * @class FluxComputeKernel * @tparam NUM_COMP number of fluid components @@ -174,6 +176,17 @@ class FluxComputeKernel : public FluxComputeKernelBase stackArray2d< real64, maxNumElems * numEqn * maxStencilSize * numDof > localFluxJacobian; }; +/** + * @brief Initialize velocity container + * @param[in] iconn the connection index + */ + GEOS_HOST_DEVICE + inline + void initVelocity( localIndex const iconn ) const + { + if constexpr ( std::is_same< CellElementStencilTPFAWrapper, STENCILWRAPPER >::value) + m_stencilWrapper.initVelocity( m_stencilWrapper, iconn, m_phaseVelocity ); + } /** * @brief Getter for the stencil size at this connection @@ -373,6 +386,23 @@ class FluxComputeKernel : public FluxComputeKernelBase phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans, compFlux, dCompFlux_dP, dCompFlux_dC, dCompFlux_dTrans ); + + if( m_kernelFlags.isSet( KernelFlags::VelocityCompute )) + { + if constexpr (std::is_same< CellElementStencilTPFAWrapper, STENCILWRAPPER >::value) { + m_stencilWrapper.computeVelocity( m_stencilWrapper, + iconn, ip, + phaseFlux, + // {m_globalCellDims[seri[0]][sesri[0]][sei[0]], + // m_globalCellDims[seri[1]][sesri[1]][sei[1]]}, + {m_ghostRank[seri[0]][sesri[0]][sei[0]], + m_ghostRank[seri[1]][sesri[1]][sei[1]]},//TODO replace once decided on globalDims + {m_ghostRank[seri[0]][sesri[0]][sei[0]], + m_ghostRank[seri[1]][sesri[1]][sei[1]]}, + m_phaseVelocity ); + } + } + // 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, m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ), @@ -477,6 +507,10 @@ class FluxComputeKernel : public FluxComputeKernelBase KERNEL_TYPE const & kernelComponent ) { GEOS_MARK_FUNCTION; + //velocity has to be reset on all faces prior computing flux. (to be specific on all stencil extends) + forAll< POLICY >( numConnections, [=] GEOS_HOST_DEVICE ( localIndex const iconn ) { + kernelComponent.initVelocity( iconn ); + } ); forAll< POLICY >( numConnections, [=] GEOS_HOST_DEVICE ( localIndex const iconn ) { typename KERNEL_TYPE::StackVariables stack( kernelComponent.stencilSize( iconn ), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.cpp index 9034e85b645..77b2b75263b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.cpp @@ -49,6 +49,7 @@ FluxComputeKernelBase::FluxComputeKernelBase( integer const numPhases, m_phaseVolFrac( compFlowAccessors.get( flow::phaseVolumeFraction {} ) ), m_dPhaseVolFrac( compFlowAccessors.get( flow::dPhaseVolumeFraction {} ) ), m_dCompFrac_dCompDens( compFlowAccessors.get( flow::dGlobalCompFraction_dGlobalCompDensity {} ) ), + m_phaseVelocity( compFlowAccessors.get( flow::phaseVelocity {} ) ), m_phaseCompFrac( multiFluidAccessors.get( multifluid::phaseCompFraction {} ) ), m_dPhaseCompFrac( multiFluidAccessors.get( multifluid::dPhaseCompFraction {} ) ), m_localMatrix( localMatrix ), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp index bc46059cec1..581252aac33 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp @@ -58,7 +58,9 @@ enum class KernelFlags /// Flag indicating whether IHU is used or not IHU = 1 << 6, // 64 /// Flag indicating whether HU 2-phase simplified version is used or not - HU2PH = 1 << 7 // 128 + HU2PH = 1 << 7, // 128 + /// Flag indicting we must compute cell-centered velocity + VelocityCompute = 1 << 8 // 256 }; /******************************** FluxComputeKernelBase ********************************/ @@ -80,6 +82,9 @@ class FluxComputeKernelBase template< typename VIEWTYPE > using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + template< typename VIEWTYPE > + using ElementView = ElementRegionManager::ElementView< VIEWTYPE >; + using DofNumberAccessor = ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > >; using CompFlowAccessors = @@ -89,6 +94,7 @@ class FluxComputeKernelBase fields::flow::dGlobalCompFraction_dGlobalCompDensity, fields::flow::phaseVolumeFraction, fields::flow::dPhaseVolumeFraction, + fields::flow::phaseVelocity, fields::flow::phaseMobility, fields::flow::dPhaseMobility >; using MultiFluidAccessors = @@ -162,10 +168,14 @@ class FluxComputeKernelBase /// Views on derivatives of comp fractions ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dCompFrac_dCompDens; + /// Views on phase velocity + ElementView< arrayView3d< real64, compflow::USD_PHASE_VELOCITY > > const m_phaseVelocity; + /// Views on phase component fractions ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const m_phaseCompFrac; ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const m_dPhaseCompFrac; + // Residual and jacobian /// View on the local CRS matrix From 3aadacabf8211fa69dd020e7cf945a2e5fb0ab41 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 17 Feb 2026 18:10:55 +0100 Subject: [PATCH 2/8] adding test --- .../finiteVolume/CellElementStencilTPFA.cpp | 14 +- .../finiteVolume/CellElementStencilTPFA.hpp | 14 +- .../compositional/FluxComputeKernel.hpp | 18 +-- .../fluidFlow/unitTests/CMakeLists.txt | 1 + .../unitTests/testVelocityReconstruction.cpp | 151 ++++++++++++++++++ 5 files changed, 175 insertions(+), 23 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp diff --git a/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp b/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp index a1e9331dfcf..c37bed37c14 100644 --- a/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp +++ b/src/coreComponents/finiteVolume/CellElementStencilTPFA.cpp @@ -122,11 +122,11 @@ CellElementStencilTPFAWrapper:: void CellElementStencilTPFAWrapper::computeVelocity( - localIndex iconn, localIndex ip, - const real64 (&phaseFlux), - arraySlice1d< real64 const > const (&cellCartDim)[2], - localIndex const (&ghostRank)[2], - ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ) + localIndex iconn, localIndex ip, + const real64 (&phaseFlux), + arraySlice1d< real64 const > const (&cellCartDim)[2], + localIndex const (&ghostRank)[2], + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ) { real64 surface[2]; @@ -170,8 +170,8 @@ CellElementStencilTPFAWrapper::computeVelocity( } void -CellElementStencilTPFAWrapper::initVelocity( localIndex iconn, - ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ) +CellElementStencilTPFAWrapper::initVelocity( localIndex iconn, + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ) { for( localIndex i = 0; i < 2; i++ ) { diff --git a/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp b/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp index 094aec0b110..01da35e8b65 100644 --- a/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp +++ b/src/coreComponents/finiteVolume/CellElementStencilTPFA.hpp @@ -34,15 +34,15 @@ class CellElementStencilTPFAWrapper : public StencilWrapperBase< TwoPointStencil public: void - computeVelocity( - localIndex iconn, localIndex ip, - const real64 (&phaseFlux), - arraySlice1d< real64 const > const (&cellCartDim)[2], - localIndex const (&ghostRank)[2], - ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ); + computeVelocity( + localIndex iconn, localIndex ip, + const real64 (&phaseFlux), + arraySlice1d< real64 const > const (&cellCartDim)[2], + localIndex const (&ghostRank)[2], + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ); void initVelocity( localIndex iconn, - ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ); + ElementRegionManager::ElementView< arrayView3d< real64 > > const & phaseVelocity ); /// Coefficient view accessory type template< typename VIEWTYPE > diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp index bc030637d61..523f40b4494 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp @@ -391,15 +391,15 @@ class FluxComputeKernel : public FluxComputeKernelBase { if constexpr (std::is_same< CellElementStencilTPFAWrapper, STENCILWRAPPER >::value) { m_stencilWrapper.computeVelocity( m_stencilWrapper, - iconn, ip, - phaseFlux, - // {m_globalCellDims[seri[0]][sesri[0]][sei[0]], - // m_globalCellDims[seri[1]][sesri[1]][sei[1]]}, - {m_ghostRank[seri[0]][sesri[0]][sei[0]], - m_ghostRank[seri[1]][sesri[1]][sei[1]]},//TODO replace once decided on globalDims - {m_ghostRank[seri[0]][sesri[0]][sei[0]], - m_ghostRank[seri[1]][sesri[1]][sei[1]]}, - m_phaseVelocity ); + iconn, ip, + phaseFlux, + // {m_globalCellDims[seri[0]][sesri[0]][sei[0]], + // m_globalCellDims[seri[1]][sesri[1]][sei[1]]}, + {m_ghostRank[seri[0]][sesri[0]][sei[0]], + m_ghostRank[seri[1]][sesri[1]][sei[1]]},//TODO replace once decided on globalDims + {m_ghostRank[seri[0]][sesri[0]][sei[0]], + m_ghostRank[seri[1]][sesri[1]][sei[1]]}, + m_phaseVelocity ); } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt index 77d628e7ae2..cdc5bae8516 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt @@ -4,6 +4,7 @@ Builds unit test for fluid flow solver statistics #]] set( fluidFlow_gtest_tests + testVelocityReconstruction.cpp testFlowStatistics.cpp ) set( dependencyList mainInterface testingUtilities gtest ${parallelDeps} ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp new file mode 100644 index 00000000000..efb05e8d265 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp @@ -0,0 +1,151 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + + +// Source includes +#include "mainInterface/initialization.hpp" +#include "finiteVolume/CellElementStencilTPFA.hpp" +// #include "testFlowKernelHelpers.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/unitTests/testFlowKernelHelpers.hpp" + +// TPL includes +#include + + +using namespace geos; + + +void setReferences( const real64 flux, const real64 dx, const real64 dy, real64 (& expectedVel)[3] ) +{ + + + expectedVel[0] = flux/dx; + expectedVel[1] = flux/dy; + expectedVel[2] = 0.0; + + +} + + +TEST( testAligned2D, Velocity_aligned2D ) +{ + constexpr int nfaces = 4; //2-hexa problem + constexpr int numElemtInStencil = 2; + + + real64 const faceNormal[nfaces][3] = {{1., 0., 0.}, + {1., 0., 0.}, + {0., 1., 0.}, + {0., 1., 0.}}; + real64 const cellToFaceVec[nfaces][numElemtInStencil][3] = { {{1., 0., 0.}, {-1., 0., 0}}, + {{1., 0., 0.}, {-1., 0., 0}}, + {{0., 1., 0.}, {0., -1., 0}}, + {{0., 1., 0.}, {0., -1., 0}} + }; + real64 const transMult[nfaces] = {1., 1., 1., 1.}; + real64 weight[] = {1, -1}; + real64 const geomStabSum = 1.; //irrelevant for now + + + localIndex elementRegionIndices[] = {0, 0}; + localIndex elementSubRegionIndices[] = {0, 0}; + localIndex ei_[4][numElemtInStencil] = { + {0, 1}, //face #1 + {2, 3}, //face #2 + {0, 2}, //face #3 + {1, 3} //face #4 + }; + + + const real64 flux = 0.001; + const int numElem = 4; + constexpr localIndex nPhases = 2; + constexpr localIndex nDirs = 3; + const real64 phaseVelocity[numElem][nPhases][nDirs] = { + {{0., 0., 0.}, {0., 0., 0.}}, + {{0., 0., 0.}, {0., 0., 0.}}, + {{0., 0., 0.}, {0., 0., 0.}}, + {{0., 0., 0.}, {0., 0., 0.}} + }; + + + CellElementStencilTPFA tpfa; + for( int kf = 0; kf < nfaces; ++kf ) + { + localIndex elementIndices[] = {ei_[kf][0], ei_[kf][1]}; + tpfa.add( 2, elementRegionIndices, elementSubRegionIndices, elementIndices, weight, kf ); + tpfa.addVectors( transMult[kf], geomStabSum, faceNormal[kf], cellToFaceVec[kf] ); + } + + + CellElementStencilTPFA::KernelWrapper wrapper = tpfa.createKernelWrapper(); + array2d< real64 > globalCellDim( 2, 3 ); + globalCellDim[0][0] = 2.; globalCellDim[1][0] = 2.; + globalCellDim[0][1] = 2.; globalCellDim[1][1] = 2.; + globalCellDim[0][2] = 0.; globalCellDim[1][2] = 0.; + arrayView2d< const real64 > const globalCellDimView = globalCellDim.toViewConst(); + + + CellElementStencilTPFA::IndexContainerViewConstType const & seri = tpfa.getElementRegionIndices(); + CellElementStencilTPFA::IndexContainerViewConstType const & sesri = tpfa.getElementSubRegionIndices(); + CellElementStencilTPFA::IndexContainerViewConstType const & sei = tpfa.getElementIndices(); + + + ElementRegionManager::ElementViewAccessor< array3d< real64 > > phaseVelocityView = AccessorHelper< true >::makeElementAccessor< 3 >( phaseVelocity[0][0], + tpfa.stencilSize( 0 ), + seri[nfaces-1], + sesri[nfaces-1], + sei[nfaces-1], + nPhases, + nDirs ); + for( int iconn = 0; iconn < nfaces; ++iconn ) + { + + + tpfa.computeVelocity( wrapper, iconn /*iconn*/, 0 /*ip*/, flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1 }, phaseVelocityView.toNestedView()); + tpfa.computeVelocity( wrapper, iconn /*iconn*/, 1 /*ip*/, 100*flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1}, phaseVelocityView.toNestedView()); + } + + + for( int ip = 0; ip < 2; ++ip ) + { + auto fip = (99*ip+1)*flux; //shorthand to get the 1:100 ratio between phases as imposed above + for( int ib = 0; ib < numElem; ++ib ) + { + real64 expectedVel[3]; + setReferences( fip, globalCellDim[0][0], globalCellDim[0][1], expectedVel ); + EXPECT_EQ( phaseVelocityView[0][0][ib][ip][0], expectedVel[0] ); + EXPECT_EQ( phaseVelocityView[0][0][ib][ip][1], expectedVel[1] ); + EXPECT_EQ( phaseVelocityView[0][0][ib][ip][2], expectedVel[2] ); + } + } +} + + +int main( int argc, char * *argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + + geos::basicSetup( argc, argv ); + + + int const result = RUN_ALL_TESTS(); + + + geos::basicCleanup(); + + + return result; +} \ No newline at end of file From 02002192c380bc47c19d6fcd4ff906ebac5d1b90 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 17 Feb 2026 20:39:07 +0100 Subject: [PATCH 3/8] restore test working --- .../fluidFlow/unitTests/testVelocityReconstruction.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp index efb05e8d265..242429c34b1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp @@ -113,8 +113,8 @@ TEST( testAligned2D, Velocity_aligned2D ) { - tpfa.computeVelocity( wrapper, iconn /*iconn*/, 0 /*ip*/, flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1 }, phaseVelocityView.toNestedView()); - tpfa.computeVelocity( wrapper, iconn /*iconn*/, 1 /*ip*/, 100*flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1}, phaseVelocityView.toNestedView()); + wrapper.computeVelocity( iconn /*iconn*/, 0 /*ip*/, flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1 }, phaseVelocityView.toNestedView()); + wrapper.computeVelocity( iconn /*iconn*/, 1 /*ip*/, 100*flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1}, phaseVelocityView.toNestedView()); } From aa2bc6f10e376274da339fe59b35b24b304b320d Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 17 Feb 2026 20:39:35 +0100 Subject: [PATCH 4/8] starting operator machinery --- .../mesh/utilities/ComputationalGeometry.hpp | 29 +++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 61de4181531..7874291afc9 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -284,6 +284,35 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, return area; } +//Compute LSQ mat per cell for lsq interpolation +template< typename MATRIX_TYPE, typename CELL_TYPE > +GEOS_HOST_DEVICE +auto PrecomputeLSQInterpolant( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, + MATRIX_TYPE&& lsq_op ) +{ + + MATRIX_TYPE ntn, n, lsq_op; + n[ 0 ][ 0 ] = normals[ 0 ][ 0 ]; + n[ 0 ][ 1 ] = normals[ 0 ][ 1 ]; + n[ 0 ][ 2 ] = normals[ 0 ][ 2 ]; + + n[ 1 ][ 0 ] = normals[ 1 ][ 0 ]; + n[ 1 ][ 1 ] = normals[ 1 ][ 1 ]; + n[ 1 ][ 2 ] = normals[ 1 ][ 2 ]; + + n[ 2 ][ 0 ] = normals[ 2 ][ 0 ]; + n[ 2 ][ 1 ] = normals[ 2 ][ 1 ]; + n[ 2 ][ 2 ] = normals[ 2 ][ 2 ]; + // N is nfaces x 3 matrices of normal + // (N^T N)^(-1)*N^T is the interpolant + LvArray::tensorOps::Rij_add_AikAjk<3, CELL_TYPE::NFACE >(ntn, n); + LvArray::tensorOps::symInvert(lsq_op_0, ntn); + LvArray::tensorOps::Rij_add_AikBjk(lsq_op,lsq_op_0,n); + + return lsq_op; +} + + /** * @brief Change the orientation of the input vector to be consistent in a global sense. * @tparam NORMAL_TYPE type of @p normal From 8cebb504254be12742126883fa292491920eac1d Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 18 Feb 2026 14:45:16 +0100 Subject: [PATCH 5/8] with tensorOps --- .../mesh/utilities/ComputationalGeometry.hpp | 75 ++++++++++++++----- 1 file changed, 56 insertions(+), 19 deletions(-) diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 7874291afc9..a3dcad62b8b 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -284,32 +284,69 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, return area; } -//Compute LSQ mat per cell for lsq interpolation -template< typename MATRIX_TYPE, typename CELL_TYPE > +// //Compute LSQ mat per cell for lsq interpolation +// template< typename CELL > +// array2d computeVelocity(arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes) +// { + +// if(normals.size(0) == 2){ +// // quad -- split cell Dim ? + +// // tri + +// } +// else if (normals.size(1) == 3) { +// // hexahedron + +// // tetra + +// // wedge +// // pyramid +// } + + + +// } + +template< int NFACES > GEOS_HOST_DEVICE -auto PrecomputeLSQInterpolant( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, - MATRIX_TYPE&& lsq_op ) +array2d computeVelocities( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes) { + array2d velocity; - MATRIX_TYPE ntn, n, lsq_op; - n[ 0 ][ 0 ] = normals[ 0 ][ 0 ]; - n[ 0 ][ 1 ] = normals[ 0 ][ 1 ]; - n[ 0 ][ 2 ] = normals[ 0 ][ 2 ]; + //check compatibility on sizes + GEOS_ERROR_IF( (normals.size(1)!=fluxes.size()), GEOS_FMT("Error in parameters: normal({}) and fluxes({}) are different in sizes.", normals.size(1), fluxes.size()) ); - n[ 1 ][ 0 ] = normals[ 1 ][ 0 ]; - n[ 1 ][ 1 ] = normals[ 1 ][ 1 ]; - n[ 1 ][ 2 ] = normals[ 1 ][ 2 ]; + //TODO dispatch as a function of dims and nfaces + real64 velocity_c[3][1];//, fluxes_c[NFACES][1]; + real64 n[ NFACES ][ 3 ], ntn[ 3 ][ 3 ], ope[3][NFACES]; + + // + n[0][0] = normals[0][0]; + n[0][1] = normals[0][1]; + n[0][2] = normals[0][2]; + + n[1][0] = normals[1][0]; + n[1][1] = normals[1][1]; + n[1][2] = normals[1][2]; + + n[2][0] = normals[2][0]; + n[2][1] = normals[2][1]; + n[2][2] = normals[2][2]; - n[ 2 ][ 0 ] = normals[ 2 ][ 0 ]; - n[ 2 ][ 1 ] = normals[ 2 ][ 1 ]; - n[ 2 ][ 2 ] = normals[ 2 ][ 2 ]; // N is nfaces x 3 matrices of normal // (N^T N)^(-1)*N^T is the interpolant - LvArray::tensorOps::Rij_add_AikAjk<3, CELL_TYPE::NFACE >(ntn, n); - LvArray::tensorOps::symInvert(lsq_op_0, ntn); - LvArray::tensorOps::Rij_add_AikBjk(lsq_op,lsq_op_0,n); - - return lsq_op; + LvArray::tensorOps::Rij_add_AikAjk<3,NFACES>(ntn, n); + assert(LvArray::tensorOps::determinant<3>(ntn) > 1e-12); // else 2D most likely + //to SymMatrix + LvArray::tensorOps::symInvert<3>(ntn); + //from Symmatrix + LvArray::tensorOps::Rij_add_AikBjk<3,NFACES>(ope, ntn,n); + LvArray::tensorOps::Ri_add_AijBj<3,NFACES>(velocity, ope, fluxes); + + LvArray::tensorOps::copy(velocity,velocity_c); + + return velocity; } From 252a2eaf756fd515232e796f92b288a17df22c91 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 18 Feb 2026 16:09:49 +0100 Subject: [PATCH 6/8] with analytics --- .../mesh/utilities/ComputationalGeometry.cpp | 52 ++++++- .../mesh/utilities/ComputationalGeometry.hpp | 140 ++++++++++++------ 2 files changed, 149 insertions(+), 43 deletions(-) diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp index c920605de80..bae67cbdec6 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp @@ -23,6 +23,56 @@ namespace geos { namespace computationalGeometry -{} /* namespace computationalGeometry */ +{ + +array2d computeVelocity(arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, ElementType& elem /*subRegion.getElementType()*/) +{ + +switch (elem) { + + case geos::ElementType::Triangle: + //TODO pre-post rotate + GEOS_ERROR("Not Implemented Yet !"); + break; + // return computeVelocities_<3>(normals,fluxes, rotation); + case geos::ElementType::Quadrilateral: + //TODO pre-post rotate + GEOS_ERROR("Not Implemented Yet !"); + break; + // return computeVelocities_<4>(normals,fluxes, rotation); + case geos::ElementType::Tetrahedron: + return computeVelocities_<4>(normals,fluxes); + case geos::ElementType::Pyramid: + case geos::ElementType::Prism5: + case geos::ElementType::Wedge: + return computeVelocities_<5>(normals,fluxes); + case geos::ElementType::Hexahedron: + case geos::ElementType::Prism6: + return computeVelocities_<6>(normals,fluxes); + case geos::ElementType::Prism7: + return computeVelocities_<7>(normals,fluxes); + case geos::ElementType::Prism8: + return computeVelocities_<8>(normals,fluxes); + case geos::ElementType::Prism9: + return computeVelocities_<9>(normals,fluxes); + case geos::ElementType::Prism10: + return computeVelocities_<10>(normals,fluxes); + case geos::ElementType::Prism11: + return computeVelocities_<11>(normals,fluxes); + + + case geos::ElementType::Polygon: + case geos::ElementType::Polyhedron: + default: + GEOS_ERROR("Velocities are not computed on 2D Polygons cell type"); + break; + +} + + return {}; +} + + +} /* namespace computationalGeometry */ } /* namespace geos */ diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index a3dcad62b8b..10fcc00de2f 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -284,71 +284,127 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, return area; } -// //Compute LSQ mat per cell for lsq interpolation -// template< typename CELL > -// array2d computeVelocity(arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes) -// { -// if(normals.size(0) == 2){ -// // quad -- split cell Dim ? -// // tri +template< int NFACES> +GEOS_HOST_DEVICE +array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, real64 const (&rotation)[3][3]) +{ + array2d velocity;//should be rotated so that the off plane is aligned with z ? + + //check compatibility on sizes + GEOS_ERROR_IF( (normals.size(1)!=fluxes.size()), GEOS_FMT("Error in parameters: normal({}) and fluxes({}) are different in sizes.", normals.size(1), fluxes.size()) ); + GEOS_ERROR_IF( (NFACES!=fluxes.size()), GEOS_FMT("Error in parameters: templatized size({}) and fluxes({}) are different in sizes.", NFACES, fluxes.size()) ); + + //TODO dispatch as a function of dims and nfaces + real64 velocity_c[2], fluxes_c[NFACES]; + real64 n[ NFACES ][ 2 ], ntn[ 2 ][ 2 ], ope[2][ NFACES ]; -// } -// else if (normals.size(1) == 3) { -// // hexahedron + // N is nfaces x 3 matrices of normal + // (N^T N)^(-1)*N^T is the interpolant + real64 a = 0., b = 0., c = 0.; + for(int i=0; i 1e-12); // else 2D most likely + // analytical inverse for symmetrical matrix + // [a c ] + // [. b ] + ntn[0][0] = 1/detA * (b); + ntn[1][1] = 1/detA * (a); + ntn[0][1] = 1/detA * (-c); + ntn[1][0] = 1/detA * (-c); + + //from Symmatrix + LvArray::tensorOps::Rij_add_AikBjk<2,NFACES>(ope, ntn, n); + LvArray::tensorOps::Ri_add_AijBj<2,NFACES>(velocity_c, ope, fluxes_c); + velocity[0][0] = velocity_c[0]; + velocity[0][1] = velocity_c[1]; + velocity[0][2] = 0.0;//off plane + + return velocity; +} -// } -template< int NFACES > +template< int NFACES> GEOS_HOST_DEVICE -array2d computeVelocities( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes) +array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes) { array2d velocity; //check compatibility on sizes GEOS_ERROR_IF( (normals.size(1)!=fluxes.size()), GEOS_FMT("Error in parameters: normal({}) and fluxes({}) are different in sizes.", normals.size(1), fluxes.size()) ); + GEOS_ERROR_IF( (NFACES!=fluxes.size()), GEOS_FMT("Error in parameters: templatized size({}) and fluxes({}) are different in sizes.", NFACES, fluxes.size()) ); //TODO dispatch as a function of dims and nfaces - real64 velocity_c[3][1];//, fluxes_c[NFACES][1]; - real64 n[ NFACES ][ 3 ], ntn[ 3 ][ 3 ], ope[3][NFACES]; + real64 velocity_c[3], fluxes_c[NFACES]; + real64 n[ NFACES ][ 3 ], ntn[ 3 ][ 3 ], ope[3][ NFACES ]; - // - n[0][0] = normals[0][0]; - n[0][1] = normals[0][1]; - n[0][2] = normals[0][2]; - - n[1][0] = normals[1][0]; - n[1][1] = normals[1][1]; - n[1][2] = normals[1][2]; - - n[2][0] = normals[2][0]; - n[2][1] = normals[2][1]; - n[2][2] = normals[2][2]; - - // N is nfaces x 3 matrices of normal - // (N^T N)^(-1)*N^T is the interpolant - LvArray::tensorOps::Rij_add_AikAjk<3,NFACES>(ntn, n); - assert(LvArray::tensorOps::determinant<3>(ntn) > 1e-12); // else 2D most likely - //to SymMatrix - LvArray::tensorOps::symInvert<3>(ntn); - //from Symmatrix - LvArray::tensorOps::Rij_add_AikBjk<3,NFACES>(ope, ntn,n); - LvArray::tensorOps::Ri_add_AijBj<3,NFACES>(velocity, ope, fluxes); + // N is nfaces x 3 matrices of normal + // (N^T N)^(-1)*N^T is the interpolant + real64 a = 0., b = 0., c = 0., d = 0., e = 0., f = 0.; + for(int i=0; i 1e-12); // else 2D most likely + // analytical inverse for symmetrical matrix + // [a d f] + // [. b e] + // [. . c] + + ntn[0][0] = 1/detA * (b*c-e*e); + ntn[1][1] = 1/detA * (a*c-f*f); + ntn[2][2] = 1/detA * (a*b-d*d); + + ntn[0][1] = 1/detA * (d*c-e*f); + ntn[0][2] = 1/detA * (d*e-f*b); + ntn[1][2] = 1/detA * (a*e-f*d); + + ntn[1][0] = 1/detA * (d*c-e*f); + ntn[2][0] = 1/detA * (d*e-f*b); + ntn[2][1] = 1/detA * (a*e-f*d); + + //from Symmatrix + LvArray::tensorOps::Rij_add_AikBjk<3,NFACES,3>(ope, ntn, n); + LvArray::tensorOps::Ri_add_AijBj<3,NFACES>(velocity_c, ope, fluxes_c); + + velocity[0][0] = velocity_c[0]; + velocity[0][1] = velocity_c[1]; + velocity[0][2] = velocity_c[2]; + return velocity; } +// //Compute LSQ mat per cell for lsq interpolation +array2d computeVelocity(arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, ElementType& elem /*subRegion.getElementType()*/); + /** * @brief Change the orientation of the input vector to be consistent in a global sense. From 6b7f24645a67ec32eb94e7c60ca3bfbb72c8393b Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 18 Feb 2026 20:28:37 +0100 Subject: [PATCH 7/8] wip --- .../mesh/utilities/ComputationalGeometry.hpp | 2 +- .../unitTests/testVelocityReconstruction.cpp | 165 +++++++++--------- 2 files changed, 88 insertions(+), 79 deletions(-) diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 10fcc00de2f..709618c75e0 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -402,7 +402,7 @@ array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_ return velocity; } -// //Compute LSQ mat per cell for lsq interpolation +// Compute Velocity interpolations array2d computeVelocity(arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, ElementType& elem /*subRegion.getElementType()*/); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp index 242429c34b1..f3a948525a9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp @@ -15,9 +15,9 @@ // Source includes #include "mainInterface/initialization.hpp" -#include "finiteVolume/CellElementStencilTPFA.hpp" -// #include "testFlowKernelHelpers.hpp" +// #include "finiteVolume/CellElementStencilTPFA.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/unitTests/testFlowKernelHelpers.hpp" +#include "mesh/utilities/ComputationalGeometry.hpp" // TPL includes #include @@ -38,98 +38,107 @@ void setReferences( const real64 flux, const real64 dx, const real64 dy, real64 } -TEST( testAligned2D, Velocity_aligned2D ) +TEST( testAligned2D, Velocity_Hexa ) { - constexpr int nfaces = 4; //2-hexa problem - constexpr int numElemtInStencil = 2; + constexpr int nfaces = 6; //2-hexa problem - - real64 const faceNormal[nfaces][3] = {{1., 0., 0.}, - {1., 0., 0.}, + real64 const normals[nfaces][3] = {{1., 0., 0.}, + {-1., 0., 0.}, {0., 1., 0.}, - {0., 1., 0.}}; - real64 const cellToFaceVec[nfaces][numElemtInStencil][3] = { {{1., 0., 0.}, {-1., 0., 0}}, - {{1., 0., 0.}, {-1., 0., 0}}, - {{0., 1., 0.}, {0., -1., 0}}, - {{0., 1., 0.}, {0., -1., 0}} - }; - real64 const transMult[nfaces] = {1., 1., 1., 1.}; - real64 weight[] = {1, -1}; - real64 const geomStabSum = 1.; //irrelevant for now - - - localIndex elementRegionIndices[] = {0, 0}; - localIndex elementSubRegionIndices[] = {0, 0}; - localIndex ei_[4][numElemtInStencil] = { - {0, 1}, //face #1 - {2, 3}, //face #2 - {0, 2}, //face #3 - {1, 3} //face #4 + {0., -1., 0.}, + {0., 0., 1.}, + {0., 0., -1.} }; + real64 const fluxes[nfaces] = {1., 1., 1., 1., 1., 1.}; + const real64 flux = 0.001; - const int numElem = 4; + const int numElem = 1; constexpr localIndex nPhases = 2; constexpr localIndex nDirs = 3; - const real64 phaseVelocity[numElem][nPhases][nDirs] = { - {{0., 0., 0.}, {0., 0., 0.}}, - {{0., 0., 0.}, {0., 0., 0.}}, - {{0., 0., 0.}, {0., 0., 0.}}, - {{0., 0., 0.}, {0., 0., 0.}} - }; - - - CellElementStencilTPFA tpfa; - for( int kf = 0; kf < nfaces; ++kf ) - { - localIndex elementIndices[] = {ei_[kf][0], ei_[kf][1]}; - tpfa.add( 2, elementRegionIndices, elementSubRegionIndices, elementIndices, weight, kf ); - tpfa.addVectors( transMult[kf], geomStabSum, faceNormal[kf], cellToFaceVec[kf] ); - } - - - CellElementStencilTPFA::KernelWrapper wrapper = tpfa.createKernelWrapper(); - array2d< real64 > globalCellDim( 2, 3 ); - globalCellDim[0][0] = 2.; globalCellDim[1][0] = 2.; - globalCellDim[0][1] = 2.; globalCellDim[1][1] = 2.; - globalCellDim[0][2] = 0.; globalCellDim[1][2] = 0.; - arrayView2d< const real64 > const globalCellDimView = globalCellDim.toViewConst(); - - - CellElementStencilTPFA::IndexContainerViewConstType const & seri = tpfa.getElementRegionIndices(); - CellElementStencilTPFA::IndexContainerViewConstType const & sesri = tpfa.getElementSubRegionIndices(); - CellElementStencilTPFA::IndexContainerViewConstType const & sei = tpfa.getElementIndices(); + // const real64 phaseVelocity[numElem][nPhases][nDirs] = { + // {{0., 0., 0.}, {0., 0., 0.}}, + // }; + ElementRegionManager::ElementViewAccessor< array2d< real64 > > normalsView = AccessorHelper< true >::makeElementAccessor< 2 >( normals[0][0], + 2, + seri[nfaces-1], + sesri[nfaces-1], + sei[nfaces-1], + nPhases, + nDirs ); - ElementRegionManager::ElementViewAccessor< array3d< real64 > > phaseVelocityView = AccessorHelper< true >::makeElementAccessor< 3 >( phaseVelocity[0][0], - tpfa.stencilSize( 0 ), + ElementRegionManager::ElementViewAccessor< array1d< real64 > > fluxesView = AccessorHelper< true >::makeElementAccessor< 1 >( fluxes[0], + 2, seri[nfaces-1], sesri[nfaces-1], sei[nfaces-1], nPhases, nDirs ); - for( int iconn = 0; iconn < nfaces; ++iconn ) - { - - - wrapper.computeVelocity( iconn /*iconn*/, 0 /*ip*/, flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1 }, phaseVelocityView.toNestedView()); - wrapper.computeVelocity( iconn /*iconn*/, 1 /*ip*/, 100*flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1}, phaseVelocityView.toNestedView()); - } - - - for( int ip = 0; ip < 2; ++ip ) - { - auto fip = (99*ip+1)*flux; //shorthand to get the 1:100 ratio between phases as imposed above - for( int ib = 0; ib < numElem; ++ib ) - { - real64 expectedVel[3]; - setReferences( fip, globalCellDim[0][0], globalCellDim[0][1], expectedVel ); - EXPECT_EQ( phaseVelocityView[0][0][ib][ip][0], expectedVel[0] ); - EXPECT_EQ( phaseVelocityView[0][0][ib][ip][1], expectedVel[1] ); - EXPECT_EQ( phaseVelocityView[0][0][ib][ip][2], expectedVel[2] ); - } - } + + // ElementRegionManager::ElementViewAccessor< array3d< real64 > > phaseVelocityView = AccessorHelper< true >::makeElementAccessor< 3 >( phaseVelocity[0][0], + // tpfa.stencilSize( 0 ), + // seri[nfaces-1], + // sesri[nfaces-1], + // sei[nfaces-1], + // nPhases, + // nDirs ); + auto phaseVelocity = computationalGeometry::computeVelocities_(normalsView,fluxesView); + + // CellElementStencilTPFA tpfa; + // for( int kf = 0; kf < nfaces; ++kf ) + // { + // localIndex elementIndices[] = {ei_[kf][0], ei_[kf][1]}; + // tpfa.add( 2, elementRegionIndices, elementSubRegionIndices, elementIndices, weight, kf ); + // tpfa.addVectors( transMult[kf], geomStabSum, faceNormal[kf], cellToFaceVec[kf] ); + // } + + + // CellElementStencilTPFA::KernelWrapper wrapper = tpfa.createKernelWrapper(); + // array2d< real64 > globalCellDim( 2, 3 ); + // globalCellDim[0][0] = 2.; globalCellDim[1][0] = 2.; + // globalCellDim[0][1] = 2.; globalCellDim[1][1] = 2.; + // globalCellDim[0][2] = 0.; globalCellDim[1][2] = 0.; + // arrayView2d< const real64 > const globalCellDimView = globalCellDim.toViewConst(); + + + // CellElementStencilTPFA::IndexContainerViewConstType const & seri = tpfa.getElementRegionIndices(); + // CellElementStencilTPFA::IndexContainerViewConstType const & sesri = tpfa.getElementSubRegionIndices(); + // CellElementStencilTPFA::IndexContainerViewConstType const & sei = tpfa.getElementIndices(); + + + // ElementRegionManager::ElementViewAccessor< array3d< real64 > > phaseVelocityView = AccessorHelper< true >::makeElementAccessor< 3 >( phaseVelocity[0][0], + // tpfa.stencilSize( 0 ), + // seri[nfaces-1], + // sesri[nfaces-1], + // sei[nfaces-1], + // nPhases, + // nDirs ); + + + + // for( int iconn = 0; iconn < nfaces; ++iconn ) + // { + + + // wrapper.computeVelocity( iconn /*iconn*/, 0 /*ip*/, flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1 }, phaseVelocityView.toNestedView()); + // wrapper.computeVelocity( iconn /*iconn*/, 1 /*ip*/, 100*flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1}, phaseVelocityView.toNestedView()); + // } + + + // for( int ip = 0; ip < 2; ++ip ) + // { + // auto fip = (99*ip+1)*flux; //shorthand to get the 1:100 ratio between phases as imposed above + // for( int ib = 0; ib < numElem; ++ib ) + // { + // real64 expectedVel[3]; + // setReferences( fip, globalCellDim[0][0], globalCellDim[0][1], expectedVel ); + // EXPECT_EQ( phaseVelocityView[0][0][ib][ip][0], expectedVel[0] ); + // EXPECT_EQ( phaseVelocityView[0][0][ib][ip][1], expectedVel[1] ); + // EXPECT_EQ( phaseVelocityView[0][0][ib][ip][2], expectedVel[2] ); + // } + // } } From a9a3ae2219173a7ef2fd24d4e9322115cd6fa456 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 19 Feb 2026 17:57:31 +0100 Subject: [PATCH 8/8] debug test --- .../mesh/utilities/ComputationalGeometry.cpp | 49 ---- .../mesh/utilities/ComputationalGeometry.hpp | 179 ++++++++++----- .../unitTests/testVelocityReconstruction.cpp | 216 +++++++++--------- 3 files changed, 228 insertions(+), 216 deletions(-) diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp index bae67cbdec6..d4a7d386f2e 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp @@ -24,55 +24,6 @@ namespace geos namespace computationalGeometry { - -array2d computeVelocity(arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, ElementType& elem /*subRegion.getElementType()*/) -{ - -switch (elem) { - - case geos::ElementType::Triangle: - //TODO pre-post rotate - GEOS_ERROR("Not Implemented Yet !"); - break; - // return computeVelocities_<3>(normals,fluxes, rotation); - case geos::ElementType::Quadrilateral: - //TODO pre-post rotate - GEOS_ERROR("Not Implemented Yet !"); - break; - // return computeVelocities_<4>(normals,fluxes, rotation); - case geos::ElementType::Tetrahedron: - return computeVelocities_<4>(normals,fluxes); - case geos::ElementType::Pyramid: - case geos::ElementType::Prism5: - case geos::ElementType::Wedge: - return computeVelocities_<5>(normals,fluxes); - case geos::ElementType::Hexahedron: - case geos::ElementType::Prism6: - return computeVelocities_<6>(normals,fluxes); - case geos::ElementType::Prism7: - return computeVelocities_<7>(normals,fluxes); - case geos::ElementType::Prism8: - return computeVelocities_<8>(normals,fluxes); - case geos::ElementType::Prism9: - return computeVelocities_<9>(normals,fluxes); - case geos::ElementType::Prism10: - return computeVelocities_<10>(normals,fluxes); - case geos::ElementType::Prism11: - return computeVelocities_<11>(normals,fluxes); - - - case geos::ElementType::Polygon: - case geos::ElementType::Polyhedron: - default: - GEOS_ERROR("Velocities are not computed on 2D Polygons cell type"); - break; - -} - - return {}; -} - - } /* namespace computationalGeometry */ } /* namespace geos */ diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 709618c75e0..92a741f13a6 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -286,66 +286,67 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, -template< int NFACES> -GEOS_HOST_DEVICE -array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, real64 const (&rotation)[3][3]) -{ - array2d velocity;//should be rotated so that the off plane is aligned with z ? - - //check compatibility on sizes - GEOS_ERROR_IF( (normals.size(1)!=fluxes.size()), GEOS_FMT("Error in parameters: normal({}) and fluxes({}) are different in sizes.", normals.size(1), fluxes.size()) ); - GEOS_ERROR_IF( (NFACES!=fluxes.size()), GEOS_FMT("Error in parameters: templatized size({}) and fluxes({}) are different in sizes.", NFACES, fluxes.size()) ); - - //TODO dispatch as a function of dims and nfaces - real64 velocity_c[2], fluxes_c[NFACES]; - real64 n[ NFACES ][ 2 ], ntn[ 2 ][ 2 ], ope[2][ NFACES ]; - - // N is nfaces x 3 matrices of normal - // (N^T N)^(-1)*N^T is the interpolant - real64 a = 0., b = 0., c = 0.; - for(int i=0; i +// GEOS_HOST_DEVICE +// array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, real64 const (&rotation)[3][3]) +// { +// array2d velocity;//should be rotated so that the off plane is aligned with z ? + +// //check compatibility on sizes +// GEOS_ERROR_IF( (normals.size(1)!=fluxes.size()), GEOS_FMT("Error in parameters: normal({}) and fluxes({}) are different in sizes.", normals.size(1), fluxes.size()) ); +// GEOS_ERROR_IF( (NFACES!=fluxes.size()), GEOS_FMT("Error in parameters: templatized size({}) and fluxes({}) are different in sizes.", NFACES, fluxes.size()) ); + +// //TODO dispatch as a function of dims and nfaces +// real64 velocity_c[2], fluxes_c[NFACES]; +// real64 n[ NFACES ][ 2 ], ntn[ 2 ][ 2 ], ope[2][ NFACES ]; + +// // N is nfaces x 3 matrices of normal +// // (N^T N)^(-1)*N^T is the interpolant +// real64 a = 0., b = 0., c = 0.; +// for(int i=0; i 1e-12); // else 2D most likely - // analytical inverse for symmetrical matrix - // [a c ] - // [. b ] - - ntn[0][0] = 1/detA * (b); - ntn[1][1] = 1/detA * (a); - ntn[0][1] = 1/detA * (-c); - ntn[1][0] = 1/detA * (-c); +// fluxes_c[i][0] = fluxes[i]; +// } + +// real64 detA = a*b - c*c; +// assert(LvArray::math::abs(detA) > 1e-12); // else 2D most likely +// // analytical inverse for symmetrical matrix +// // [a c ] +// // [. b ] + +// ntn[0][0] = 1/detA * (b); +// ntn[1][1] = 1/detA * (a); +// ntn[0][1] = 1/detA * (-c); +// ntn[1][0] = 1/detA * (-c); - //from Symmatrix - LvArray::tensorOps::Rij_add_AikBjk<2,NFACES>(ope, ntn, n); - LvArray::tensorOps::Ri_add_AijBj<2,NFACES>(velocity_c, ope, fluxes_c); +// //from Symmatrix +// LvArray::tensorOps::Rij_add_AikBjk<2,NFACES>(ope, ntn, n); +// LvArray::tensorOps::Ri_add_AijBj<2,NFACES>(velocity_c, ope, fluxes_c); - velocity[0][0] = velocity_c[0]; - velocity[0][1] = velocity_c[1]; - velocity[0][2] = 0.0;//off plane +// velocity[0][0] = velocity_c[0]; +// velocity[0][1] = velocity_c[1]; +// velocity[0][2] = 0.0;//off plane - return velocity; -} +// return velocity; +// } -template< int NFACES> +template< int NFACES, typename NORMAL_TYPE > GEOS_HOST_DEVICE -array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes) +array2d computeVelocities_( NORMAL_TYPE const & normals, arrayView1d< real64 const > const & fluxes) { array2d velocity; + velocity.resize(1,3); //check compatibility on sizes - GEOS_ERROR_IF( (normals.size(1)!=fluxes.size()), GEOS_FMT("Error in parameters: normal({}) and fluxes({}) are different in sizes.", normals.size(1), fluxes.size()) ); + // GEOS_ERROR_IF( (normals.size(1)!=fluxes.size()), GEOS_FMT("Error in parameters: normal({}) and fluxes({}) are different in sizes.", normals.size(1), fluxes.size()) ); GEOS_ERROR_IF( (NFACES!=fluxes.size()), GEOS_FMT("Error in parameters: templatized size({}) and fluxes({}) are different in sizes.", NFACES, fluxes.size()) ); //TODO dispatch as a function of dims and nfaces @@ -362,17 +363,29 @@ array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_ n[i][2] = normals[i][2]; a += normals[i][0]*normals[i][0]; - b += normals[i][1]*normals[i][1]; - c += normals[i][2]*normals[i][2]; + b += normals[i][1]*normals[i][1]; + c += normals[i][2]*normals[i][2]; - d += normals[i][0]*normals[i][1]; - f += normals[i][0]*normals[i][2]; - e += normals[i][1]*normals[i][2]; + d += normals[i][0]*normals[i][1]; + f += normals[i][0]*normals[i][2]; + e += normals[i][1]*normals[i][2]; - fluxes_c[i] = fluxes[i]; + fluxes_c[i] = fluxes[i]; + ope[0][i] = 0.; ope[1][i] = 0.; ope[2][i] = 0.; + + std::cout << i << " : " << fluxes[i] << "==" << fluxes_c[i] << std::endl; } real64 detA = a*b*c + 2*d*e*f - b*f*f - c*d*d - a*e*e; + std::cout << "n : \n" + << n[0][0] << " " << n[0][1] << " " << n[0][2] << "\n" + << n[1][0] << " " << n[1][1] << " " << n[1][2] << "\n" + << n[2][0] << " " << n[2][1] << " " << n[2][2] << "\n" + << n[3][0] << " " << n[3][1] << " " << n[3][2] << "\n" + << n[4][0] << " " << n[4][1] << " " << n[4][2] << "\n" + << n[5][0] << " " << n[5][1] << " " << n[5][2] << "\n"; + std::cout << "determinant " << detA << "\n"; + assert(LvArray::math::abs(detA) > 1e-12); // else 2D most likely // analytical inverse for symmetrical matrix // [a d f] @@ -391,10 +404,21 @@ array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_ ntn[2][0] = 1/detA * (d*e-f*b); ntn[2][1] = 1/detA * (a*e-f*d); + std::cout << "ntn : \n" + << ntn[0][0] << " " << ntn[0][1] << " " << ntn[0][2] << "\n" + << ntn[1][0] << " " << ntn[1][1] << " " << ntn[1][2] << "\n" + << ntn[2][0] << " " << ntn[2][1] << " " << ntn[2][2] << "\n"; + + //from Symmatrix LvArray::tensorOps::Rij_add_AikBjk<3,NFACES,3>(ope, ntn, n); LvArray::tensorOps::Ri_add_AijBj<3,NFACES>(velocity_c, ope, fluxes_c); + std::cout << "ope : \n" + << ope[0][0] << " " << ope[0][1] << " " << ope[0][2] << " " << ope[0][3]<< " " << ope[0][4]<< " " << ope[0][5]<< "\n" + << ope[1][0] << " " << ope[1][1] << " " << ope[1][2] << " " << ope[1][3]<< " " << ope[1][4]<< " " << ope[1][5]<< "\n" + << ope[2][0] << " " << ope[2][1] << " " << ope[2][2] << " " << ope[2][3]<< " " << ope[2][4]<< " " << ope[2][5]<< "\n"; + velocity[0][0] = velocity_c[0]; velocity[0][1] = velocity_c[1]; velocity[0][2] = velocity_c[2]; @@ -403,8 +427,53 @@ array2d computeVelocities_( arrayView2d< real64 const, nodes::REFERENCE_ } // Compute Velocity interpolations -array2d computeVelocity(arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & normals, arrayView1d< real64 const > const & fluxes, ElementType& elem /*subRegion.getElementType()*/); +template< typename NORMAL_TYPE > +array2d computeVelocity( NORMAL_TYPE const & normals, arrayView1d< real64 const > const & fluxes, ElementType& elem /*subRegion.getElementType()*/) +{ + +switch (elem) { + + case geos::ElementType::Triangle: + //TODO pre-post rotate + GEOS_ERROR("Not Implemented Yet !"); + break; + // return computeVelocities_<3>(normals,fluxes, rotation); + case geos::ElementType::Quadrilateral: + //TODO pre-post rotate + GEOS_ERROR("Not Implemented Yet !"); + break; + // return computeVelocities_<4>(normals,fluxes, rotation); + case geos::ElementType::Tetrahedron: + return computeVelocities_<4>(normals,fluxes); + case geos::ElementType::Pyramid: + case geos::ElementType::Prism5: + case geos::ElementType::Wedge: + return computeVelocities_<5>(normals,fluxes); + case geos::ElementType::Hexahedron: + case geos::ElementType::Prism6: + return computeVelocities_<6>(normals,fluxes); + case geos::ElementType::Prism7: + return computeVelocities_<7>(normals,fluxes); + case geos::ElementType::Prism8: + return computeVelocities_<8>(normals,fluxes); + case geos::ElementType::Prism9: + return computeVelocities_<9>(normals,fluxes); + case geos::ElementType::Prism10: + return computeVelocities_<10>(normals,fluxes); + case geos::ElementType::Prism11: + return computeVelocities_<11>(normals,fluxes); + + case geos::ElementType::Polygon: + case geos::ElementType::Polyhedron: + default: + GEOS_ERROR("Velocities are not computed on 2D Polygons cell type"); + break; + +} + + return {}; +} /** * @brief Change the orientation of the input vector to be consistent in a global sense. diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp index f3a948525a9..2e53d4c2fdb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp @@ -14,147 +14,139 @@ // Source includes -#include "mainInterface/initialization.hpp" +// #include "mainInterface/initialization.hpp" // #include "finiteVolume/CellElementStencilTPFA.hpp" -#include "physicsSolvers/fluidFlow/kernels/singlePhase/unitTests/testFlowKernelHelpers.hpp" +// #include "physicsSolvers/fluidFlow/kernels/singlePhase/unitTests/testFlowKernelHelpers.hpp" #include "mesh/utilities/ComputationalGeometry.hpp" // TPL includes +#include #include using namespace geos; +constexpr real64 EPS = 1e-16; -void setReferences( const real64 flux, const real64 dx, const real64 dy, real64 (& expectedVel)[3] ) +TEST(testIsoTAligned, Velocity_Tetra) { + constexpr int NFACES = 4; //2-hexa problem + constexpr int NTESTS = 2; + real64 tests[NTESTS][NFACES] = {{1.,1.,1.,1.}, + {2.,-1.,2.,-1.}}; + real64 expectedVel[NTESTS][3] = { + {0.,0.,0.}, + {0.,3.*LvArray::math::sqrt(3)/2.,0.} + }; + real64 const normals[NFACES][3] = { + {1./LvArray::math::sqrt(3), 1./LvArray::math::sqrt(3), 1./LvArray::math::sqrt(3)}, + {1./LvArray::math::sqrt(3), -1./LvArray::math::sqrt(3), -1./LvArray::math::sqrt(3)}, + {-1./LvArray::math::sqrt(3), 1./LvArray::math::sqrt(3), -1./LvArray::math::sqrt(3)}, + {-1./LvArray::math::sqrt(3), -1./LvArray::math::sqrt(3), 1./LvArray::math::sqrt(3)} + }; - expectedVel[0] = flux/dx; - expectedVel[1] = flux/dy; - expectedVel[2] = 0.0; +for(int itest=0; itest fluxes; + fluxes.resize(NFACES); + fluxes[0] = tests[itest][0]; + fluxes[1] = tests[itest][1]; + fluxes[2] = tests[itest][2]; + fluxes[3] = tests[itest][3]; -} + auto phaseVelocity = computationalGeometry::computeVelocities_(normals, fluxes.toViewConst()); + GEOS_LOG(GEOS_FMT("phaseVelocity {}", phaseVelocity)); + EXPECT_LT(phaseVelocity[0][0]-expectedVel[itest][0], EPS); + EXPECT_LT(phaseVelocity[0][1]-expectedVel[itest][1], EPS); + EXPECT_LT(phaseVelocity[0][2]-expectedVel[itest][2], EPS); + } -TEST( testAligned2D, Velocity_Hexa ) -{ - constexpr int nfaces = 6; //2-hexa problem - real64 const normals[nfaces][3] = {{1., 0., 0.}, - {-1., 0., 0.}, - {0., 1., 0.}, - {0., -1., 0.}, +} + +TEST( testRotated3D, Velocity_Hexa ) +{ + constexpr int NFACES = 6; //2-hexa problem + constexpr int NTESTS = 2; + real64 tests[NTESTS][NFACES] = {{1.,1.,1.,1.,1.,1.}, + {2.,-1.,2.,-1.,2.,-1.}}; + real64 expectedVel[NTESTS][3] = { + {0.,0.,0.}, + {0.,1.5*LvArray::math::sqrt(2),1.5} + }; + real64 const normals[NFACES][3] = {{1./LvArray::math::sqrt(2), 1./LvArray::math::sqrt(2), 0.}, + {-1./LvArray::math::sqrt(2), -1./LvArray::math::sqrt(2), 0.}, + {-1./LvArray::math::sqrt(2), 1./LvArray::math::sqrt(2), 0.}, + {1./LvArray::math::sqrt(2), -1./LvArray::math::sqrt(2), 0.}, {0., 0., 1.}, {0., 0., -1.} }; - real64 const fluxes[nfaces] = {1., 1., 1., 1., 1., 1.}; - - - const real64 flux = 0.001; - const int numElem = 1; - constexpr localIndex nPhases = 2; - constexpr localIndex nDirs = 3; - // const real64 phaseVelocity[numElem][nPhases][nDirs] = { - // {{0., 0., 0.}, {0., 0., 0.}}, - // }; - - ElementRegionManager::ElementViewAccessor< array2d< real64 > > normalsView = AccessorHelper< true >::makeElementAccessor< 2 >( normals[0][0], - 2, - seri[nfaces-1], - sesri[nfaces-1], - sei[nfaces-1], - nPhases, - nDirs ); - - ElementRegionManager::ElementViewAccessor< array1d< real64 > > fluxesView = AccessorHelper< true >::makeElementAccessor< 1 >( fluxes[0], - 2, - seri[nfaces-1], - sesri[nfaces-1], - sei[nfaces-1], - nPhases, - nDirs ); - - // ElementRegionManager::ElementViewAccessor< array3d< real64 > > phaseVelocityView = AccessorHelper< true >::makeElementAccessor< 3 >( phaseVelocity[0][0], - // tpfa.stencilSize( 0 ), - // seri[nfaces-1], - // sesri[nfaces-1], - // sei[nfaces-1], - // nPhases, - // nDirs ); - auto phaseVelocity = computationalGeometry::computeVelocities_(normalsView,fluxesView); - - // CellElementStencilTPFA tpfa; - // for( int kf = 0; kf < nfaces; ++kf ) - // { - // localIndex elementIndices[] = {ei_[kf][0], ei_[kf][1]}; - // tpfa.add( 2, elementRegionIndices, elementSubRegionIndices, elementIndices, weight, kf ); - // tpfa.addVectors( transMult[kf], geomStabSum, faceNormal[kf], cellToFaceVec[kf] ); - // } - - - // CellElementStencilTPFA::KernelWrapper wrapper = tpfa.createKernelWrapper(); - // array2d< real64 > globalCellDim( 2, 3 ); - // globalCellDim[0][0] = 2.; globalCellDim[1][0] = 2.; - // globalCellDim[0][1] = 2.; globalCellDim[1][1] = 2.; - // globalCellDim[0][2] = 0.; globalCellDim[1][2] = 0.; - // arrayView2d< const real64 > const globalCellDimView = globalCellDim.toViewConst(); - - - // CellElementStencilTPFA::IndexContainerViewConstType const & seri = tpfa.getElementRegionIndices(); - // CellElementStencilTPFA::IndexContainerViewConstType const & sesri = tpfa.getElementSubRegionIndices(); - // CellElementStencilTPFA::IndexContainerViewConstType const & sei = tpfa.getElementIndices(); - - - // ElementRegionManager::ElementViewAccessor< array3d< real64 > > phaseVelocityView = AccessorHelper< true >::makeElementAccessor< 3 >( phaseVelocity[0][0], - // tpfa.stencilSize( 0 ), - // seri[nfaces-1], - // sesri[nfaces-1], - // sei[nfaces-1], - // nPhases, - // nDirs ); - - - - // for( int iconn = 0; iconn < nfaces; ++iconn ) - // { - - - // wrapper.computeVelocity( iconn /*iconn*/, 0 /*ip*/, flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1 }, phaseVelocityView.toNestedView()); - // wrapper.computeVelocity( iconn /*iconn*/, 1 /*ip*/, 100*flux, {globalCellDimView[0], globalCellDimView[1]}, {-1, -1}, phaseVelocityView.toNestedView()); - // } - - - // for( int ip = 0; ip < 2; ++ip ) - // { - // auto fip = (99*ip+1)*flux; //shorthand to get the 1:100 ratio between phases as imposed above - // for( int ib = 0; ib < numElem; ++ib ) - // { - // real64 expectedVel[3]; - // setReferences( fip, globalCellDim[0][0], globalCellDim[0][1], expectedVel ); - // EXPECT_EQ( phaseVelocityView[0][0][ib][ip][0], expectedVel[0] ); - // EXPECT_EQ( phaseVelocityView[0][0][ib][ip][1], expectedVel[1] ); - // EXPECT_EQ( phaseVelocityView[0][0][ib][ip][2], expectedVel[2] ); - // } - // } -} +for(int itest=0; itest fluxes; + fluxes.resize(NFACES); + //{1., 1., 1., 1., 1., 1.}; + fluxes[0] = tests[itest][0]; + fluxes[1] = tests[itest][1]; + fluxes[2] = tests[itest][2]; + fluxes[3] = tests[itest][3]; + fluxes[4] = tests[itest][4]; + fluxes[5] = tests[itest][5]; + // real64 const fluxes[nfaces] = -int main( int argc, char * *argv ) -{ - ::testing::InitGoogleTest( &argc, argv ); + auto phaseVelocity = computationalGeometry::computeVelocities_(normals, fluxes.toViewConst()); + GEOS_LOG(GEOS_FMT("phaseVelocity {}", phaseVelocity)); + EXPECT_LT(phaseVelocity[0][0]-expectedVel[itest][0], EPS); + EXPECT_LT(phaseVelocity[0][1]-expectedVel[itest][1], EPS); + EXPECT_LT(phaseVelocity[0][2]-expectedVel[itest][2], EPS); + + } +} - geos::basicSetup( argc, argv ); +TEST( testAligned3D, Velocity_Hexa ) +{ + constexpr int NFACES = 6; //2-hexa problem + constexpr int NTESTS = 2; + real64 tests[NTESTS][NFACES] = {{1.,1.,1.,1.,1.,1.}, + {2.,-1.,2.,-1.,2.,-1.}}; + real64 expectedVel[NTESTS][3] = { + {0.,0.,0.}, + {1.5,1.5,1.5} + }; + + real64 const normals[NFACES][3] = + {{1., 0., 0.}, + {-1., 0., 0.}, + {0., 1., 0.}, + {0., -1., 0.}, + {0., 0., 1.}, + {0., 0., -1.} + };//ROT 0 - int const result = RUN_ALL_TESTS(); + for(int itest=0; itest fluxes; + fluxes.resize(NFACES); + //{1., 1., 1., 1., 1., 1.}; + fluxes[0] = tests[itest][0]; + fluxes[1] = tests[itest][1]; + fluxes[2] = tests[itest][2]; + fluxes[3] = tests[itest][3]; + fluxes[4] = tests[itest][4]; + fluxes[5] = tests[itest][5]; + // real64 const fluxes[nfaces] = - geos::basicCleanup(); + auto phaseVelocity = computationalGeometry::computeVelocities_(normals, fluxes.toViewConst()); + GEOS_LOG(GEOS_FMT("phaseVelocity {}", phaseVelocity)); + EXPECT_LT(phaseVelocity[0][0]-expectedVel[itest][0], EPS); + EXPECT_LT(phaseVelocity[0][1]-expectedVel[itest][1], EPS); + EXPECT_LT(phaseVelocity[0][2]-expectedVel[itest][2], EPS); + } - return result; -} \ No newline at end of file +}