diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 98934988d8..daca87f5ee 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -638,8 +638,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ { saveConvergedState( subRegion ); - arrayView1d< real64 > const & dVol = subRegion.template getField< fields::flow::deltaVolume >(); - dVol.zero(); + applyDeltaVolume( subRegion ); // This should fix NaN density in newly created fracture elements updatePorosityAndPermeability( subRegion ); @@ -703,14 +702,7 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, singlePhaseBaseKernels::StatisticsKernel:: saveDeltaPressure( subRegion.size(), pres, initPres, deltaPres ); - arrayView1d< real64 > const dVol = subRegion.getField< fields::flow::deltaVolume >(); - arrayView1d< real64 > const vol = subRegion.getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString() ); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) - { - vol[ei] += dVol[ei]; - dVol[ei] = 0.0; - } ); + applyDeltaVolume( subRegion ); SingleFluidBase const & fluid = getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); @@ -1350,4 +1342,15 @@ void SinglePhaseBase::saveConvergedState( ElementSubRegionBase & subRegion ) con mass_n.setValues< parallelDevicePolicy<> >( mass ); } +void SinglePhaseBase::applyDeltaVolume( ElementSubRegionBase & subRegion ) +{ + arrayView1d< real64 > const dVol = subRegion.template getField< fields::flow::deltaVolume >(); + arrayView1d< real64 > const vol = subRegion.template getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString()); + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + vol[ei] += dVol[ei]; + dVol[ei] = 0.0; + } ); +} + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index db72beaa75..944722d2bb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -391,6 +391,8 @@ class SinglePhaseBase : public FlowSolverBase */ virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override; + void applyDeltaVolume( ElementSubRegionBase & subRegion ); + /** * @brief Structure holding views into fluid properties used by the base solver. */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index 6c6e2be5bc..b9fc5986eb 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -26,6 +26,7 @@ #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" #include "finiteVolume/FluxApproximationBase.hpp" @@ -520,10 +521,10 @@ void SinglePhasePoromechanicsEmbeddedFractures::updateState( DomainPartition & d using HydraulicApertureModelType = TYPEOFREF( castedHydraulicApertureModel ); typename HydraulicApertureModelType::KernelWrapper hydraulicApertureModelWrapper = castedHydraulicApertureModel.createKernelWrapper(); - poromechanicsEFEMKernels::StateUpdateKernel:: + poromechanicsFracturesKernels::StateUpdateKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), - hydraulicApertureModelWrapper, porousMaterialWrapper, + hydraulicApertureModelWrapper, dispJump, pressure, area, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp index ac266ec21d..cba9e722fa 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp @@ -250,7 +250,6 @@ class SinglePhasePoromechanicsEFEM : /// The rank global densities arrayView2d< real64 const > const m_solidDensity; arrayView2d< real64 const > const m_fluidDensity; - arrayView2d< real64 const > const m_fluidDensity_n; arrayView2d< real64 const > const m_dFluidDensity_dPressure; /// The rank-global fluid pressure array. @@ -259,9 +258,6 @@ class SinglePhasePoromechanicsEFEM : /// The rank-global fluid pressure array. arrayView1d< real64 const > const m_fracturePressure; - /// The rank-global delta-fluid pressure array. - arrayView2d< real64 const > const m_porosity_n; - arrayView2d< real64 const > const m_tractionVec; arrayView3d< real64 const > const m_dTraction_dJump; @@ -284,6 +280,8 @@ class SinglePhasePoromechanicsEFEM : arrayView1d< real64 const > const m_deltaVolume; + arrayView1d< real64 const > const m_mass_n; + SortedArrayView< localIndex const > const m_fracturedElems; ArrayOfArraysView< localIndex const > const m_cellsToEmbeddedSurfaces; @@ -307,70 +305,6 @@ using SinglePhaseKernelFactory = finiteElement::KernelFactory< SinglePhasePorome real64 const (&)[3], string const >; -/** - * @brief A struct to perform volume, aperture and fracture traction updates - */ -struct StateUpdateKernel -{ - - /** - * @brief Launch the kernel function doing volume, aperture and fracture traction updates - * @tparam POLICY the type of policy used in the kernel launch - * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates - * @param[in] size the size of the subregion - * @param[in] contactWrapper the wrapper implementing the contact relationship - * @param[in] dispJump the displacement jump - * @param[in] pressure the pressure - * @param[in] area the area - * @param[in] volume the volume - * @param[out] deltaVolume the change in volume - * @param[out] aperture the aperture - * @param[out] hydraulicAperture the effecture aperture - * @param[out] fractureContactTraction the fracture contact traction - */ - template< typename POLICY, typename POROUS_WRAPPER, typename CONTACT_WRAPPER > - static void - launch( localIndex const size, - CONTACT_WRAPPER const & contactWrapper, - POROUS_WRAPPER const & porousMaterialWrapper, - arrayView2d< real64 const > const & dispJump, - arrayView1d< real64 const > const & pressure, - arrayView1d< real64 const > const & area, - arrayView1d< real64 const > const & volume, - arrayView1d< real64 > const & deltaVolume, - arrayView1d< real64 > const & aperture, - arrayView1d< real64 const > const & oldHydraulicAperture, - arrayView1d< real64 > const & hydraulicAperture, - arrayView2d< real64 > const & fractureEffectiveTraction ) - { - forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - // update aperture to be equal to the normal displacement jump - aperture[k] = dispJump[k][0]; // the first component of the jump is the normal one. - - real64 dHydraulicAperture_dNormalJump = 0.0; - real64 dHydraulicAperture_dNormalTraction = 0.0; - hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k], - fractureEffectiveTraction[k][0], - dHydraulicAperture_dNormalJump, - dHydraulicAperture_dNormalTraction ); - - deltaVolume[k] = hydraulicAperture[k] * area[k] - volume[k]; - - real64 const jump[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( dispJump[k] ); - real64 const effectiveTraction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fractureEffectiveTraction[k] ); - - // all perm update models below should need effective traction instead of total traction - // (total traction is combined forces of fluid pressure and effective traction) - porousMaterialWrapper.updateStateFromPressureApertureJumpAndTraction( k, 0, pressure[k], - oldHydraulicAperture[k], hydraulicAperture[k], - dHydraulicAperture_dNormalJump, - jump, effectiveTraction ); - - } ); - } -}; - } // namespace poromechanicsEFEMKernels } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp index efc8ae09ea..3be28e0804 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp @@ -77,12 +77,10 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, m_wDofNumber( jumpDofNumber ), m_solidDensity( inputConstitutiveType.getDensity() ), m_fluidDensity( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ), - m_fluidDensity_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ), m_dFluidDensity_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() ), m_matrixPressure( elementSubRegion.template getField< fields::flow::pressure >() ), m_fracturePressure( embeddedSurfSubRegion.template getField< fields::flow::pressure >() ), - m_porosity_n( inputConstitutiveType.getPorosity_n() ), m_tractionVec( embeddedSurfSubRegion.getField< fields::contact::traction >() ), m_dTraction_dJump( embeddedSurfSubRegion.getField< fields::contact::dTraction_dJump >() ), m_dTraction_dPressure( embeddedSurfSubRegion.getField< fields::contact::dTraction_dPressure >() ), @@ -94,6 +92,7 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, m_elementVolumeCell( elementSubRegion.getElementVolume() ), m_elementVolumeFrac( embeddedSurfSubRegion.getElementVolume() ), m_deltaVolume( embeddedSurfSubRegion.template getField< fields::flow::deltaVolume >() ), + m_mass_n( embeddedSurfSubRegion.template getField< fields::flow::mass_n >() ), m_fracturedElems( elementSubRegion.fracturedElementsList() ), m_cellsToEmbeddedSurfaces( elementSubRegion.embeddedSurfacesList().toViewConst() ), m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] }, @@ -322,9 +321,7 @@ complete( localIndex const k, // Mass balance accumulation real64 const newVolume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex ); - real64 const newMass = m_fluidDensity( embSurfIndex, 0 ) * newVolume; - real64 const oldMass = m_fluidDensity_n( embSurfIndex, 0 ) * m_elementVolumeFrac( embSurfIndex ); - real64 const localFlowResidual = ( newMass - oldMass ); + real64 const localFlowResidual = m_fluidDensity( embSurfIndex, 0 ) * newVolume - m_mass_n[embSurfIndex]; real64 const localFlowJumpJacobian = m_fluidDensity( embSurfIndex, 0 ) * m_surfaceArea[ embSurfIndex ]; real64 const localFlowFlowJacobian = m_dFluidDensity_dPressure( embSurfIndex, 0 ) * newVolume; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp index f996e34f0e..5b1d5c07d5 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp @@ -57,9 +57,7 @@ class ThermalSinglePhasePoromechanicsEFEM : using Base::m_matrixPresDofNumber; using Base::m_wDofNumber; using Base::m_fluidDensity; - using Base::m_fluidDensity_n; using Base::m_dFluidDensity_dPressure; - using Base::m_porosity_n; using Base::m_surfaceArea; using Base::m_elementVolumeFrac; using Base::m_deltaVolume; @@ -164,17 +162,18 @@ class ThermalSinglePhasePoromechanicsEFEM : arrayView2d< real64 const > const m_dFluidDensity_dTemperature; /// Views on fluid internal energy - arrayView2d< real64 const > const m_fluidInternalEnergy_n; arrayView2d< real64 const > const m_fluidInternalEnergy; arrayView2d< real64 const > const m_dFluidInternalEnergy_dPressure; arrayView2d< real64 const > const m_dFluidInternalEnergy_dTemperature; /// Views on temperature - arrayView1d< real64 const > const m_temperature_n; arrayView1d< real64 const > const m_temperature; /// The rank-global fluid pressure array. arrayView1d< real64 const > const m_matrixTemperature; + + /// Views on energy + arrayView1d< real64 const > const m_energy_n; }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp index c538093411..883ed5e6cc 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp @@ -68,15 +68,14 @@ ThermalSinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, fluidModelKey ), m_dFluidDensity_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dTemperature() ), - m_fluidInternalEnergy_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy_n() ), m_fluidInternalEnergy( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ), m_dFluidInternalEnergy_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dInternalEnergy_dPressure() ), m_dFluidInternalEnergy_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dInternalEnergy_dTemperature() ), - m_temperature_n( embeddedSurfSubRegion.template getField< fields::flow::temperature_n >() ), m_temperature( embeddedSurfSubRegion.template getField< fields::flow::temperature >() ), - m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() ) + m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() ), + m_energy_n( embeddedSurfSubRegion.template getField< fields::flow::energy_n >() ) {} @@ -176,9 +175,8 @@ complete( localIndex const k, localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0]; // Energy balance accumulation real64 const volume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex ); - real64 const volume_n = m_elementVolumeFrac( embSurfIndex ); real64 const fluidEnergy = m_fluidDensity( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * volume; - real64 const fluidEnergy_n = m_fluidDensity_n( embSurfIndex, 0 ) * m_fluidInternalEnergy_n( embSurfIndex, 0 ) * volume_n; + real64 const fluidEnergy_n = m_energy_n[embSurfIndex]; stack.dFluidMassIncrement_dTemperature = m_dFluidDensity_dTemperature( embSurfIndex, 0 ) * volume;