Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: add analytical leakoff #3525

Open
wants to merge 5 commits into
base: develop
Choose a base branch
from
Open

feat: add analytical leakoff #3525

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
3931974
add carter's law
Guotong-Ren Jan 28, 2025
43fa250
Update HydrofractureSolver.cpp
paveltomin Jan 28, 2025
47ff846
Merge branch 'develop' into analytical_leakoff
paveltomin Jan 28, 2025
2566758
add a field for fractureCreationTime and some clean-ups
Guotong-Ren Jan 29, 2025
9cbee76
Merge branch 'analytical_leakoff' of https://github.com/GEOS-DEV/GEOS…
Guotong-Ren Jan 29, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,8 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies )
subRegion.registerField< fields::flow::hydraulicAperture >( getName() ).
setApplyDefaultValue( faceRegion.getDefaultAperture() );

subRegion.registerField< fields::flow::fractureCreationTime>( getName() ).
setApplyDefaultValue( 0.0 );
} );

FaceManager & faceManager = mesh.getFaceManager();
Expand Down
7 changes: 7 additions & 0 deletions src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -281,6 +281,13 @@ DECLARE_FIELD( massCreated,
WRITE_AND_READ,
"The amount of remaining mass that was introduced when the SurfaceElement was created." );

DECLARE_FIELD( fractureCreationTime,
"fractureCreationTime",
array1d< real64 >,
0,
LEVEL_1,
WRITE_AND_READ,
"The creation time for the fracture cell." );
}

}
Expand Down
84 changes: 81 additions & 3 deletions src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,12 @@ HydrofractureSolver< POROMECHANICS_SOLVER >::HydrofractureSolver( const string &
setApplyDefaultValue( 0 ).
setInputFlag( InputFlags::OPTIONAL ).
setDescription( "Flag to determine whether or not to apply lagging update for the fracture stencil weights. " );

registerWrapper( viewKeyStruct::leakoffConstString(), &m_leakoffCoefficient ).
setApplyDefaultValue( -1.0 ).
setInputFlag( InputFlags::OPTIONAL ).
setDescription( "Analytical leakoff coefficient." );


m_numResolves[0] = 0;
}
Expand Down Expand Up @@ -190,7 +196,7 @@ real64 HydrofractureSolver< POROMECHANICS_SOLVER >::fullyCoupledSolverStep( real
DomainPartition & domain )
{
// for initial fracture initialization in case when surface generator was called outside of the solver
initializeNewFractureFields( domain );
initializeNewFractureFields( time_n, domain );

real64 dtReturn = dt;

Expand Down Expand Up @@ -239,7 +245,7 @@ real64 HydrofractureSolver< POROMECHANICS_SOLVER >::fullyCoupledSolverStep( real
else
{
// We initialize the fields for new fracture cells
initializeNewFractureFields( domain );
initializeNewFractureFields( time_n, domain );

FieldIdentifiers fieldsToBeSync;

Expand Down Expand Up @@ -678,9 +684,77 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::assembleSystem( real64 const t

assembleFluidMassResidualDerivativeWrtDisplacement( domain, localMatrix );

if (m_leakoffCoefficient > -1.0 && !m_isMatrixPoroelastic)
assembleFluidLeakSource(time, dt, domain, dofManager, localMatrix, localRhs);

this->getRefDerivativeFluxResidual_dAperture()->zero();
}

template< typename POROMECHANICS_SOLVER>
void HydrofractureSolver< POROMECHANICS_SOLVER >::
assembleFluidLeakSource( double time,
double dt,
DomainPartition & domain,
DofManager const & dofManager,
CRSMatrixView< real64, globalIndex const > const & localMatrix,
arrayView1d< real64 > const & localRhs ) const
{

GEOS_MARK_FUNCTION;

string const presDofKey = dofManager.getKey( SinglePhaseBase::viewKeyStruct::elemDofFieldString() );

globalIndex const rankOffset = dofManager.rankOffset();

forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel const & mesh,
arrayView1d< string const > const & regionNames )
{
ElementRegionManager const & elemManager = mesh.getElemManager();

elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames,
[&]( localIndex const,
FaceElementSubRegion const & subRegion )
{
localIndex regionSize = subRegion.size();

string const & fluidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::fluidNamesString() );
SingleFluidBase const & fluid = this->template getConstitutiveModel< SingleFluidBase >( subRegion, fluidName );

arrayView1d< globalIndex const > const presDofNumber = subRegion.getReference< array1d< globalIndex > >( presDofKey );

arrayView2d< real64 const > const dens = fluid.density();
arrayView2d< real64 const > const dDens_dPressure = fluid.dDensity_dPressure();

arrayView1d< real64 const > const area = subRegion.getElementArea();
arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank();

arrayView1d< real64 const > const fractureCreationTime = subRegion.getField< fields::flow::fractureCreationTime>();

constexpr integer numDof = 1;
constexpr integer numEqn = 1;
globalIndex dofIndices[numDof]{};
real64 localJacobian[numEqn][numDof]{};
forAll< serialPolicy >( regionSize,
[&] GEOS_HOST_DEVICE ( localIndex const kfe )
{
if ( elemGhostRank[kfe] >= 0 || fabs(time - fractureCreationTime[kfe]) < 1e-12) return;
const globalIndex localRow = presDofNumber[kfe] - rankOffset;
for( integer idof = 0; idof < numDof; ++idof )
dofIndices[idof] = presDofNumber[kfe] + idof;
localJacobian[0][0] = m_leakoffCoefficient * area[kfe] * dDens_dPressure[kfe][0] / LvArray::math::sqrt(time + dt / 2.0 - fractureCreationTime[kfe]);
localMatrix.template addToRow< serialAtomic > ( localRow,
dofIndices,
localJacobian[0],
numDof );
// mid-point rule in time
localRhs[localRow] += m_leakoffCoefficient * area[kfe] * dens[kfe][0] / LvArray::math::sqrt(time + dt / 2.0 - fractureCreationTime[kfe]);
} );
} );
} );
}


