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..c37bed37c14 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..01da35e8b65 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/mesh/utilities/ComputationalGeometry.cpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp index c920605de80..d4a7d386f2e 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.cpp @@ -23,6 +23,7 @@ namespace geos { namespace computationalGeometry -{} /* namespace computationalGeometry */ +{ +} /* namespace computationalGeometry */ } /* namespace geos */ diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index 61de4181531..92a741f13a6 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -284,6 +284,197 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, return area; } + + +// 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 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, typename NORMAL_TYPE > +GEOS_HOST_DEVICE +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( (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], fluxes_c[NFACES]; + real64 n[ NFACES ][ 3 ], ntn[ 3 ][ 3 ], ope[3][ 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., 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); + + 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]; + + return velocity; +} + +// Compute Velocity interpolations +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. * @tparam NORMAL_TYPE type of @p normal 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..523f40b4494 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 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..2e53d4c2fdb --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testVelocityReconstruction.cpp @@ -0,0 +1,152 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "physicsSolvers/fluidFlow/kernels/singlePhase/unitTests/testFlowKernelHelpers.hpp" +#include "mesh/utilities/ComputationalGeometry.hpp" + +// TPL includes +#include +#include + + +using namespace geos; +constexpr real64 EPS = 1e-16; + + +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)} + }; + + +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( 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.} + }; + +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] = + + 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( 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 + + + 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] = + + 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); + + } + +}