Skip to content

Commit

Permalink
Initialization of fields in new fracture cells in HydroFracture solver (
Browse files Browse the repository at this point in the history 
#2960)
  • Loading branch information
@frankfeifan
frankfeifan authored Mar 20, 2024
1 parent 9a94d80 commit ce219dc
Show file tree
Hide file tree
Showing 82 changed files with 1,423 additions and 1,491 deletions.
2 changes: 1 addition & 1 deletion integratedTests
Submodule integratedTests updated 150 files
4 changes: 1 addition & 3 deletions src/coreComponents/finiteVolume/FluxApproximationBase.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,11 +132,9 @@ class FluxApproximationBase : public dataRepository::Group
* @brief Add a new fracture stencil.
* @param[in,out] mesh the mesh on which to add the fracture stencil
* @param[in] faceElementRegionName the face element region name
* @param[in] initFields if true initialize physical fields, like pressure
*/
virtual void addToFractureStencil( MeshLevel & mesh,
string const & faceElementRegionName,
bool const initFields ) const = 0;
string const & faceElementRegionName ) const = 0;

/**
* @brief Add a new embedded fracture stencil.
Expand Down
19 changes: 0 additions & 19 deletions src/coreComponents/finiteVolume/SurfaceElementStencil.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,25 +24,6 @@
namespace geos
{

/// @cond DO_NOT_DOCUMENT
// TODO remove! This option allows for the creation of new mass inside a newly
// created FaceElement. The new mass will be equal to:
// creationMass = defaultDensity * defaultAperture * faceArea.
// If 0, then the beginning of step density is artificially set to zero...which
// may cause some newton convergence problems.
#define ALLOW_CREATION_MASS 1


// TODO remove! This option sets the pressure in a newly created FaceElement to
// be the lowest value of all attached non-new FaceElements.
#define SET_CREATION_PRESSURE 1

// TODO remove! This option sets the nodal displacements attached a newly
// created FaceElement to some scalar fraction of the aperture of the
// lowest attached non-new FaceElements.
#define SET_CREATION_DISPLACEMENT 0
/// @endcond

/**
* @brief Describes properties of SurfaceElementStencil.
*
Expand Down
244 changes: 1 addition & 243 deletions src/coreComponents/finiteVolume/TwoPointFluxApproximation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -282,11 +282,6 @@ void TwoPointFluxApproximation::addFractureFractureConnectionsDFM( MeshLevel & m
hydraulicAperture,
fractureRegionIndex,
elemGhostRank,
#if SET_CREATION_DISPLACEMENT==1
faceToNodesMap,
totalDisplacement,
aperture,
#endif
&fractureStencil]
( localIndex const k )
{
Expand Down Expand Up @@ -365,237 +360,6 @@ void TwoPointFluxApproximation::addFractureFractureConnectionsDFM( MeshLevel & m
} );
}

void TwoPointFluxApproximation::initNewFractureFieldsDFM( MeshLevel & mesh,
string const & faceElementRegionName ) const
{
// TODO Note that all of this initialization should be performed elsewhere.
// This is just here because it was convenient, but it is not appropriate
// to have physics based initialization in the flux approximator.

#if !defined(SET_CREATION_DISPLACEMENT)
static_assert( true, "must have SET_CREATION_DISPLACEMENT defined" );
#endif

#if !defined(ALLOW_CREATION_MASS)
static_assert( true, "must have ALLOW_CREATION_MASS defined" );
#endif

#if !defined(SET_CREATION_PRESSURE)
static_assert( true, "must have SET_CREATION_PRESSURE defined" );
#endif

ElementRegionManager & elemManager = mesh.getElemManager();
ElementRegionManager::ElementViewAccessor< arrayView1d< integer const > > const elemGhostRank =
elemManager.constructArrayViewAccessor< integer, 1 >( ObjectManagerBase::viewKeyStruct::ghostRankString() );

SurfaceElementRegion & fractureRegion = elemManager.getRegion< SurfaceElementRegion >( faceElementRegionName );
localIndex const fractureRegionIndex = fractureRegion.getIndexInParent();
FaceElementSubRegion & fractureSubRegion = fractureRegion.getUniqueSubRegion< FaceElementSubRegion >();
ArrayOfArraysView< localIndex const > const & fractureConnectorsToFaceElements = fractureSubRegion.m_2dFaceTo2dElems.toViewConst();
#if SET_CREATION_DISPLACEMENT==1
NodeManager & nodeManager = mesh.getNodeManager();
FaceManager const & faceManager = mesh.getFaceManager();
ArrayOfArraysView< localIndex const > const & faceToNodesMap = faceManager.nodeList();
FaceElementSubRegion::FaceMapType const & faceMap = fractureSubRegion.faceList();
array2dLayoutIncrDisplacementConst const incrementalDisplacement =
nodeManager.getField< fields::solidMechanics::incrementalDisplacement >();
array2dLayoutTotalDisplacementConst const totalDisplacement =
nodeManager.getField< fields::solidMechanics::totalDisplacement >();
arrayView1d< real64 > const aperture = fractureSubRegion.getReference< array1d< real64 > >( "elementAperture" );
#endif

#ifdef GEOSX_USE_SEPARATION_COEFFICIENT
arrayView1d< real64 > const apertureF = fractureSubRegion.getReference< array1d< real64 > >( "apertureAtFailure" );
#endif

#if ALLOW_CREATION_MASS==0
arrayView1d< real64 > const dens = fractureSubRegion.getReference< array1d< real64 > >( "density_n" );
#endif

#if SET_CREATION_PRESSURE==1
arrayView1d< real64 > const fluidPressure_n = fractureSubRegion.getField< fields::flow::pressure_n >();
arrayView1d< real64 > const fluidPressure = fractureSubRegion.getField< fields::flow::pressure >();
// Set the new face elements to some unphysical numbers to make sure they get set by the following routines.
SortedArrayView< localIndex const > const newFaceElements = fractureSubRegion.m_newFaceElements.toViewConst();

forAll< serialPolicy >( fractureSubRegion.m_newFaceElements.size(), [=]( localIndex const k )
{
localIndex const kfe = newFaceElements[k];
fluidPressure[kfe] = 1.0e99;
#ifdef GEOSX_USE_SEPARATION_COEFFICIENT
apertureF[kfe] = aperture[kfe];
#endif
#if SET_CREATION_DISPLACEMENT==1
aperture[kfe] = 1.0e99;
#endif
} );

#endif // SET_CREATION_PRESSURE

SortedArray< localIndex > allNewElems;
allNewElems.insert( fractureSubRegion.m_newFaceElements.begin(),
fractureSubRegion.m_newFaceElements.end() );
SortedArrayView< localIndex const > const recalculateFractureConnectorEdges = fractureSubRegion.m_recalculateConnectionsFor2dFaces.toViewConst();

// add new connectors/connections between face elements to the fracture stencil
forAll< serialPolicy >( recalculateFractureConnectorEdges.size(),
[ &allNewElems,
recalculateFractureConnectorEdges,
fractureConnectorsToFaceElements,
fractureRegionIndex,
elemGhostRank,
fluidPressure,
fluidPressure_n,
#if SET_CREATION_DISPLACEMENT==1
faceToNodesMap,
totalDisplacement,
aperture,
#endif
&fractureSubRegion
]
( localIndex const k )
{
localIndex const fci = recalculateFractureConnectorEdges[k];
localIndex const numElems = fractureConnectorsToFaceElements.sizeOfArray( fci );
#if SET_CREATION_PRESSURE==1
real64 initialPressure = 1.0e99;
#endif
#if SET_CREATION_DISPLACEMENT==1
real64 initialAperture = 1.0e99;
#endif
SortedArray< localIndex > newElems;
bool containsLocalElement = false;

// loop over all face elements attached to the connector and add them to the stencil
for( localIndex kfe=0; kfe<numElems; ++kfe )
{
localIndex const fractureElementIndex = fractureConnectorsToFaceElements[fci][kfe];
containsLocalElement = containsLocalElement || elemGhostRank[fractureRegionIndex][0][fractureElementIndex] < 0;

// code to initialize new face elements with pressures from neighbors
if( fractureSubRegion.m_newFaceElements.count( fractureElementIndex ) == 0 )
{
#if SET_CREATION_PRESSURE==1
initialPressure = std::min( initialPressure, fluidPressure_n[fractureElementIndex] );
#endif
#if SET_CREATION_DISPLACEMENT==1
initialAperture = std::min( initialAperture, aperture[fractureElementIndex] );
#endif
}
else
{
newElems.insert( fractureElementIndex );
allNewElems.insert( fractureElementIndex );
}
}

if( !containsLocalElement )
{
return;
}

// loop over new face elements attached to this connector
for( localIndex const newElemIndex : newElems )
{
#if SET_CREATION_PRESSURE==1
fluidPressure[newElemIndex] = std::min( fluidPressure[newElemIndex], initialPressure );
#endif
#if SET_CREATION_DISPLACEMENT==1
// Set the aperture/fluid pressure for the new face element to be the minimum
// of the existing value, smallest aperture/pressure from a connected face element.
//aperture[newElemIndex] = std::min(aperture[newElemIndex], initialAperture);

localIndex const faceIndex0 = faceMap( newElemIndex, 0 );
localIndex const faceIndex1 = faceMap( newElemIndex, 1 );

localIndex const numNodesPerFace = faceToNodesMap.sizeOfArray( faceIndex0 );

bool zeroDisp = true;

for( localIndex a=0; a<numNodesPerFace; ++a )
{
localIndex const node0 = faceToNodesMap( faceIndex0, a );
localIndex const node1 = faceToNodesMap( faceIndex1, a==0 ? a : numNodesPerFace-a );
if( LvArray::math::abs( LvArray::tensorOps::l2Norm< 3 >( totalDisplacement[node0] ) ) > 1.0e-99 &&
LvArray::math::abs( LvArray::tensorOps::l2Norm< 3 >( totalDisplacement[node1] ) ) > 1.0e-99 )
{
zeroDisp = false;
}
}
if( zeroDisp )
{
aperture[newElemIndex] = 0;
}
#endif
}
} );

SortedArray< localIndex > touchedNodes;
forAll< serialPolicy >( allNewElems.size(),
[ &allNewElems
, fluidPressure
, fluidPressure_n
#if SET_CREATION_DISPLACEMENT==1
, aperture
, faceMap
, faceNormal
, faceToNodesMap
, &touchedNodes
, incrementalDisplacement
, totalDisplacement
, this
#endif
]( localIndex const k )
{
localIndex const newElemIndex = allNewElems[k];
// if the value of pressure was not set, then set it to zero and punt.
if( fluidPressure[newElemIndex] > 1.0e98 )
{
fluidPressure[newElemIndex] = 0.0;
}
fluidPressure_n[newElemIndex] = fluidPressure[newElemIndex];
#if ALLOW_CREATION_MASS==0
// set the initial density of the face element to 0 to enforce mass conservation ( i.e. no creation of mass)
dens[newElemIndex] = 0.0;
#endif
#if SET_CREATION_DISPLACEMENT==1
// If the aperture has been set, then we can set the estimate of displacements.
if( aperture[newElemIndex] < 1e98 )
{
localIndex const faceIndex0 = faceMap( newElemIndex, 0 );
localIndex const faceIndex1 = faceMap( newElemIndex, 1 );

real64 newDisp[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( faceNormal[ faceIndex0 ] );
LvArray::tensorOps::scale< 3 >( newDisp, -aperture[newElemIndex] );
localIndex const numNodesPerFace = faceToNodesMap.sizeOfArray( faceIndex0 );
for( localIndex a=0; a<numNodesPerFace; ++a )
{
localIndex const node0 = faceToNodesMap( faceIndex0, a );
localIndex const node1 = faceToNodesMap( faceIndex1, a==0 ? a : numNodesPerFace-a );

touchedNodes.insert( node0 );
touchedNodes.insert( node1 );

if( node0 != node1 && touchedNodes.count( node0 )==0 )
{
LvArray::tensorOps::add< 3 >( incrementalDisplacement[node0], newDisp );
LvArray::tensorOps::add< 3 >( totalDisplacement[node0], newDisp );
LvArray::tensorOps::subtract< 3 >( incrementalDisplacement[node1], newDisp );
LvArray::tensorOps::subtract< 3 >( totalDisplacement[node1], newDisp );
}
}
}
if( this->getLogLevel() > 1 )
{
printf( "New elem index, init aper, init press = %4ld, %4.2e, %4.2e \n",
newElemIndex,
aperture[newElemIndex],
fluidPressure[newElemIndex] );
}
#endif
} );
}

void TwoPointFluxApproximation::cleanMatrixMatrixConnectionsDFM( MeshLevel & mesh,
string const & faceElementRegionName ) const
{
Expand Down Expand Up @@ -761,17 +525,11 @@ void TwoPointFluxApproximation::addFractureMatrixConnectionsDFM( MeshLevel & mes
}

void TwoPointFluxApproximation::addToFractureStencil( MeshLevel & mesh,
string const & faceElementRegionName,
bool const initFields ) const
string const & faceElementRegionName ) const
{
this->addFractureFractureConnectionsDFM( mesh, faceElementRegionName );
this->cleanMatrixMatrixConnectionsDFM( mesh, faceElementRegionName );
this->addFractureMatrixConnectionsDFM( mesh, faceElementRegionName );

if( initFields )
{
this->initNewFractureFieldsDFM( mesh, faceElementRegionName );
}
}

void TwoPointFluxApproximation::addFractureMatrixConnectionsEDFM( MeshLevel & mesh,
Expand Down
3 changes: 1 addition & 2 deletions src/coreComponents/finiteVolume/TwoPointFluxApproximation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,8 +73,7 @@ class TwoPointFluxApproximation : public FluxApproximationBase
virtual void registerFractureStencil( Group & stencilGroup ) const override;

virtual void addToFractureStencil( MeshLevel & mesh,
string const & faceElementRegionName,
bool const initFields ) const override;
string const & faceElementRegionName ) const override;

virtual void registerBoundaryStencil( Group & stencilGroup,
string const & setName ) const override;
Expand Down
12 changes: 1 addition & 11 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -359,14 +359,8 @@ class SurfaceElementBasedAssemblyKernel : public ElementBasedAssemblyKernel< Sur
CRSMatrixView< real64, globalIndex const > const & localMatrix,
arrayView1d< real64 > const & localRhs )
: Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs )
#if ALLOW_CREATION_MASS
, m_creationMass( subRegion.getReference< array1d< real64 > >( SurfaceElementSubRegion::viewKeyStruct::creationMassString() ) )
#endif
{
#if !defined(ALLOW_CREATION_MASS)
static_assert( true, "must have ALLOW_CREATION_MASS defined" );
#endif
}
{}

/**
* @brief Compute the local accumulation contributions to the residual and Jacobian
Expand All @@ -380,20 +374,16 @@ class SurfaceElementBasedAssemblyKernel : public ElementBasedAssemblyKernel< Sur
{
Base::computeAccumulation( ei, stack, [&] ()
{
#if ALLOW_CREATION_MASS
if( Base::m_volume[ei] * Base::m_density_n[ei][0] > 1.1 * m_creationMass[ei] )
{
stack.localResidual[0] += m_creationMass[ei] * 0.25;
}
#endif
} );
}

protected:

#if ALLOW_CREATION_MASS
arrayView1d< real64 const > const m_creationMass;
#endif

};

Expand Down
17 changes: 3 additions & 14 deletions src/coreComponents/physicsSolvers/fluidFlow/ThermalSinglePhaseBaseKernels.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -337,15 +337,9 @@ class SurfaceElementBasedAssemblyKernel : public ElementBasedAssemblyKernel< Sur
constitutive::CoupledSolidBase const & solid,
CRSMatrixView< real64, globalIndex const > const & localMatrix,
arrayView1d< real64 > const & localRhs )
: Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs )
#if ALLOW_CREATION_MASS
, m_creationMass( subRegion.getReference< array1d< real64 > >( SurfaceElementSubRegion::viewKeyStruct::creationMassString() ) )
#endif
{
#if !defined(ALLOW_CREATION_MASS)
static_assert( true, "must have ALLOW_CREATION_MASS defined" );
#endif
}
: Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ),
m_creationMass( subRegion.getReference< array1d< real64 > >( SurfaceElementSubRegion::viewKeyStruct::creationMassString() ) )
{}

/**
* @brief Compute the local accumulation contributions to the residual and Jacobian
Expand All @@ -358,20 +352,15 @@ class SurfaceElementBasedAssemblyKernel : public ElementBasedAssemblyKernel< Sur
Base::StackVariables & stack ) const
{
Base::computeAccumulation( ei, stack );

#if ALLOW_CREATION_MASS
if( Base::m_volume[ei] * Base::m_density_n[ei][0] > 1.1 * m_creationMass[ei] )
{
stack.localResidual[0] += m_creationMass[ei] * 0.25;
}
#endif
}

protected:

#if ALLOW_CREATION_MASS
arrayView1d< real64 const > const m_creationMass;
#endif

};

Expand Down
Loading

0 comments on commit ce219dc

Please sign in to comment.