template< typename POROMECHANICS_SOLVER >
void HydrofractureSolver< POROMECHANICS_SOLVER >::
assembleForceResidualDerivativeWrtPressure( DomainPartition & domain,
Expand Down Expand Up @@ -1002,7 +1076,8 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setUpDflux_dApertureMatrix( Do
}

template< typename POROMECHANICS_SOLVER >
void HydrofractureSolver< POROMECHANICS_SOLVER >::initializeNewFractureFields( DomainPartition & domain )
void HydrofractureSolver< POROMECHANICS_SOLVER >::initializeNewFractureFields( double time,
DomainPartition & domain )
{
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & meshLevel,
Expand Down Expand Up @@ -1038,6 +1113,7 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::initializeNewFractureFields( D
SingleFluidBase const & fluid = subRegion.getConstitutiveModel< SingleFluidBase >( fluidName );
real64 const defaultDensity = fluid.defaultDensity();
arrayView1d< real64 > const massCreated = subRegion.getField< fields::flow::massCreated >();
arrayView1d< real64 > const fractureCreationTime = subRegion.getField< fields::flow::fractureCreationTime>();


arrayView1d< real64 > const aperture = subRegion.getField< fields::elementAperture >();
Expand Down Expand Up @@ -1124,6 +1200,8 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::initializeNewFractureFields( D
aperture[newElemIndex] = 0;
}
}
// add creation time
fractureCreationTime[newElemIndex] = time;
GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator,
GEOS_FMT( "New elem index = {:4d} , init aper = {:4.2e}, init press = {:4.2e} ",
newElemIndex, aperture[newElemIndex], fluidPressure[newElemIndex] ) );
Expand Down
12 changes: 11 additions & 1 deletion src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,6 +179,8 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER

constexpr static char const * isLaggingFractureStencilWeightsUpdateString() { return "isLaggingFractureStencilWeightsUpdate"; }

constexpr static char const * leakoffConstString() {return "leakoffCoefficient"; }

#ifdef GEOS_USE_SEPARATION_COEFFICIENT
constexpr static char const * separationCoeff0String() { return "separationCoeff0"; }
constexpr static char const * apertureAtFailureString() { return "apertureAtFailure"; }
Expand Down Expand Up @@ -225,11 +227,17 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER
DomainPartition & domain ) override final;


void assembleFluidLeakSource( double time,
double dt,
DomainPartition & domain,
DofManager const & dofManager,
CRSMatrixView< real64, globalIndex const > const & localMatrix,
arrayView1d< real64 > const & localRhs ) const;
/**
* @brief Initialize fields on the newly created elements of the fracture.
* @param domain the physical domain object
*/
void initializeNewFractureFields( DomainPartition & domain );
void initializeNewFractureFields( double time, DomainPartition & domain );

// name of the contact relation
string m_contactRelationName;
Expand All @@ -256,6 +264,8 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER
// flag to determine whether or not to apply lagging update for the fracture stencil weights
integer m_isLaggingFractureStencilWeightsUpdate;

// analytical leakoff coefficient
real64 m_leakoffCoefficient;
};

ENUM_STRINGS( HydrofractureSolver< SinglePhasePoromechanics< SinglePhaseBase > >::InitializationType,
Expand Down
Loading