diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 352b2b972a..89ea94460e 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3006-5860-947a907 + baseline: integratedTests/baseline_integratedTests-pr3194-6060-e5ba0ce allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 7d44e61640..f2b0f19dec 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,18 +6,26 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3194 (2024-07-10) +====================== +Use aperture table in poromechanics with conforming fractures. Rebaseline the corresponding cases. + + PR #3006 (2024-07-01) ====================== Added baselines for new tests. Relaxing tolerances for singlePhasePoromechanics_FaultModel_smoke. + PR #3196 (2024-06-28) ====================== Added isLaggingFractureStencilWeightsUpdate to hydrofracture solve. Rebaseline because of the new input. + PR #3177 (2024-06-28) ====================== Added logLevel to TimeHistoryOutput. Rebaseline because of the new input flag. + PR #3181 (2024-06-25) ====================== Decouple debug matrix output from logLevel. Rebaseline because of the new input flag. diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp index 0ad5e5f865..d45ddd77e9 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp @@ -122,19 +122,20 @@ void ContactSolverBase::setFractureRegions( dataRepository::Group const & meshBo void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh, globalIndex & numStick, + globalIndex & numNewSlip, globalIndex & numSlip, globalIndex & numOpen ) const { ElementRegionManager const & elemManager = mesh.getElemManager(); - array1d< globalIndex > localCounter( 3 ); + array1d< globalIndex > localCounter( 4 ); elemManager.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion const & subRegion ) { arrayView1d< integer const > const & ghostRank = subRegion.ghostRank(); arrayView1d< integer const > const & fractureState = subRegion.getField< fields::contact::fractureState >(); - RAJA::ReduceSum< parallelHostReduce, localIndex > stickCount( 0 ), slipCount( 0 ), openCount( 0 ); + RAJA::ReduceSum< parallelHostReduce, localIndex > stickCount( 0 ), newSlipCount( 0 ), slipCount( 0 ), openCount( 0 ); forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { if( ghostRank[kfe] < 0 ) @@ -147,6 +148,10 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh, break; } case FractureState::NewSlip: + { + newSlipCount += 1; + break; + } case FractureState::Slip: { slipCount += 1; @@ -162,21 +167,23 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh, } ); localCounter[0] += stickCount.get(); - localCounter[1] += slipCount.get(); - localCounter[2] += openCount.get(); + localCounter[1] += newSlipCount.get(); + localCounter[2] += slipCount.get(); + localCounter[3] += openCount.get(); } ); - array1d< globalIndex > totalCounter( 3 ); + array1d< globalIndex > totalCounter( 4 ); MpiWrapper::allReduce( localCounter.data(), totalCounter.data(), - 3, + 4, MPI_SUM, MPI_COMM_GEOSX ); - numStick = totalCounter[0]; - numSlip = totalCounter[1]; - numOpen = totalCounter[2]; + numStick = totalCounter[0]; + numNewSlip = totalCounter[1]; + numSlip = totalCounter[2]; + numOpen = totalCounter[3]; } void ContactSolverBase::outputConfigurationStatistics( DomainPartition const & domain ) const @@ -184,6 +191,7 @@ void ContactSolverBase::outputConfigurationStatistics( DomainPartition const & d if( getLogLevel() >=1 ) { globalIndex numStick = 0; + globalIndex numNewSlip = 0; globalIndex numSlip = 0; globalIndex numOpen = 0; @@ -191,11 +199,11 @@ void ContactSolverBase::outputConfigurationStatistics( DomainPartition const & d MeshLevel const & mesh, arrayView1d< string const > const & ) { - computeFractureStateStatistics( mesh, numStick, numSlip, numOpen ); + computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen ); GEOS_LOG_RANK_0( GEOS_FMT( " Number of element for each fracture state:" - " stick: {:12} | slip: {:12} | open: {:12}", - numStick, numSlip, numOpen ) ); + " stick: {:12} | new slip: {:12} | slip: {:12} | open: {:12}", + numStick, numNewSlip, numSlip, numOpen ) ); } ); } } diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp index 16c16403a4..eaf0f35501 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp @@ -62,6 +62,7 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM void computeFractureStateStatistics( MeshLevel const & mesh, globalIndex & numStick, + globalIndex & numNewSlip, globalIndex & numSlip, globalIndex & numOpen ) const; diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp index 7c291b0f89..c313ef6570 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp @@ -406,8 +406,8 @@ struct ComputeRotationMatricesKernel forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { - localIndex const & f0 = elemsToFaces[k][0]; - localIndex const & f1 = elemsToFaces[k][1]; + localIndex const f0 = elemsToFaces[k][0]; + localIndex const f1 = elemsToFaces[k][1]; real64 Nbar[3]; Nbar[0] = faceNormal[f0][0] - faceNormal[f1][0]; diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp index 5418db4725..4e5bd37fc4 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp @@ -253,6 +253,13 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain { GEOS_MARK_FUNCTION; + real64 minNormalTractionTolerance( 1e10 ); + real64 maxNormalTractionTolerance( -1e10 ); + real64 minNormalDisplacementTolerance( 1e10 ); + real64 maxNormalDisplacementTolerance( -1e10 ); + real64 minSlidingTolerance( 1e10 ); + real64 maxSlidingTolerance( -1e10 ); + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & ) @@ -302,6 +309,13 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain arrayView1d< real64 > const & slidingTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::slidingToleranceString() ); + RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalTractionTolerance( 1e10 ); + RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalTractionTolerance( -1e10 ); + RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalDisplacementTolerance( 1e10 ); + RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalDisplacementTolerance( -1e10 ); + RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionSlidingTolerance( 1e10 ); + RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionSlidingTolerance( -1e10 ); + forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { if( ghostRank[kfe] < 0 ) @@ -385,15 +399,36 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain LvArray::tensorOps::scale< 3, 3 >( rotatedInvStiffApprox, area ); // Finally, compute tolerances for the given fracture element + normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus / 2.e+7; + minSubRegionNormalDisplacementTolerance.min( normalDisplacementTolerance[kfe] ); + maxSubRegionNormalDisplacementTolerance.max( normalDisplacementTolerance[kfe] ); + slidingTolerance[kfe] = sqrt( rotatedInvStiffApprox[ 1 ][ 1 ] * rotatedInvStiffApprox[ 1 ][ 1 ] + rotatedInvStiffApprox[ 2 ][ 2 ] * rotatedInvStiffApprox[ 2 ][ 2 ] ) * averageYoungModulus / 2.e+7; + minSubRegionSlidingTolerance.min( slidingTolerance[kfe] ); + maxSubRegionSlidingTolerance.max( slidingTolerance[kfe] ); + normalTractionTolerance[kfe] = 1.0 / 2.0 * averageConstrainedModulus / averageBoxSize0 * normalDisplacementTolerance[kfe]; + minSubRegionNormalTractionTolerance.min( normalTractionTolerance[kfe] ); + maxSubRegionNormalTractionTolerance.max( normalTractionTolerance[kfe] ); } } ); + + minNormalDisplacementTolerance = std::min( minNormalDisplacementTolerance, minSubRegionNormalDisplacementTolerance.get() ); + maxNormalDisplacementTolerance = std::max( maxNormalDisplacementTolerance, maxSubRegionNormalDisplacementTolerance.get() ); + minSlidingTolerance = std::min( minSlidingTolerance, minSubRegionSlidingTolerance.get() ); + maxSlidingTolerance = std::max( maxSlidingTolerance, maxSubRegionSlidingTolerance.get() ); + minNormalTractionTolerance = std::min( minNormalTractionTolerance, minSubRegionNormalTractionTolerance.get() ); + maxNormalTractionTolerance = std::max( maxNormalTractionTolerance, maxSubRegionNormalTractionTolerance.get() ); } } ); } ); + + GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{}: normal displacement tolerance = [{}, {}], sliding tolerance = [{}, {}], normal traction tolerance = [{}, {}]", + this->getName(), minNormalDisplacementTolerance, maxNormalDisplacementTolerance, + minSlidingTolerance, maxSlidingTolerance, + minNormalTractionTolerance, maxNormalTractionTolerance ) ); } void SolidMechanicsLagrangeContact::resetStateToBeginningOfStep( DomainPartition & domain ) @@ -586,7 +621,7 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time, assembleContact( domain, dofManager, localMatrix, localRhs ); - // for sequenatial: add (fixed) pressure force contribution into residual (no derivatives) + // for sequential: add (fixed) pressure force contribution into residual (no derivatives) if( m_isFixedStressPoromechanicsUpdate ) { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, @@ -2221,20 +2256,24 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai if( ghostRank[kfe] < 0 ) { integer const originalFractureState = fractureState[kfe]; - if( originalFractureState == contact::FractureState::Open ) + if( originalFractureState == FractureState::Open ) { - if( dispJump[kfe][0] > -normalDisplacementTolerance[kfe] ) + if( dispJump[kfe][0] <= -normalDisplacementTolerance[kfe] ) { - fractureState[kfe] = contact::FractureState::Open; - } - else - { - fractureState[kfe] = contact::FractureState::Stick; + fractureState[kfe] = FractureState::Stick; + if( getLogLevel() >= 10 ) + GEOS_LOG( GEOS_FMT( "{}: {} -> {}: dispJump = {}, normalDisplacementTolerance = {}", + kfe, originalFractureState, fractureState[kfe], + dispJump[kfe][0], normalDisplacementTolerance[kfe] ) ); } } else if( traction[kfe][0] > normalTractionTolerance[kfe] ) { - fractureState[kfe] = contact::FractureState::Open; + fractureState[kfe] = FractureState::Open; + if( getLogLevel() >= 10 ) + GEOS_LOG( GEOS_FMT( "{}: {} -> {}: traction = {}, normalTractionTolerance = {}", + kfe, originalFractureState, fractureState[kfe], + traction[kfe][0], normalTractionTolerance[kfe] ) ); } else { @@ -2245,29 +2284,33 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai contactWrapper.computeLimitTangentialTractionNorm( traction[kfe][0], dLimitTangentialTractionNorm_dTraction ); - if( originalFractureState == contact::FractureState::Stick && currentTau >= limitTau ) + if( originalFractureState == FractureState::Stick && currentTau >= limitTau ) { currentTau *= (1.0 - m_slidingCheckTolerance); } - else if( originalFractureState != contact::FractureState::Stick && currentTau <= limitTau ) + else if( originalFractureState != FractureState::Stick && currentTau <= limitTau ) { currentTau *= (1.0 + m_slidingCheckTolerance); } if( currentTau > limitTau ) { - if( originalFractureState == contact::FractureState::Stick ) + if( originalFractureState == FractureState::Stick ) { - fractureState[kfe] = contact::FractureState::NewSlip; + fractureState[kfe] = FractureState::NewSlip; } else { - fractureState[kfe] = contact::FractureState::Slip; + fractureState[kfe] = FractureState::Slip; } } else { - fractureState[kfe] = contact::FractureState::Stick; + fractureState[kfe] = FractureState::Stick; } + if( getLogLevel() >= 10 && fractureState[kfe] != originalFractureState ) + GEOS_LOG( GEOS_FMT( "{}: {} -> {}: currentTau = {}, limitTau = {}", + kfe, originalFractureState, fractureState[kfe], + currentTau, limitTau ) ); } changed += faceArea[kfe] * !compareFractureStates( originalFractureState, fractureState[kfe] ); @@ -2289,21 +2332,23 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai // and total area of fracture elements totalArea = MpiWrapper::sum( totalArea ); + GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( " {}: changed area {} out of {}", getName(), changedArea, totalArea ) ); + // Assume converged if changed area is below certain fraction of total area return changedArea <= m_nonlinearSolverParameters.m_configurationTolerance * totalArea; } bool SolidMechanicsLagrangeContact::isFractureAllInStickCondition( DomainPartition const & domain ) const { - globalIndex numStick, numSlip, numOpen; + globalIndex numStick, numNewSlip, numSlip, numOpen; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel const & mesh, arrayView1d< string const > const & ) { - computeFractureStateStatistics( mesh, numStick, numSlip, numOpen ); + computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen ); } ); - return ( ( numSlip + numOpen ) == 0 ); + return ( ( numNewSlip + numSlip + numOpen ) == 0 ); } real64 SolidMechanicsLagrangeContact::setNextDt( real64 const & currentDt, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index f117160c68..7a3ac6e99f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -2585,9 +2585,6 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; - if( m_keepFlowVariablesConstantDuringInitStep ) - return; - real64 maxDeltaPhaseVolFrac = 0.0; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 3c32bc821e..3c8554d6c4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -185,9 +185,6 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< fields::flow::hydraulicAperture >( getName() ). setApplyDefaultValue( faceRegion.getDefaultAperture() ); - subRegion.registerField< fields::flow::minimumHydraulicAperture >( getName() ). - setApplyDefaultValue( faceRegion.getDefaultAperture() ); - } ); FaceManager & faceManager = mesh.getFaceManager(); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp index 53f7ea8704..59b9975dde 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp @@ -168,14 +168,6 @@ DECLARE_FIELD( hydraulicAperture, WRITE_AND_READ, "Hydraulic aperture" ); -DECLARE_FIELD( minimumHydraulicAperture, - "minimumHydraulicAperture", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "minimum value of the hydraulic aperture" ); - DECLARE_FIELD( gravityCoefficient, "gravityCoefficient", array1d< real64 >, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 2c32fc1022..3fddc4300b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -750,12 +750,13 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, singlePhaseBaseKernels::StatisticsKernel:: saveDeltaPressure( subRegion.size(), pres, initPres, deltaPres ); - arrayView1d< real64 const > const dVol = subRegion.getField< fields::flow::deltaVolume >(); + 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; } ); SingleFluidBase const & fluid = @@ -803,7 +804,7 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, if( volume[ei] * density_n[ei][0] > 1.1 * creationMass[ei] ) { creationMass[ei] *= 0.75; - if( creationMass[ei]<1.0e-20 ) + if( creationMass[ei] < 1.0e-20 ) { creationMass[ei] = 0.0; } @@ -1255,9 +1256,6 @@ void SinglePhaseBase::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; - if( m_keepFlowVariablesConstantDuringInitStep ) - return; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp index 661318fdc4..e4591b67a1 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp @@ -19,6 +19,7 @@ #include "SinglePhasePoromechanicsConformingFractures.hpp" #include "constitutive/solid/PorousSolid.hpp" +#include "constitutive/contact/ContactSelector.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "linearAlgebra/solvers/BlockPreconditioner.hpp" #include "linearAlgebra/solvers/SeparateComponentPreconditioner.hpp" @@ -759,7 +760,6 @@ void SinglePhasePoromechanicsConformingFractures< FLOW_SOLVER >::updateHydraulic arrayView2d< real64 const > const fractureTraction = subRegion.getField< fields::contact::traction >(); arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >(); arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< fields::flow::aperture0 >(); - arrayView1d< real64 const > const minimumHydraulicAperture = subRegion.getField< flow::minimumHydraulicAperture >(); arrayView1d< real64 > const aperture = subRegion.getElementAperture(); arrayView1d< real64 > const hydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); @@ -768,25 +768,33 @@ void SinglePhasePoromechanicsConformingFractures< FLOW_SOLVER >::updateHydraulic string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( porousSolidName ); - constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) - { + string const & contactRelationName = subRegion.template getReference< string >( SolidMechanicsLagrangianFEM::viewKeyStruct::contactRelationNameString() ); + ContactBase const & contact = subRegion.getConstitutiveModel< ContactBase >( contactRelationName ); - typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); - - poromechanicsFracturesKernels::StateUpdateKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - porousMaterialWrapper, - dispJump, - pressure, - area, - volume, - deltaVolume, - aperture, - minimumHydraulicAperture, - oldHydraulicAperture, - hydraulicAperture, - fractureTraction ); + constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + { + using ContactType = TYPEOFREF( castedContact ); + typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + { + typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); + + poromechanicsFracturesKernels::StateUpdateKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + porousMaterialWrapper, + contactWrapper, + dispJump, + pressure, + area, + volume, + deltaVolume, + aperture, + oldHydraulicAperture, + hydraulicAperture, + fractureTraction ); + + } ); } ); } ); } ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp index 7c098a010d..c48eceed31 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp @@ -43,22 +43,21 @@ struct StateUpdateKernel * @param[in] volume the volume * @param[out] deltaVolume the change in volume * @param[out] aperture the aperture - * @param[in] minimumHydraulicAperture the * @param[in] oldHydraulicAperture the old hydraulic aperture * @param[out] hydraulicAperture the effecture aperture * @param[in] fractureTraction the fracture traction */ - template< typename POLICY, typename POROUS_WRAPPER > + template< typename POLICY, typename POROUS_WRAPPER, typename CONTACT_WRAPPER > static void launch( localIndex const size, POROUS_WRAPPER const & porousMaterialWrapper, + CONTACT_WRAPPER const & contactWrapper, 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 & minimumHydraulicAperture, arrayView1d< real64 const > const & oldHydraulicAperture, arrayView1d< real64 > const & hydraulicAperture, arrayView2d< real64 const > const & fractureTraction ) @@ -69,8 +68,8 @@ struct StateUpdateKernel // 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. - hydraulicAperture[k] = minimumHydraulicAperture[k] + aperture[k]; - real64 const dHydraulicAperture_dNormalJump = 1.0; + real64 dHydraulicAperture_dNormalJump = 0.0; + hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k], dHydraulicAperture_dNormalJump ); deltaVolume[k] = hydraulicAperture[k] * area[k] - volume[k]; diff --git a/src/coreComponents/schema/docs/Hydrofracture.rst b/src/coreComponents/schema/docs/Hydrofracture.rst index 250756df11..896a5cb77b 100644 --- a/src/coreComponents/schema/docs/Hydrofracture.rst +++ b/src/coreComponents/schema/docs/Hydrofracture.rst @@ -1,27 +1,32 @@ -============================= =================================================================================================================================== ======== ====================================================================================================================================================================================================================================================================================================================== -Name Type Default Description -============================= =================================================================================================================================== ======== ====================================================================================================================================================================================================================================================================================================================== -cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] -contactRelationName groupNameRef required Name of contact relation to enforce constraints on fracture boundary. -flowSolverName groupNameRef required Name of the flow solver used by the coupled solver -initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. -isMatrixPoroelastic integer 0 (no description available) -isThermal integer 0 Flag indicating whether the problem is thermal or not. Set isThermal="1" to enable the thermal coupling -logLevel integer 0 Log level -maxNumResolves integer 10 Value to indicate how many resolves may be executed to perform surface generation after the execution of flow and mechanics solver. -name groupName required A name is required for any non-unique nodes -newFractureInitializationType geos_HydrofractureSolver >_InitializationType Pressure Type of new fracture element initialization. Can be Pressure or Displacement. -solidSolverName groupNameRef required Name of the solid solver used by the coupled solver -stabilizationMultiplier real64 1 Constant multiplier of stabilization strength -stabilizationRegionNames groupNameRef_array {} Regions where stabilization is applied. -stabilizationType geos_stabilization_StabilizationType None | StabilizationType. Options are: - | None- Add no stabilization to mass equation - | Global- Add jump stabilization to all faces - | Local- Add jump stabilization on interior of macro elements -surfaceGeneratorName groupNameRef required Name of the surface generator to use in the hydrofracture solver -targetRegions groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. -useQuasiNewton integer 0 (no description available) -writeLinearSystem integer 0 Write matrix, rhs, solution to screen ( = 1) or file ( = 2). -LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` -NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` -============================= =================================================================================================================================== ======== ====================================================================================================================================================================================================================================================================================================================== \ No newline at end of file + + +===================================== =================================================================================================================================== ======== ====================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +===================================== =================================================================================================================================== ======== ====================================================================================================================================================================================================================================================================================================================== +cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +contactRelationName groupNameRef required Name of contact relation to enforce constraints on fracture boundary. +flowSolverName groupNameRef required Name of the flow solver used by the coupled solver +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +isLaggingFractureStencilWeightsUpdate integer 0 Flag to determine whether or not to apply lagging update for the fracture stencil weights. +isMatrixPoroelastic integer 0 (no description available) +isThermal integer 0 Flag indicating whether the problem is thermal or not. Set isThermal="1" to enable the thermal coupling +logLevel integer 0 Log level +maxNumResolves integer 10 Value to indicate how many resolves may be executed to perform surface generation after the execution of flow and mechanics solver. +name groupName required A name is required for any non-unique nodes +newFractureInitializationType geos_HydrofractureSolver >_InitializationType Pressure Type of new fracture element initialization. Can be Pressure or Displacement. +solidSolverName groupNameRef required Name of the solid solver used by the coupled solver +stabilizationMultiplier real64 1 Constant multiplier of stabilization strength +stabilizationRegionNames groupNameRef_array {} Regions where stabilization is applied. +stabilizationType geos_stabilization_StabilizationType None | StabilizationType. Options are: + | None- Add no stabilization to mass equation + | Global- Add jump stabilization to all faces + | Local- Add jump stabilization on interior of macro elements +surfaceGeneratorName groupNameRef required Name of the surface generator to use in the hydrofracture solver +targetRegions groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager. +useQuasiNewton integer 0 (no description available) +writeLinearSystem integer 0 Write matrix, rhs, solution to screen ( = 1) or file ( = 2). +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +===================================== =================================================================================================================================== ======== ====================================================================================================================================================================================================================================================================================================================== + +