Skip to content

Fixed bug related to matching perforations to reservoir elements #1206

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

Merged
merged 3 commits into from
Jan 11, 2021
Merged

Fixed bug related to matching perforations to reservoir elements #1206

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
24ad3f4
fixed bug with perforations
francoishamon Oct 16, 2020
2145b4f
Merge branch 'develop' into bugfix/hamon/perforations
francoishamon Jan 7, 2021
46e0a62
removed dead code
francoishamon Jan 7, 2021
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
5 changes: 4 additions & 1 deletion src/coreComponents/mesh/WellElementRegion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,10 @@ void WellElementRegion::GenerateWell( MeshLevel & mesh,

globalIndex const matchedPerforations = MpiWrapper::Sum( perforationData->size() );
GEOSX_ERROR_IF( matchedPerforations != numPerforationsGlobal,
"Invalid mapping perforation-to-element in well " << this->getName() );
"Invalid mapping perforation-to-element in well " << this->getName() << "." <<
" This happens when GEOSX cannot match a perforation with a reservoir element." <<
" The most common reason for this error is that a perforation is on a section of " <<
" the well polyline located outside the domain." );


// 2) classify well elements based on connectivity to local mesh partition
Expand Down
28 changes: 10 additions & 18 deletions src/coreComponents/mesh/WellElementSubRegion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -538,7 +538,10 @@ void WellElementSubRegion::CheckPartitioningValidity( InternalWellGenerator cons
globalIndex const prevGlobal = prevElemIdsGlobal[iwelemGlobal][numBranches-1];

GEOSX_ERROR_IF( prevGlobal <= iwelemGlobal || prevGlobal < 0,
"Invalid partitioning in well " << getName() );
"The structure of well " << getName() << " is invalid. " <<
" The main reason for this error is that there may be no perforation" <<
" in the bottom well element of the well, which is required to have" <<
" a well-posed problem." );

if( elemStatusGlobal[prevGlobal] == WellElemStatus::LOCAL )
{
Expand Down Expand Up @@ -790,23 +793,12 @@ void WellElementSubRegion::ConnectPerforationsToMeshElements( MeshLevel & mesh,

// for each perforation, we have to find the reservoir element that contains the perforation

if( iperfLocal > 0 )
{
// get the info of the element matched with the previous perforation
// this will be used next to search around this reservoir element
erInit = m_perforationData.GetMeshElements().m_toElementRegion[iperfLocal-1];
esrInit = m_perforationData.GetMeshElements().m_toElementSubRegion[iperfLocal-1];
eiInit = m_perforationData.GetMeshElements().m_toElementIndex[iperfLocal-1];
}
else
{
// Step 1: first, we search for the reservoir element that is the *closest* from the center of well element
// note that this reservoir element does not necessarily contain the center of the well element
// this "init" reservoir element will be used in SearchLocalElements to find the reservoir element that
// contains the well element
InitializeLocalSearch( mesh, location,
erInit, esrInit, eiInit );
}
// Step 1: first, we search for the reservoir element that is the *closest* from the center of well element
// note that this reservoir element does not necessarily contain the center of the well element
// this "init" reservoir element will be used in SearchLocalElements to find the reservoir element that
// contains the well element
InitializeLocalSearch( mesh, location,
erInit, esrInit, eiInit );

// Step 2: then, search for the reservoir element that contains the well element
// to do that, we loop over the reservoir elements that are in the neighborhood of (erInit,esrInit,eiInit)
Expand Down
118 changes: 0 additions & 118 deletions src/coreComponents/meshUtilities/PerforationData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -225,84 +225,6 @@ void PerforationData::GetReservoirElementDimensions( MeshLevel const & mesh,

}


void PerforationData::ConnectToMeshElements( MeshLevel const & mesh,
InternalWellGenerator const & wellGeometry )
{
ElementRegionManager const * const elemManager = mesh.getElemManager();
NodeManager const * const nodeManager = mesh.getNodeManager();

ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > > elemCenter =
elemManager->ConstructArrayViewAccessor< real64, 2 >( ElementSubRegionBase::viewKeyStruct::elementCenterString );

arrayView2d< real64 const > const & perfCoordsGlobal = wellGeometry.GetPerfCoords();
arrayView1d< real64 const > const & perfTransGlobal = wellGeometry.GetPerfTransmissibility();

resize( perfCoordsGlobal.size( 0 ) );
localIndex iperfLocal = 0;

// loop over all the perforations
for( globalIndex iperfGlobal = 0; iperfGlobal < perfCoordsGlobal.size( 0 ); ++iperfGlobal )
{
R1Tensor const coord = LVARRAY_TENSOROPS_INIT_LOCAL_3( perfCoordsGlobal[iperfGlobal] );

// TODO actually trace coords
// TODO what if a fracture element is located

// find the closest reservoir element
auto ret = minLocOverElemsInMesh( &mesh, [&] ( localIndex const er,
localIndex const esr,
localIndex const ei ) -> real64
{
R1Tensor v = coord;
LvArray::tensorOps::subtract< 3 >( v, elemCenter[er][esr][ei] );
return LvArray::tensorOps::l2Norm< 3 >( v );
} );

// save the region, subregion and index
localIndex const er = std::get< 0 >( ret.second );
localIndex const esr = std::get< 1 >( ret.second );
localIndex const ei = std::get< 2 >( ret.second );

// a ghostRank check may be inserted if needed
// (not needed in the current implementation)

// check if perforation location is indeed inside the element
CellBlock const * const cellBlock = elemManager->GetRegion( er )->GetSubRegion< CellElementSubRegion >( esr );
array1d< array1d< localIndex > > faceNodes( cellBlock->numFacesPerElement() );

for( localIndex kf = 0; kf < cellBlock->numFacesPerElement(); ++kf )
{
cellBlock->GetFaceNodes( ei, kf, faceNodes[kf] );
}

if( !computationalGeometry::IsPointInsidePolyhedron( nodeManager->referencePosition(), faceNodes, coord ))
{
continue;
}

// TODO: what happens when the boundary is at the boundary of the MPI domain??

// now construct the local data

// store the indices of the mesh element
m_toMeshElements.m_toElementRegion [iperfLocal] = er;
m_toMeshElements.m_toElementSubRegion[iperfLocal] = esr;
m_toMeshElements.m_toElementIndex [iperfLocal] = ei;

// construct the local transmissibility and location maps
m_wellTransmissibility[iperfLocal] = perfTransGlobal[iperfGlobal];
LvArray::tensorOps::copy< 3 >( m_location[iperfLocal], coord );
m_localToGlobalMap[iperfLocal++] = iperfGlobal;
}

// set the size based on the number of perforations matched with local reservoir elements
resize( iperfLocal );
ConstructGlobalToLocalMap();

}


void PerforationData::ConnectToWellElements( InternalWellGenerator const & wellGeometry,
unordered_map< globalIndex, localIndex > const & globalToLocalWellElemMap,
globalIndex elemOffsetGlobal )
Expand All @@ -317,46 +239,6 @@ void PerforationData::ConnectToWellElements( InternalWellGenerator const & wellG
m_wellElementIndex[iperfLocal] = globalToLocalWellElemMap.at( ielemGlobal );
}

//DebugLocalPerforations();
}

void PerforationData::DebugLocalPerforations() const
{
if( size() == 0 )
{
return;
}

if( MpiWrapper::Comm_rank( MPI_COMM_GEOSX ) != 1 )
{
return;
}

std::cout << std::endl;
std::cout << "++++++++++++++++++++++++++" << std::endl;
std::cout << "PerforationData = " << getName() << " of " << getParent()->getName() << std::endl;
std::cout << "MPI rank = " << MpiWrapper::Comm_rank( MPI_COMM_GEOSX ) << std::endl;
std::cout << "Number of local perforations = " << size() << std::endl;

for( localIndex iperf = 0; iperf < size(); ++iperf )
{
std::cout << "m_toMeshElements.m_toElementRegion[" << iperf << "] = "
<< m_toMeshElements.m_toElementRegion[iperf]
<< std::endl;
std::cout << "m_toMeshElements.m_toElementSubRegion[" << iperf << "] = "
<< m_toMeshElements.m_toElementSubRegion[iperf]
<< std::endl;
std::cout << "m_toMeshElements.m_toElementIndex[" << iperf << "] = "
<< m_toMeshElements.m_toElementIndex[iperf]
<< std::endl;

std::cout << "m_wellElementIndexLocal[" << iperf << "] = " << m_wellElementIndex[iperf]
<< std::endl;
std::cout << "m_location[" << iperf << "] = " << m_location[iperf]
<< std::endl;
std::cout << "m_wellTransmissibility[" << iperf << "] = " << m_wellTransmissibility[iperf]
<< std::endl;
}
}

} //namespace geosx
12 changes: 0 additions & 12 deletions src/coreComponents/meshUtilities/PerforationData.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,15 +202,6 @@ class PerforationData : public ObjectManagerBase
*/
///@{

/**
* @brief Locate connected local mesh elements and resizes current object appropriately.
* @param[in] mesh target mesh level
* @param[in] wellGeometry InternalWellGenerator containing the global well topology
*/
void ConnectToMeshElements( MeshLevel const & mesh,
InternalWellGenerator const & wellGeometry );


/**
* @brief Connect each perforation to a local wellbore element.
* @param[in] wellGeometry InternalWellGenerator containing the global well topology
Expand Down Expand Up @@ -296,9 +287,6 @@ class PerforationData : public ObjectManagerBase

///@}

/// @cond DO_NOT_DOCUMENT
void DebugLocalPerforations() const;
/// @endcond

/// Global number of perforations
globalIndex m_numPerforationsGlobal;
Expand Down