From 406a9f31bda0219cd9b191eb4a318ac73983decc Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 26 Jan 2024 14:09:08 -0600 Subject: [PATCH 1/5] Fix the black oil case by tightening the nonlinear tolerance (#2955) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Pavel Tomin <“paveltomin@users.noreply.github.com”> --- .../black_oil_wells_unsaturated_3d.xml | 2 +- integratedTests | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/inputFiles/compositionalMultiphaseWell/black_oil_wells_unsaturated_3d.xml b/inputFiles/compositionalMultiphaseWell/black_oil_wells_unsaturated_3d.xml index 97e883224e8..142561d196e 100644 --- a/inputFiles/compositionalMultiphaseWell/black_oil_wells_unsaturated_3d.xml +++ b/inputFiles/compositionalMultiphaseWell/black_oil_wells_unsaturated_3d.xml @@ -12,7 +12,7 @@ targetRegions="{ region, wellRegion1, wellRegion2 }"> diff --git a/integratedTests b/integratedTests index 933508f06ce..402c5a9f7ab 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit 933508f06ce8290168af6c9f919e81e567e78c40 +Subproject commit 402c5a9f7ab2f3c5f123603c9db93e3feed51c12 From 3eea6dec610f4ba764906139bda48f245c23b7f9 Mon Sep 17 00:00:00 2001 From: Lionel Untereiner Date: Fri, 26 Jan 2024 21:09:52 +0100 Subject: [PATCH 2/5] Add a VTK well import (#2551) --- src/coreComponents/mesh/CMakeLists.txt | 11 +- .../mesh/generators/InternalWellGenerator.cpp | 534 +----------------- .../mesh/generators/InternalWellGenerator.hpp | 289 +--------- .../mesh/generators/VTKUtilities.cpp | 22 +- .../mesh/generators/VTKWellGenerator.cpp | 105 ++++ .../mesh/generators/VTKWellGenerator.hpp | 88 +++ .../mesh/generators/WellGeneratorABC.hpp | 204 +++++++ .../mesh/generators/WellGeneratorBase.cpp | 531 ++++++++++++++++- .../mesh/generators/WellGeneratorBase.hpp | 218 ++++++- .../schema/docs/AcousticElasticSEM.rst | 18 + .../schema/docs/AcousticElasticSEM_other.rst | 13 + .../schema/docs/InternalMesh.rst | 1 + .../schema/docs/InternalMesh_other.rst | 1 + .../schema/docs/InternalWell.rst | 4 +- .../schema/docs/InternalWellbore.rst | 1 + .../schema/docs/InternalWellbore_other.rst | 1 + .../docs/PorousDamageElasticIsotropic.rst | 13 + .../PorousDamageElasticIsotropic_other.rst | 9 + .../PorousDamageSpectralElasticIsotropic.rst | 13 + ...usDamageSpectralElasticIsotropic_other.rst | 9 + .../PorousDamageVolDevElasticIsotropic.rst | 13 + ...rousDamageVolDevElasticIsotropic_other.rst | 9 + src/coreComponents/schema/docs/VTKMesh.rst | 1 + .../schema/docs/VTKMesh_other.rst | 1 + src/coreComponents/schema/docs/VTKWell.rst | 17 + .../schema/docs/VTKWell_other.rst | 9 + src/coreComponents/schema/schema.xsd | 55 +- src/coreComponents/schema/schema.xsd.other | 8 + .../unitTests/meshTests/well1.vtk | 11 + .../unitTests/meshTests/well2.vtk | 11 + .../testReservoirSinglePhaseMSWells.cpp | 257 ++++++--- src/docs/sphinx/CompleteXMLSchema.rst | 14 + 32 files changed, 1565 insertions(+), 926 deletions(-) create mode 100644 src/coreComponents/mesh/generators/VTKWellGenerator.cpp create mode 100644 src/coreComponents/mesh/generators/VTKWellGenerator.hpp create mode 100644 src/coreComponents/mesh/generators/WellGeneratorABC.hpp create mode 100644 src/coreComponents/schema/docs/AcousticElasticSEM.rst create mode 100644 src/coreComponents/schema/docs/AcousticElasticSEM_other.rst create mode 100644 src/coreComponents/schema/docs/PorousDamageElasticIsotropic.rst create mode 100644 src/coreComponents/schema/docs/PorousDamageElasticIsotropic_other.rst create mode 100644 src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic.rst create mode 100644 src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic_other.rst create mode 100644 src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic.rst create mode 100644 src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic_other.rst create mode 100644 src/coreComponents/schema/docs/VTKWell.rst create mode 100644 src/coreComponents/schema/docs/VTKWell_other.rst create mode 100644 src/coreComponents/unitTests/meshTests/well1.vtk create mode 100644 src/coreComponents/unitTests/meshTests/well2.vtk diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 52b2fd69b60..6065ef14dae 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -59,6 +59,7 @@ set( mesh_headers generators/ParticleMeshGenerator.hpp generators/PartitionDescriptor.hpp generators/PrismUtilities.hpp + generators/WellGeneratorABC.hpp generators/WellGeneratorBase.hpp mpiCommunications/CommID.hpp mpiCommunications/CommunicationTools.hpp @@ -150,19 +151,23 @@ set( mesh_sources set( dependencyList ${parallelDeps} schema dataRepository constitutive finiteElement parmetis metis ) if( ENABLE_VTK ) - message(STATUS "Adding VTKMeshGenerator sources and headers") + message(STATUS "Adding VTK readers") set( mesh_headers ${mesh_headers} generators/CollocatedNodes.hpp generators/VTKFaceBlockUtilities.hpp generators/VTKMeshGenerator.hpp generators/VTKMeshGeneratorTools.hpp - generators/VTKUtilities.hpp ) + generators/VTKWellGenerator.hpp + generators/VTKUtilities.hpp + ) set( mesh_sources ${mesh_sources} generators/CollocatedNodes.cpp generators/VTKFaceBlockUtilities.cpp generators/VTKMeshGenerator.cpp generators/VTKMeshGeneratorTools.cpp - generators/VTKUtilities.cpp ) + generators/VTKWellGenerator.cpp + generators/VTKUtilities.cpp + ) list( APPEND dependencyList VTK::IOLegacy VTK::FiltersParallelDIY2 ) if( ENABLE_MPI ) list( APPEND dependencyList VTK::IOParallelXML VTK::ParallelMPI ) diff --git a/src/coreComponents/mesh/generators/InternalWellGenerator.cpp b/src/coreComponents/mesh/generators/InternalWellGenerator.cpp index 22060669605..5e5d9211a58 100644 --- a/src/coreComponents/mesh/generators/InternalWellGenerator.cpp +++ b/src/coreComponents/mesh/generators/InternalWellGenerator.cpp @@ -19,29 +19,12 @@ #include "InternalWellGenerator.hpp" -#include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "mesh/Perforation.hpp" -#include "mesh/generators/LineBlockABC.hpp" -#include "LvArray/src/genericTensorOps.hpp" -#include "physicsSolvers/fluidFlow/wells/WellControls.hpp" - namespace geos { using namespace dataRepository; InternalWellGenerator::InternalWellGenerator( string const & name, Group * const parent ): - WellGeneratorBase( name, parent ), - m_numElemsPerSegment( 0 ), - m_minSegmentLength( 1e-2 ), - m_minElemLength( 1e-3 ), - m_radius( 0 ), - m_wellRegionName( "" ), - m_wellControlsName( "" ), - m_numElems( 0 ), - m_numNodesPerElem( 2 ), - m_numNodes( 0 ), - m_nDims( 3 ), - m_polylineHeadNodeId( -1 ) + WellGeneratorBase( name, parent ) { enableLogLevelInput(); @@ -54,40 +37,6 @@ InternalWellGenerator::InternalWellGenerator( string const & name, Group * const setInputFlag( InputFlags::REQUIRED ). setSizedFromParent( 0 ). setDescription( "Connectivity of the polyline segments" ); - - registerWrapper( viewKeyStruct::radiusString(), &m_radius ). - setInputFlag( InputFlags::REQUIRED ). - setSizedFromParent( 0 ). - setDescription( "Radius of the well [m]" ); - - registerWrapper( viewKeyStruct::numElementsPerSegmentString(), &m_numElemsPerSegment ). - setInputFlag( InputFlags::REQUIRED ). - setSizedFromParent( 0 ). - setDescription( "Number of well elements per polyline segment" ); - - registerWrapper( viewKeyStruct::minSegmentLengthString(), &m_minSegmentLength ). - setApplyDefaultValue( 1e-2 ). - setInputFlag( InputFlags::OPTIONAL ). - setRestartFlags( RestartFlags::NO_WRITE ). - setDescription( "Minimum length of a well segment [m]" ); - - registerWrapper( viewKeyStruct::minElementLengthString(), &m_minElemLength ). - setApplyDefaultValue( 1e-3 ). - setInputFlag( InputFlags::OPTIONAL ). - setRestartFlags( RestartFlags::NO_WRITE ). - setDescription( "Minimum length of a well element, computed as (segment length / number of elements per segment ) [m]" ); - - registerWrapper( viewKeyStruct::wellRegionNameString(), &m_wellRegionName ). - setRTTypeName( rtTypes::CustomTypes::groupNameRef ). - setInputFlag( InputFlags::REQUIRED ). - setSizedFromParent( 0 ). - setDescription( "Name of the well element region" ); - - registerWrapper( viewKeyStruct::wellControlsNameString(), &m_wellControlsName ). - setRTTypeName( rtTypes::CustomTypes::groupNameRef ). - setInputFlag( InputFlags::REQUIRED ). - setSizedFromParent( 0 ). - setDescription( "Name of the set of constraints associated with this well" ); } void InternalWellGenerator::postProcessInput() @@ -107,492 +56,11 @@ void InternalWellGenerator::postProcessInput() " and " << getWrapperDataContext( viewKeyStruct::polylineSegmentConnString() ), InputError ); - GEOS_THROW_IF( m_radius <= 0, - "InternalWell " << getWrapperDataContext( viewKeyStruct::radiusString() ) << - ": Radius value must be greater that 0.", - InputError ); - - GEOS_THROW_IF( m_wellRegionName.empty(), - "InternalWell " << getWrapperDataContext( viewKeyStruct::wellRegionNameString() ) << - ": Empty well region name.", - InputError ); - - GEOS_THROW_IF( m_wellControlsName.empty(), - "InternalWell " << getWrapperDataContext( viewKeyStruct::wellControlsNameString() ) << - ": Empty well constraint name.", - InputError ); - // TODO: add more checks here // TODO: check that the connectivity of the well is valid // TODO: check that with no branching we can go from top to bottom and touch all the elements } -void InternalWellGenerator::generateWellGeometry( ) -{ - // count the number of well elements to create - m_numElems = m_numElemsPerSegment * m_segmentToPolyNodeMap.size( 0 ); - m_numNodes = m_numElems + 1; - - // resize the well element, node, and perforation arrays - m_elemCenterCoords.resize( m_numElems, 3 ); - m_nextElemId.resize( m_numElems ); - m_prevElemId.resize( m_numElems ); - m_elemToNodesMap.resizeDimension< 0 >( m_numElems ); - m_elemToNodesMap.resizeDimension< 1 >( m_numNodesPerElem ); - m_elemVolume.resize( m_numElems ); - - m_nodeDistFromHead.resize( m_numNodes ); - m_nodeCoords.resize( m_numNodes, 3 ); - - m_perfCoords.resize( m_numPerforations, 3 ); - m_perfDistFromHead.resize( m_numPerforations ); - m_perfTransmissibility.resize( m_numPerforations ); - m_perfSkinFactor.resize( m_numPerforations ); - m_perfElemId.resize( m_numPerforations ); - - // construct a reverse map from the polyline nodes to the segments - constructPolylineNodeToSegmentMap(); - - // detect the head polyline node based on depth - findPolylineHeadNodeIndex(); - - // compute the location and distance from top of the well elements - discretizePolyline(); - - // map the perforations to the well elements - connectPerforationsToWellElements(); - - // make sure that the perforation locations are valid - checkPerforationLocationsValidity(); - - if( getLogLevel() >= 1 ) - { - debugWellGeometry(); - } - -} - -void InternalWellGenerator::constructPolylineNodeToSegmentMap() -{ - m_polyNodeToSegmentMap.resize( m_polyNodeCoords.size( 0 ) ); - - bool foundSmallElem = false; - - // loop over the segments - for( globalIndex iseg = 0; iseg < m_segmentToPolyNodeMap.size( 0 ); ++iseg ) - { - globalIndex const ipolyNode_a = m_segmentToPolyNodeMap[iseg][0]; - globalIndex const ipolyNode_b = m_segmentToPolyNodeMap[iseg][1]; - - real64 vSeg[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_polyNodeCoords[ipolyNode_a] ); - LvArray::tensorOps::subtract< 3 >( vSeg, m_polyNodeCoords[ipolyNode_b] ); - - real64 const segmentLength = LvArray::tensorOps::l2Norm< 3 >( vSeg ); - - // various checks and warnings on the segment and element length - GEOS_THROW_IF( segmentLength < m_minSegmentLength, - getDataContext() << ": Error in the topology of the well. We detected a" << - " polyline segment measuring less than " << m_minSegmentLength << "m. \n" << - "You can change the minimum segment length using the field " << viewKeyStruct::minSegmentLengthString(), - InputError ); - - GEOS_THROW_IF( m_polyNodeCoords[ipolyNode_a][2] < m_polyNodeCoords[ipolyNode_b][2], - getDataContext() << ": Error in the topology of the well. In the polyline"<< - ", each segment must be going down. \n" << - "This is not the case between polyline nodes " << m_polyNodeCoords[ipolyNode_a] << " and " << m_polyNodeCoords[ipolyNode_b], - InputError ); - - if( segmentLength / m_numElemsPerSegment < m_minElemLength ) - { - foundSmallElem = true; - } - - // map the polyline node ids to the polyline segment ids - m_polyNodeToSegmentMap[ipolyNode_a].insert( iseg ); - m_polyNodeToSegmentMap[ipolyNode_b].insert( iseg ); - } - - if( foundSmallElem ) - { - GEOS_LOG_RANK_0( "\nWarning for well " << getDataContext() << ": the chosen number of well" << - " elements per polyline segment (" << m_numElemsPerSegment << ") leads to well" << - " elements measuring less than " << m_minElemLength << "m in the topology.\n" << - "The simulation can proceed like that, but small well elements may cause" << - " numerical issues, so it is something to keep an eye on.\n" << - "You can get rid of this message by changing the field " << viewKeyStruct::minElementLengthString() ); - } -} - -void InternalWellGenerator::findPolylineHeadNodeIndex() -{ - // we assume that the first *polyline segment* in the XML input is the well head segment - globalIndex const wellHeadSegId = 0; - - // get the corresponding node indices - globalIndex const ipolyNode_a = m_segmentToPolyNodeMap[wellHeadSegId][0]; - globalIndex const ipolyNode_b = m_segmentToPolyNodeMap[wellHeadSegId][1]; - - // we determine which node is the head node based on depth - // therefore here we throw an error if the well head segment is horizontal - GEOS_THROW_IF( !(m_polyNodeCoords[ipolyNode_a][2] < m_polyNodeCoords[ipolyNode_b][2]) - && !(m_polyNodeCoords[ipolyNode_a][2] > m_polyNodeCoords[ipolyNode_b][2]), - getDataContext() << ": The head polyline segment cannot be horizontal in the well" << - " since we use depth to determine which of its nodes is to head node of the well.\n" << - "If you are trying to set up a horizontal well, please simply add a non-horizontal" << - " segment at the top of the well, and this error will go away.", - InputError ); - - // detect the top node, assuming z oriented upwards - m_polylineHeadNodeId = - ( m_polyNodeCoords[ipolyNode_a][2] > m_polyNodeCoords[ipolyNode_b][2] ) - ? ipolyNode_a - : ipolyNode_b; - - real64 const headZcoord = m_polyNodeCoords[m_polylineHeadNodeId][2]; - for( globalIndex inode = 0; inode < m_polyNodeCoords.size( 0 ); ++inode ) - { - if( inode == m_polylineHeadNodeId ) - { - continue; - } - real64 const currentZcoord = m_polyNodeCoords[inode][2]; - - GEOS_THROW_IF( !(currentZcoord < headZcoord), - getDataContext() << ": Error in the topology since we found a well node that" << - " is above the head node", - InputError ); - } -} - -void InternalWellGenerator::discretizePolyline() -{ - // initialize well elements and node ids - globalIndex ipolyNodeTop = m_polylineHeadNodeId; - globalIndex iwelemCurrent = 0; - globalIndex isegCurrent = 0; - - // set the location of the first well node and distance from well head - LvArray::tensorOps::copy< 3 >( m_nodeCoords[0], m_polyNodeCoords[ipolyNodeTop] ); - m_nodeDistFromHead[0] = 0.0; - - // note: this part of the code does not support well branching - // TODO: check that there is only one branch - // TODO: read wells with branching (already supported elsewhere in the code) - - // go through the well from top to bottom - for( globalIndex is = 0; is < m_segmentToPolyNodeMap.size( 0 ); ++is ) - { - - GEOS_THROW_IF( isegCurrent == -1, - getWrapperDataContext( viewKeyStruct::polylineSegmentConnString() ) << - ": Invalid map.", - InputError ); - - globalIndex const ipolyNodeBottom = ( ipolyNodeTop == m_segmentToPolyNodeMap[isegCurrent][0] ) - ? m_segmentToPolyNodeMap[isegCurrent][1] - : m_segmentToPolyNodeMap[isegCurrent][0]; - - real64 vPoly[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_polyNodeCoords[ipolyNodeBottom] ); - LvArray::tensorOps::subtract< 3 >( vPoly, m_polyNodeCoords[ipolyNodeTop] ); - - // add the well elements and well nodes corresponding to this polyline segment - for( localIndex iw = 0; iw < m_numElemsPerSegment; ++iw ) - { - - // 1) set the element location, connectivity - real64 const scaleCenter = (iw + 0.5) / static_cast< real64 >(m_numElemsPerSegment); - LvArray::tensorOps::copy< 3 >( m_elemCenterCoords[iwelemCurrent], vPoly ); - LvArray::tensorOps::scale< 3 >( m_elemCenterCoords[iwelemCurrent], scaleCenter ); - LvArray::tensorOps::add< 3 >( m_elemCenterCoords[iwelemCurrent], m_polyNodeCoords[ipolyNodeTop] ); - - GEOS_THROW_IF( iwelemCurrent >= m_numElems, - getDataContext() << ": Invalid well topology", - InputError ); - - globalIndex const iwelemTop = iwelemCurrent - 1; - globalIndex const iwelemBottom = iwelemCurrent + 1; - m_nextElemId[iwelemCurrent] = iwelemTop; - m_prevElemId[iwelemCurrent].resize( 1 ); - m_prevElemId[iwelemCurrent][0] = (iwelemBottom < m_numElems) - ? iwelemBottom - : -1; - - // 2) set the node properties - globalIndex const iwellNodeTop = iwelemCurrent; - globalIndex const iwellNodeBottom = iwelemCurrent+1; - m_elemToNodesMap[iwelemCurrent][LineBlockABC::NodeLocation::TOP] = iwellNodeTop; - m_elemToNodesMap[iwelemCurrent][LineBlockABC::NodeLocation::BOTTOM] = iwellNodeBottom; - - real64 const scaleBottom = (iw + 1.0) / m_numElemsPerSegment; - LvArray::tensorOps::copy< 3 >( m_nodeCoords[iwellNodeBottom], vPoly ); - LvArray::tensorOps::scale< 3 >( m_nodeCoords[iwellNodeBottom], scaleBottom ); - LvArray::tensorOps::add< 3 >( m_nodeCoords[iwellNodeBottom], m_polyNodeCoords[ipolyNodeTop] ); - - // 3) increment the distance from the well head to the bottom of current element - real64 vWellElem[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_nodeCoords[iwellNodeBottom] ); - LvArray::tensorOps::subtract< 3 >( vWellElem, m_nodeCoords[iwellNodeTop] ); - m_nodeDistFromHead[iwellNodeBottom] = LvArray::tensorOps::l2Norm< 3 >( vWellElem ); - m_nodeDistFromHead[iwellNodeBottom] += m_nodeDistFromHead[iwellNodeTop]; - - // 4) set element volume - m_elemVolume[iwelemCurrent] = LvArray::tensorOps::l2Norm< 3 >( vWellElem ) * M_PI * m_radius * m_radius; - - // 4) increment the element counter - ++iwelemCurrent; - } - - // then consider the next polyline segment - ipolyNodeTop = ipolyNodeBottom; - isegCurrent = getNextSegmentIndex( isegCurrent, ipolyNodeTop ); - } - - // set the previous index for the bottom segment - m_prevElemId[iwelemCurrent-1].resize( 1 ); - m_prevElemId[iwelemCurrent-1][0] = -1; - -} - - -void InternalWellGenerator::connectPerforationsToWellElements() -{ - - // assign a well element to each perforation - for( globalIndex iperf = 0; iperf < m_numPerforations; ++iperf ) - { - - // get the perforation and its properties - Perforation const & perf = this->getGroup< Perforation >( m_perforationList[iperf] ); - m_perfDistFromHead[iperf] = perf.getDistanceFromWellHead(); - m_perfTransmissibility[iperf] = perf.getWellTransmissibility(); - m_perfSkinFactor[iperf] = perf.getWellSkinFactor(); - - // search in all the elements of this well between head and bottom - globalIndex iwelemTop = 0; - globalIndex iwelemBottom = m_numElems - 1; - - // check the validity of the perforation before starting - real64 const wellLength = m_nodeDistFromHead[m_elemToNodesMap[iwelemBottom][LineBlockABC::NodeLocation::BOTTOM]]; - - GEOS_THROW_IF( m_perfDistFromHead[iperf] > wellLength, - perf.getWrapperDataContext( Perforation::viewKeyStruct::distanceFromHeadString() ) << - ": Distance from well perforation to head (" << - Perforation::viewKeyStruct::distanceFromHeadString() << " = " << m_perfDistFromHead[iperf] << - ") is larger than well polyline length (" << wellLength << - ")\n \n You should check the following values:" << - "\n 1 - " << perf.getWrapperDataContext( Perforation::viewKeyStruct::distanceFromHeadString() ) << - "\n 2 - " << getWrapperDataContext( viewKeyStruct::polylineNodeCoordsString() ) << ", Z values", - InputError ); - - // start binary search - const globalIndex maxNumSteps = m_numElems + 1; - globalIndex currentNumSteps = 0; - while( iwelemTop < iwelemBottom ) - { - globalIndex iwelemMid = - static_cast< globalIndex >(floor( static_cast< real64 >(iwelemTop + iwelemBottom) / 2.0 )); - real64 const headToBottomDist = - m_nodeDistFromHead[m_elemToNodesMap[iwelemMid][LineBlockABC::NodeLocation::BOTTOM]]; - - if( headToBottomDist < m_perfDistFromHead[iperf] ) - { - iwelemTop = iwelemMid + 1; - } - else - { - iwelemBottom = iwelemMid; - } - - GEOS_THROW_IF( currentNumSteps > maxNumSteps, - perf.getDataContext() << ": Perforation cannot be mapped to a well element.", - InputError ); - - currentNumSteps++; - } - - // set the index of the matched element - globalIndex iwelemMatched = iwelemTop; - GEOS_THROW_IF( iwelemMatched >= m_numElems, - getDataContext() << ": Invalid well topology.", - InputError ); - - m_perfElemId[iperf] = iwelemMatched; - - // compute the physical location of the perforation - globalIndex const inodeTop = m_elemToNodesMap[iwelemMatched][LineBlockABC::NodeLocation::TOP]; - globalIndex const inodeBottom = m_elemToNodesMap[iwelemMatched][LineBlockABC::NodeLocation::BOTTOM]; - real64 const elemLength = m_nodeDistFromHead[inodeBottom] - m_nodeDistFromHead[inodeTop]; - real64 const topToPerfDist = m_perfDistFromHead[iperf] - m_nodeDistFromHead[inodeTop]; - - GEOS_THROW_IF( (elemLength <= 0) || (topToPerfDist < 0), - getDataContext() << ": Invalid well topology.", - InputError ); - - LvArray::tensorOps::copy< 3 >( m_perfCoords[iperf], m_nodeCoords[inodeBottom] ); - LvArray::tensorOps::subtract< 3 >( m_perfCoords[iperf], m_nodeCoords[inodeTop] ); - LvArray::tensorOps::scale< 3 >( m_perfCoords[iperf], topToPerfDist / elemLength ); - LvArray::tensorOps::add< 3 >( m_perfCoords[iperf], m_nodeCoords[inodeTop] ); - - } -} - - -globalIndex InternalWellGenerator::getNextSegmentIndex( globalIndex topSegId, - globalIndex currentPolyNodeId ) const -{ - globalIndex nextSegId = -1; - - // get the index of the two segments sharing this node - for( auto is : m_polyNodeToSegmentMap[currentPolyNodeId] ) - { - // if this is not equal to the current segId, we found the next segId - if( is != topSegId ) - { - nextSegId = is; - break; - } - } - - return nextSegId; -} - -void InternalWellGenerator::checkPerforationLocationsValidity() -{ - array1d< array1d< localIndex > > elemToPerfMap; - elemToPerfMap.resize( m_numElems ); - - for( globalIndex iperf = 0; iperf < m_numPerforations; ++iperf ) - { - elemToPerfMap[m_perfElemId[iperf]].emplace_back( iperf ); - } - - // merge perforations to make sure that no well element is shared between two MPI domains - // TODO: instead of merging perforations, split the well elements and do not change the physical location of the - // perforation - int const mpiSize = MpiWrapper::commSize( MPI_COMM_GEOSX ); - if( mpiSize > 1 ) - { - mergePerforations( elemToPerfMap ); - } - - for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) - { - // check that there is always a perforation in the last well element (otherwise, the problem is not well posed) - for( localIndex iwelemPrev = 0; iwelemPrev < m_prevElemId[iwelem].size(); ++iwelemPrev ) - { - GEOS_THROW_IF( m_prevElemId[iwelem][iwelemPrev] == -1 && elemToPerfMap[iwelem].size() == 0, - getDataContext() << "The bottom element of the well does not have a perforation. " << - "This is needed to have a well-posed problem. \n\n" << - "Here are the two possible ways to solve this problem: \n\n" << - "1) Adding a perforation located close to the bottom of the well. " << - "To do that, compute the total length of the well polyline (by summing the length of the well segments defined by the keywords \"polylineNodeCoords\" and \"polylineSegmentConn\") " << - "and place a perforation whose \"distanceFromHead\" is slightly smaller than this total length. \n \n" << - "2) Shorten the well polyline. " << - "To do that, reduce the length of the well polyline by shortening the segments defined by the keywords \"polylineNodeCoords\" and \"polylineSegmentConn\", or by removing a segment.", - InputError ); - } - } -} - -void InternalWellGenerator::mergePerforations( array1d< array1d< localIndex > > const & elemToPerfMap ) -{ - - for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) - { - // collect the indices of the elems with more that one perforation - if( elemToPerfMap[iwelem].size() > 1 ) - { - // find the perforation with the largest Peaceman index and keep its location - globalIndex iperfMaxTransmissibility = elemToPerfMap[iwelem][0]; - real64 maxTransmissibility = m_perfTransmissibility[iperfMaxTransmissibility]; - for( localIndex ip = 1; ip < elemToPerfMap[iwelem].size(); ++ip ) - { - if( m_perfTransmissibility[elemToPerfMap[iwelem][ip]] > maxTransmissibility ) - { - iperfMaxTransmissibility = elemToPerfMap[iwelem][ip]; - maxTransmissibility = m_perfTransmissibility[iperfMaxTransmissibility]; - } - } - - // assign the coordinates of the perf with the largest trans to the other perfs on this elem - for( localIndex ip = 0; ip < elemToPerfMap[iwelem].size(); ++ip ) - { - if( elemToPerfMap[iwelem][ip] == iperfMaxTransmissibility ) - { - continue; - } - - GEOS_LOG_RANK_0( "\n \nThe GEOSX wells currently have the following limitation in parallel: \n" << - "We cannot allow an element of the well mesh to have two or more perforations associated with it. \n" << - "So, in the present simulation, perforation #" << elemToPerfMap[iwelem][ip] << - " of well " << getDataContext() << - " is moved from " << m_perfCoords[elemToPerfMap[iwelem][ip]] << - " to " << m_perfCoords[iperfMaxTransmissibility] << - " to make sure that no element of the well mesh has two perforations associated with it. \n" << - "To circumvent this issue, please increase the value of \"numElementsPerSegment\" that controls the number of (uniformly distributed) well mesh elements per segment of the well polyline. \n" << - "Our recommendation is to choose \"numElementsPerSegment\" such that each well mesh element has at most one perforation. \n\n" ); - LvArray::tensorOps::copy< 3 >( m_perfCoords[elemToPerfMap[iwelem][ip]], m_perfCoords[iperfMaxTransmissibility] ); - } - } - } -} - -void InternalWellGenerator::debugWellGeometry() const -{ - if( MpiWrapper::commRank( MPI_COMM_GEOSX ) != 0 ) - { - return; - } - - std::cout << std::endl; - std::cout << "++++++++++++++++++++++++++" << std::endl; - std::cout << "InternalWellGenerator = " << getName() << std::endl; - std::cout << "MPI rank = " << MpiWrapper::commRank( MPI_COMM_GEOSX ) << std::endl << std::endl; - std::cout << "Number of well elements = " << m_numElems << std::endl; - - for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) - { - std::cout << "Well element #" << iwelem << std::endl; - std::cout << "Coordinates of the element center: " << m_elemCenterCoords[iwelem] << std::endl; - if( m_nextElemId[iwelem] < 0 ) - { - std::cout << "No next well element" << std::endl; - } - else - { - std::cout << "Next well element # = " << m_nextElemId[iwelem] << std::endl; - } - if( m_prevElemId[iwelem][0] < 0 ) - { - std::cout << "No previous well element" << std::endl; - } - else - { - std::cout << "Previous well element #" << m_prevElemId[iwelem][0] << std::endl; - } - for( globalIndex inode = 0; inode < m_numNodesPerElem; ++inode ) - { - if( inode == 0 ) - { - std::cout << "First well node: #" << m_elemToNodesMap[iwelem][inode] << std::endl; - } - else - { - std::cout << "Second well node: #" << m_elemToNodesMap[iwelem][inode] << std::endl; - } - } - } - - std::cout << std::endl << "Number of perforations = " << m_numPerforations << std::endl; - - for( globalIndex iperf = 0; iperf < m_numPerforations; ++iperf ) - { - std::cout << "Perforation #" << iperf << std::endl; - std::cout << "Coordinates of the perforation: " << m_perfCoords[iperf] << std::endl; - std::cout << "Is connected to well element #" << m_perfElemId[iperf] << std::endl; - } - std::cout << std::endl; - -} REGISTER_CATALOG_ENTRY( WellGeneratorBase, InternalWellGenerator, string const &, Group * const ) } diff --git a/src/coreComponents/mesh/generators/InternalWellGenerator.hpp b/src/coreComponents/mesh/generators/InternalWellGenerator.hpp index 1093124dc4e..16cdeee3553 100644 --- a/src/coreComponents/mesh/generators/InternalWellGenerator.hpp +++ b/src/coreComponents/mesh/generators/InternalWellGenerator.hpp @@ -59,301 +59,14 @@ class InternalWellGenerator : public WellGeneratorBase */ static string catalogName() { return "InternalWell"; } - - /** - * @brief Main function of the class that generates the well geometry - */ - void generateWellGeometry( ) override; - - - ///@} - - /** - * @name Getters / Setters - */ - ///@{ - - // getters for element data - - /** - * @brief Get the global number of well elements. - * @return the global number of elements - */ - globalIndex numElements() const override { return m_numElems; } - - /** - * @brief Getter to the Segment to PolyNode mapping - * @return The Segment to PolyNode mapping as a 2D array - */ - const array2d< globalIndex > & getSegmentToPolyNodeMap() const override { return m_segmentToPolyNodeMap; }; - - /** - * @brief Get the number of nodes per well element - * @return the number of nodes per well element - */ - globalIndex numNodesPerElement() const override { return m_numNodesPerElem; } - - /** - * @brief Get the Coordinates of the polyline nodes - * @return the Coordinates of the polyline nodes - */ - const array2d< real64 > & getPolyNodeCoord() const override { return m_polyNodeCoords; } - - /** - * @return The minimum segment length - */ - real64 getMinSegmentLength() const override { return m_minSegmentLength; } - - /** - * @return The minimum element length - */ - real64 getMinElemLength() const override { return m_minElemLength; } - - /** - * @brief Get the physical location of the centers of well elements. - * @return list of center locations of the well elements - */ - arrayView2d< real64 const > getElemCoords() const override { return m_elemCenterCoords; } - - /** - * @brief Get the global indices mapping an element to the next. - * @return list providing the global index of the next element for each element - */ - arrayView1d< globalIndex const > getNextElemIndex() const override { return m_nextElemId; } - - /** - * @brief Get the global indices mapping an element to the previous ones. - * @return list providing the global indices of the previous elements for each element - */ - arrayView1d< arrayView1d< globalIndex const > const > getPrevElemIndices() const override { return m_prevElemId.toNestedViewConst(); } - - /** - * @brief Get the global indices of the well nodes nodes connected to each element. - * @return list providing the global index of the well nodes for each well element - */ - arrayView2d< globalIndex const > getElemToNodesMap() const override { return m_elemToNodesMap; } - - /** - * @brief Get the volume of the well elements. - * @return list of volumes of the well elements - */ - arrayView1d< real64 const > getElemVolume() const override { return m_elemVolume; } - - /** - * @brief Get the radius in the well. - * @return the radius in the well - */ - real64 getElementRadius() const override { return m_radius; } - - // getters for node data - - /** - * @brief Get the global number of well nodes. - * @return the global number of nodes - */ - globalIndex numNodes() const override { return m_numNodes; } - - /** - * @brief Get the physical location of the centers of well elements. - * @return list of center locations of the well elements - */ - arrayView2d< real64 const > getNodeCoords() const override { return m_nodeCoords; } - - // getters for perforation data - - /** - * @brief Get the locations of the perforations. - * @return list of locations of all the perforations on the well - */ - arrayView2d< real64 const > getPerfCoords() const override { return m_perfCoords; } - - /** - * @brief Get the well transmissibility at the perforations. - * @return list of well transmissibility at all the perforations on the well - */ - arrayView1d< real64 const > getPerfTransmissibility() const override { return m_perfTransmissibility; } - - /** - * @brief Get the well skin factor at the perforations. - * @return list of well skin factor at all the perforations on the well - */ - arrayView1d< real64 const > getPerfSkinFactor() const override { return m_perfSkinFactor; } - - /** - * @brief Get the global indices of the well elements connected to each perforation. - * @return list providing the global index of the connected well element for each perforation - */ - arrayView1d< globalIndex const > getPerfElemIndex() const override { return m_perfElemId; } - - /** - * @returns The number of physical dimensions - */ - int getPhysicalDimensionsNumber() const override { return m_nDims; } - ///@} - const string getWellRegionName() const override { return m_wellRegionName; } - const string getWellControlsName() const override { return m_wellControlsName; } - - - protected: - /** * @brief This function provides capability to post process input values prior to * any other initialization operations. */ - void postProcessInput() override final; - -private: - - /** - * @name Helper functions to construct the geometry of the well - */ - ///@{ - - /** - * @brief Map each polyline node to the polyline segment(s) it is connected to. - */ - void constructPolylineNodeToSegmentMap(); - - /** - * @brief Find the head node of the well (i.e., top node of the polyline). - */ - void findPolylineHeadNodeIndex(); - - /** - * @brief Discretize the polyline by placing well elements. - */ - void discretizePolyline(); - - /** - * @brief Map each perforation to a well element. - */ - void connectPerforationsToWellElements(); - - /** - * @brief Make sure that the perforation locations are valid: - * - for partitioning purposes - * - to have a well-posed problem - */ - void checkPerforationLocationsValidity(); - - /** - * @brief Merge perforations on the elements with multiple perforations. - */ - void mergePerforations( array1d< array1d< localIndex > > const & elemToPerfMap ); - - /** - * @brief At a given node, find the next segment going in the direction of the bottom of the well. - * @param[in] topSegId index of the top segment - * @param[in] currentNodeId index of the current node - */ - globalIndex getNextSegmentIndex( globalIndex topSegId, - globalIndex currentNodeId ) const; - - ///@} - - /// @cond DO_NOT_DOCUMENT - void debugWellGeometry() const; - /// @endcond - - // XML Input - - /// Connectivity between the polyline nodes - array2d< globalIndex > m_segmentToPolyNodeMap; - - /// Number of well elements per polyline interval - int m_numElemsPerSegment; - - /// Min segment length - real64 m_minSegmentLength; - - /// Min well element length - real64 m_minElemLength; - - /// Radius area of the well (assumed to be valid for the entire well) - real64 m_radius; - - /// Name of the corresponding well region - string m_wellRegionName; - - /// Name of the constraints associated with this well - string m_wellControlsName; - - - - // Geometry of the well (later passed to the WellElementSubRegion) - - // well element data - - /// Global number of well elements - globalIndex m_numElems; - - /// Physical location of the center of the well element - array2d< real64 > m_elemCenterCoords; - - /// Global index of the next well element - array1d< globalIndex > m_nextElemId; - - /// Global indices of the prev well elements (maybe need multiple prevs for branching) - array1d< array1d< globalIndex > > m_prevElemId; - - /// Connectivity between elements and nodes - array2d< globalIndex > m_elemToNodesMap; - - /// Volume of well elements - array1d< real64 > m_elemVolume; - - - // well node data - - /// Number of nodes per well element - globalIndex const m_numNodesPerElem; - - /// Global number of well nodes - globalIndex m_numNodes; - - /// Physical location of the nodes - array2d< real64 > m_nodeCoords; - - // perforation data - - /// Absolute physical location of the perforation - array2d< real64 > m_perfCoords; - - /// Well Peaceman index at the perforation - array1d< real64 > m_perfTransmissibility; - - /// Skin Factor at the perforation - array1d< real64 > m_perfSkinFactor; - - /// Global index of the well element - array1d< globalIndex > m_perfElemId; - - - - // Auxiliary data - - // Number of physical dimensions - const int m_nDims; - - /// Coordinates of the polyline nodes - array2d< real64 > m_polyNodeCoords; - - /// Map from the polyline nodes to the polyline nodes - array1d< SortedArray< globalIndex > > m_polyNodeToSegmentMap; - - /// Index of the node at the well head - globalIndex m_polylineHeadNodeId; - - /// Physical location of the polyline node wrt to well head - array1d< real64 > m_nodeDistFromHead; - - // Perforation data - - /// Physical location of the perforation wrt to well head - array1d< real64 > m_perfDistFromHead; + void postProcessInput() override; }; } diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 25be52ebe7a..47bce1e7aa0 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -37,6 +37,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -48,6 +51,8 @@ #include #include #include +#include +#include #include #include #include @@ -81,10 +86,12 @@ enum class VTKMeshExtension : integer vtr, ///< XML serial vtkRectilinearGrid (structured) vts, ///< XML serial vtkStructuredGrid (structured) vti, ///< XML serial vtkImageData (structured) + vtp, ///< XML serial vtkPolyData pvtu, ///< XML parallel vtkUnstructuredGrid (unstructured) pvtr, ///< XML parallel vtkRectilinearGrid (structured) pvts, ///< XML parallel vtkStructuredGrid (structured) pvti, ///< XML parallel vtkImageData (structured) + pvtp, ///< XML parallel vtkPolyData }; /// Strings for VTKMeshGenerator::VTKMeshExtension enumeration @@ -95,10 +102,12 @@ ENUM_STRINGS( VTKMeshExtension, "vtr", "vts", "vti", + "vtp", "pvtu", "pvtr", "pvts", - "pvti" ); + "pvti", + "pvtp" ); /** * @brief Supported VTK legacy dataset types @@ -109,6 +118,7 @@ enum class VTKLegacyDatasetType : integer structuredGrid, ///< Structured grid (structured) unstructuredGrid, ///< Unstructured grid (unstructured) rectilinearGrid, ///< Rectilinear grid (structured) + polyData, ///< PolyData }; /// Strings for VTKMeshGenerator::VTKLegacyDatasetType enumeration @@ -116,7 +126,8 @@ ENUM_STRINGS( VTKLegacyDatasetType, "structuredPoints", "structuredGrid", "unstructuredGrid", - "rectilinearGrid" ); + "rectilinearGrid", + "polyData" ); /** * @brief Gathers all the data from all ranks, merge them, sort them, and remove duplicates. @@ -423,6 +434,10 @@ VTKLegacyDatasetType getVTKLegacyDatasetType( vtkSmartPointer< vtkDataSetReader { return VTKLegacyDatasetType::rectilinearGrid; } + else if( vtkGridReader->IsFilePolyData()) + { + return VTKLegacyDatasetType::polyData; + } else { GEOS_ERROR( "Unsupported legacy VTK dataset format.\nLegacy supported formats are: " << @@ -515,6 +530,7 @@ loadMesh( Path const & filePath, case VTKLegacyDatasetType::structuredGrid: return serialRead( vtkSmartPointer< vtkStructuredGridReader >::New() ); case VTKLegacyDatasetType::unstructuredGrid: return serialRead( vtkSmartPointer< vtkUnstructuredGridReader >::New() ); case VTKLegacyDatasetType::rectilinearGrid: return serialRead( vtkSmartPointer< vtkRectilinearGridReader >::New() ); + case VTKLegacyDatasetType::polyData: return serialRead( vtkSmartPointer< vtkPolyDataReader >::New() ); } break; } @@ -522,6 +538,7 @@ loadMesh( Path const & filePath, case VTKMeshExtension::vtr: return serialRead( vtkSmartPointer< vtkXMLRectilinearGridReader >::New() ); case VTKMeshExtension::vts: return serialRead( vtkSmartPointer< vtkXMLStructuredGridReader >::New() ); case VTKMeshExtension::vti: return serialRead( vtkSmartPointer< vtkXMLImageDataReader >::New() ); + case VTKMeshExtension::vtp: return serialRead( vtkSmartPointer< vtkXMLPolyDataReader >::New() ); case VTKMeshExtension::pvtu: { return parallelRead( vtkSmartPointer< vtkXMLPUnstructuredGridReader >::New() ); @@ -532,6 +549,7 @@ loadMesh( Path const & filePath, case VTKMeshExtension::pvts: return parallelRead( vtkSmartPointer< vtkXMLPStructuredGridReader >::New() ); case VTKMeshExtension::pvtr: return parallelRead( vtkSmartPointer< vtkXMLPRectilinearGridReader >::New() ); case VTKMeshExtension::pvti: return parallelRead( vtkSmartPointer< vtkXMLPImageDataReader >::New() ); + case VTKMeshExtension::pvtp: return parallelRead( vtkSmartPointer< vtkXMLPPolyDataReader >::New() ); default: { GEOS_ERROR( extension << " is not a recognized extension for VTKMesh. Please use ." << EnumStrings< VTKMeshExtension >::concat( ", ." ) ); diff --git a/src/coreComponents/mesh/generators/VTKWellGenerator.cpp b/src/coreComponents/mesh/generators/VTKWellGenerator.cpp new file mode 100644 index 00000000000..3dc7fe343bf --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKWellGenerator.cpp @@ -0,0 +1,105 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/* + * @file VTKWellGenerator.cpp + * + */ + +#include "VTKWellGenerator.hpp" + +#include "mesh/generators/VTKUtilities.hpp" +#include +#include +#include +#include +#include + +namespace geos +{ +using namespace dataRepository; + +VTKWellGenerator::VTKWellGenerator( string const & name, Group * const parent ): + WellGeneratorBase( name, parent ) +{ + registerWrapper( viewKeyStruct::filePathString(), &m_filePath ). + setInputFlag( InputFlags::REQUIRED ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( "Path to the well file" ); + +} + +void VTKWellGenerator::fillPolylineDataStructure( ) +{ + GEOS_MARK_FUNCTION; + + vtkSmartPointer< vtkMultiProcessController > controller = vtk::getController(); + vtkMultiProcessController::SetGlobalController( controller ); + + GEOS_LOG_RANK_0( GEOS_FMT( "{} '{}': reading well from {}", catalogName(), getName(), m_filePath ) ); + { + GEOS_LOG_LEVEL_RANK_0( 2, " reading the dataset..." ); + vtk::AllMeshes allMeshes = vtk::loadAllMeshes( m_filePath, "main", array1d< string >()); + vtkSmartPointer< vtkDataSet > loadedMesh = allMeshes.getMainMesh(); + controller->Broadcast( loadedMesh, 0 ); + + vtkSmartPointer< vtkPolyData > polyData = vtkPolyData::SafeDownCast( loadedMesh ); + + // load points + vtkPoints * points = polyData->GetPoints(); + m_polyNodeCoords.resize( points->GetNumberOfPoints(), 3 ); + globalIndex ipoint = 0; + for( vtkIdType c = 0; c < points->GetNumberOfPoints(); ++c, ++ipoint ) + { + real64 point[3]; + points->GetPoint( c, point ); + LvArray::tensorOps::copy< 3 >( m_polyNodeCoords[ipoint], point ); + } + + GEOS_ERROR_IF( polyData->GetLines()->GetNumberOfCells() == 0, GEOS_FMT( "{}: Error! Your VTK file {} doesn't contain any well", + this->getName(), m_filePath )); + + GEOS_LOG_RANK_0_IF( polyData->GetLines()->GetNumberOfCells() > 1, GEOS_FMT( "{}: Warning! Your VTK file {} contains multiple wells. Only the first one will be read", + this->getName(), m_filePath )); + + // load edges + // polyData->GetLines()->InitTraversal(); + // vtkNew< vtkIdList > idList; + // polyData->GetLines()->GetNextCell( idList ); + + // vtkNew< vtkIdList > idList2; + // polyData->GetLines()->GetCell( 0, idList2 ); + + //newer version of vtk prefer local thread-safe iterator + //TODO deal with multiple lines and add a for loop to handle them instead of accessing line 0 + vtkNew< vtkIdList > cellPts; + auto iter = ::vtk::TakeSmartPointer( polyData->GetLines()->NewIterator()); + iter->GetCellAtId( 0, cellPts ); + + const globalIndex nbSegments = cellPts->GetNumberOfIds() - 1; + m_segmentToPolyNodeMap.resizeDimension< 0 >( nbSegments ); + m_segmentToPolyNodeMap.resizeDimension< 1 >( m_numNodesPerElem ); + + globalIndex iseg = 0; + for( vtkIdType pointId = 0; pointId < nbSegments; ++pointId ) + { + m_segmentToPolyNodeMap[iseg][0] = cellPts->GetId( pointId ); + m_segmentToPolyNodeMap[iseg][1] = cellPts->GetId( pointId + 1 ); + ++iseg; + } + } +} + +REGISTER_CATALOG_ENTRY( WellGeneratorBase, VTKWellGenerator, string const &, Group * const ) +} diff --git a/src/coreComponents/mesh/generators/VTKWellGenerator.hpp b/src/coreComponents/mesh/generators/VTKWellGenerator.hpp new file mode 100644 index 00000000000..1290db3f7fe --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKWellGenerator.hpp @@ -0,0 +1,88 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file VTKWellGenerator.hpp + */ + +#ifndef GEOS_MESH_GENERATORS_VTKWELLGENERATOR_HPP +#define GEOS_MESH_GENERATORS_VTKWELLGENERATOR_HPP + +#include "mesh/generators/WellGeneratorBase.hpp" +#include "mesh/Perforation.hpp" +#include "dataRepository/Group.hpp" +#include "codingUtilities/Utilities.hpp" +#include "common/DataTypes.hpp" + +#include "mesh/generators/VTKUtilities.hpp" + +#include + +namespace geos +{ + +/** + * @class VTKWellGenerator + * @brief The VTKWellGenerator class provides a class implementation of VTK generated well. + */ +class VTKWellGenerator final : public WellGeneratorBase +{ +public: + + /** + * @name Constructor / Destructor + */ + ///@{ + + /** + * @brief Main constructor for VTKWellGenerator base class. + * @param[in] name of the VTKWellGenerator object + * @param[in] parent the parent Group pointer for the WellGenerator object + */ + VTKWellGenerator( const string & name, + Group * const parent ); + +/** + * @brief Return the name of the VTKWellGenerator in object Catalog. + * @return string that contains the key name to VTKWellGenerator in the Catalog + */ + static string catalogName() { return "VTKWell"; } + + ///@} + + +private: + + ///@cond DO_NOT_DOCUMENT + struct viewKeyStruct + { + constexpr static char const * filePathString() { return "file"; } + }; + /// @endcond + + + /** + * @brief + */ + void fillPolylineDataStructure( ) override; + + + /// Path to the mesh file + Path m_filePath; + +}; + +} // namespace geos + +#endif /* GEOS_MESH_GENERATORS_VTKWELLGENERATOR_HPP */ diff --git a/src/coreComponents/mesh/generators/WellGeneratorABC.hpp b/src/coreComponents/mesh/generators/WellGeneratorABC.hpp new file mode 100644 index 00000000000..e1b879de296 --- /dev/null +++ b/src/coreComponents/mesh/generators/WellGeneratorABC.hpp @@ -0,0 +1,204 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/* + * @file WellGeneratorABC.hpp + * + */ + +#ifndef GEOS_MESH_GENERATORS_WELLGENERATORABC_HPP_ +#define GEOS_MESH_GENERATORS_WELLGENERATORABC_HPP_ + +#include "dataRepository/Group.hpp" +#include "codingUtilities/Utilities.hpp" +#include "common/DataTypes.hpp" + + +namespace geos +{ + +/** + * @class WellGeneratorABC + * + * Abstract base class defining the information provided by any the well generator class. + */ +class WellGeneratorABC : public dataRepository::Group +{ +public: + + /** + * @brief Constructor. + * @param name name of the object in the data hierarchy. + * @param parent pointer to the parent group in the data hierarchy. + */ + WellGeneratorABC( const string & name, + Group * const parent ) + : + Group( name, parent ) + { } + + /** + * @brief Main function of the class that generates the well geometry + */ + virtual void generateWellGeometry( ) = 0; + + /** + * @name Getters / Setters + */ + ///@{ + + // getters for element data + + /** + * @brief Get the global number of well elements. + * @return the global number of elements + */ + virtual globalIndex numElements() const = 0; + + /** + * @brief Getter to the Segment to PolyNode mapping + * @return The Segment to PolyNode mapping as a 2D array + */ + virtual const array2d< globalIndex > & getSegmentToPolyNodeMap() const = 0; + + /** + * @brief Get the number of nodes per well element + * @return the number of nodes per well element + */ + virtual globalIndex numNodesPerElement() const = 0; + + /** + * @brief Get the Coordinates of the polyline nodes + * @return the Coordinates of the polyline nodes + */ + virtual const array2d< real64 > & getPolyNodeCoord() const = 0; + + /** + * @return The minimum segment length + */ + virtual real64 getMinSegmentLength() const = 0; + + /** + * @return The minimum element length + */ + virtual real64 getMinElemLength() const = 0; + + /** + * @return The list of perforation names + */ + virtual const string_array & getPerforationList() const = 0; + + /** + * @brief Get the physical location of the centers of well elements. + * @return list of center locations of the well elements + */ + virtual arrayView2d< real64 const > getElemCoords() const = 0; + + /** + * @brief Get the global indices mapping an element to the next. + * @return list providing the global index of the next element for each element + */ + virtual arrayView1d< globalIndex const > getNextElemIndex() const = 0; + + /** + * @brief Get the global indices mapping an element to the previous ones. + * @return list providing the global indices of the previous elements for each element + */ + virtual arrayView1d< arrayView1d< globalIndex const > const > getPrevElemIndices() const = 0; + + /** + * @brief Get the global indices of the well nodes nodes connected to each element. + * @return list providing the global index of the well nodes for each well element + */ + virtual arrayView2d< globalIndex const > getElemToNodesMap() const = 0; + + /** + * @brief Get the volume of the well elements. + * @return list of volumes of the well elements + */ + virtual arrayView1d< real64 const > getElemVolume() const = 0; + + /** + * @brief Get the radius in the well. + * @return the radius in the well + */ + virtual real64 getElementRadius() const = 0; + + // getters for node data + + /** + * @brief Get the global number of well nodes. + * @return the global number of nodes + */ + virtual globalIndex numNodes() const = 0; + + /** + * @brief Get the physical location of the centers of well elements. + * @return list of center locations of the well elements + */ + virtual arrayView2d< real64 const > getNodeCoords() const = 0; + + + + // getters for perforation data + /** + * @brief Get the global number of perforations on this well. + * @return the global number of elements + */ + virtual globalIndex numPerforations() const = 0; + + /** + * @brief Get the locations of the perforations. + * @return list of locations of all the perforations on the well + */ + virtual arrayView2d< real64 const > getPerfCoords() const = 0; + + /** + * @brief Get the well transmissibility at the perforations. + * @return list of well transmissibility at all the perforations on the well + */ + virtual arrayView1d< real64 const > getPerfTransmissibility() const = 0; + + /** + * @brief Get the skin factor at a perforation. + * @return the skin factor at a perforation + */ + virtual arrayView1d< real64 const > getPerfSkinFactor() const = 0; + + /** + * @brief Get the global indices of the well elements connected to each perforation. + * @return list providing the global index of the connected well element for each perforation + */ + virtual arrayView1d< globalIndex const > getPerfElemIndex() const = 0; + + /** + * @returns The number of physical dimensions + */ + virtual int getPhysicalDimensionsNumber() const = 0; + + /** + * Getter for the associated well region name + * @return the associated well region name + */ + virtual const string getWellRegionName() const = 0; + + /** + * Getter for the associated well control name + * @return the associated well control name + */ + virtual const string getWellControlsName() const = 0; + ///@} +}; +} +#endif /* GEOS_MESH_GENERATORS_WELLGENERATORABC_HPP_ */ diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp index 98f93ecef7f..554f4530bf7 100644 --- a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp +++ b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp @@ -13,17 +13,63 @@ */ #include "WellGeneratorBase.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" #include "mesh/Perforation.hpp" +#include "mesh/generators/LineBlockABC.hpp" +#include "LvArray/src/genericTensorOps.hpp" namespace geos { using namespace dataRepository; WellGeneratorBase::WellGeneratorBase( string const & name, Group * const parent ): - Group( name, parent ), - m_numPerforations( 0 ) + WellGeneratorABC( name, parent ) + , m_numPerforations( 0 ) + , m_numElemsPerSegment( 0 ) + , m_minSegmentLength( 1e-2 ) + , m_minElemLength( 1e-3 ) + , m_radius( 0 ) + , m_wellRegionName( "" ) + , m_wellControlsName( "" ) + , m_numElems( 0 ) + , m_numNodesPerElem( 2 ) + , m_numNodes( 0 ) + , m_nDims( 3 ) + , m_polylineHeadNodeId( -1 ) { setInputFlags( InputFlags::OPTIONAL_NONUNIQUE ); + + registerWrapper( viewKeyStruct::radiusString(), &m_radius ). + setInputFlag( InputFlags::REQUIRED ). + setSizedFromParent( 0 ). + setDescription( "Radius of the well [m]" ); + + registerWrapper( viewKeyStruct::numElementsPerSegmentString(), &m_numElemsPerSegment ). + setInputFlag( InputFlags::REQUIRED ). + setSizedFromParent( 0 ). + setDescription( "Number of well elements per polyline segment" ); + + registerWrapper( viewKeyStruct::minSegmentLengthString(), &m_minSegmentLength ). + setApplyDefaultValue( 1e-2 ). + setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( "Minimum length of a well segment [m]" ); + + registerWrapper( viewKeyStruct::minElementLengthString(), &m_minElemLength ). + setApplyDefaultValue( 1e-3 ). + setInputFlag( InputFlags::OPTIONAL ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( "Minimum length of a well element, computed as (segment length / number of elements per segment ) [m]" ); + + registerWrapper( viewKeyStruct::wellRegionNameString(), &m_wellRegionName ). + setInputFlag( InputFlags::REQUIRED ). + setSizedFromParent( 0 ). + setDescription( "Name of the well element region" ); + + registerWrapper( viewKeyStruct::wellControlsNameString(), &m_wellControlsName ). + setInputFlag( InputFlags::REQUIRED ). + setSizedFromParent( 0 ). + setDescription( "Name of the set of constraints associated with this well" ); } Group * WellGeneratorBase::createChild( string const & childKey, string const & childName ) @@ -34,7 +80,7 @@ Group * WellGeneratorBase::createChild( string const & childKey, string const & // keep track of the perforations that have been added m_perforationList.emplace_back( childName ); - + GEOS_LOG_RANK_0( "Adding Well attribute: " << childKey << ", " << childName ); return ®isterGroup< Perforation >( childName ); } else @@ -44,7 +90,6 @@ Group * WellGeneratorBase::createChild( string const & childKey, string const & return nullptr; } - void WellGeneratorBase::expandObjectCatalogs() { createChild( viewKeyStruct::perforationString(), viewKeyStruct::perforationString() ); @@ -56,4 +101,482 @@ WellGeneratorBase::CatalogInterface::CatalogType & WellGeneratorBase::getCatalog return catalog; } +void WellGeneratorBase::generateWellGeometry( ) +{ + fillPolylineDataStructure(); + + // count the number of well elements to create + m_numElems = m_numElemsPerSegment * m_segmentToPolyNodeMap.size( 0 ); + m_numNodes = m_numElems + 1; + + // resize the well element, node, and perforation arrays + m_elemCenterCoords.resize( m_numElems, 3 ); + m_nextElemId.resize( m_numElems ); + m_prevElemId.resize( m_numElems ); + m_elemToNodesMap.resizeDimension< 0 >( m_numElems ); + m_elemToNodesMap.resizeDimension< 1 >( m_numNodesPerElem ); + m_elemVolume.resize( m_numElems ); + + m_nodeDistFromHead.resize( m_numNodes ); + m_nodeCoords.resize( m_numNodes, 3 ); + + m_perfCoords.resize( m_numPerforations, 3 ); + m_perfDistFromHead.resize( m_numPerforations ); + m_perfTransmissibility.resize( m_numPerforations ); + m_perfSkinFactor.resize( m_numPerforations ); + m_perfElemId.resize( m_numPerforations ); + + // construct a reverse map from the polyline nodes to the segments + constructPolylineNodeToSegmentMap(); + + // detect the head polyline node based on depth + findPolylineHeadNodeIndex(); + + // compute the location and distance from top of the well elements + discretizePolyline(); + + // map the perforations to the well elements + connectPerforationsToWellElements(); + + // make sure that the perforation locations are valid + checkPerforationLocationsValidity(); + + if( getLogLevel() >= 1 ) + { + debugWellGeometry(); + } + +} + +void WellGeneratorBase::postProcessInput() +{ + GEOS_THROW_IF( m_radius <= 0, + "Invalid " << viewKeyStruct::radiusString() << " in well " << getName(), + InputError ); + + GEOS_THROW_IF( m_wellRegionName.empty(), + "Invalid well region name in well " << getName(), + InputError ); + + GEOS_THROW_IF( m_wellControlsName.empty(), + "Invalid well constraint name in well " << getName(), + InputError ); +} + +void WellGeneratorBase::constructPolylineNodeToSegmentMap() +{ + m_polyNodeToSegmentMap.resize( m_polyNodeCoords.size( 0 ) ); + + bool foundSmallElem = false; + + // loop over the segments + for( globalIndex iseg = 0; iseg < m_segmentToPolyNodeMap.size( 0 ); ++iseg ) + { + globalIndex const ipolyNode_a = m_segmentToPolyNodeMap[iseg][0]; + globalIndex const ipolyNode_b = m_segmentToPolyNodeMap[iseg][1]; + + real64 vSeg[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_polyNodeCoords[ipolyNode_a] ); + LvArray::tensorOps::subtract< 3 >( vSeg, m_polyNodeCoords[ipolyNode_b] ); + + real64 const segmentLength = LvArray::tensorOps::l2Norm< 3 >( vSeg ); + + // various checks and warnings on the segment and element length + GEOS_THROW_IF( segmentLength < m_minSegmentLength, + "Error in the topology of well '" << getName() << + "': we detected a polyline segment measuring less than " << m_minSegmentLength << "m. \n" << + "You can change the minimum segment length using the field " << viewKeyStruct::minSegmentLengthString(), + InputError ); + + GEOS_THROW_IF( m_polyNodeCoords[ipolyNode_a][2] < m_polyNodeCoords[ipolyNode_b][2], + "Error in the topology of well '" << getName() << + "': in the polyline, each segment must be going down. \n" << + "This is not the case between polyline nodes " << m_polyNodeCoords[ipolyNode_a] << " and " << m_polyNodeCoords[ipolyNode_b], + InputError ); + + if( segmentLength / m_numElemsPerSegment < m_minElemLength ) + { + foundSmallElem = true; + } + + // map the polyline node ids to the polyline segment ids + m_polyNodeToSegmentMap[ipolyNode_a].insert( iseg ); + m_polyNodeToSegmentMap[ipolyNode_b].insert( iseg ); + } + + if( foundSmallElem ) + { + GEOS_LOG_RANK_0( "\nWarning: the chosen number of well elements per polyline segment (" << m_numElemsPerSegment << + ") leads to well elements measuring less than " << m_minElemLength << "m in the topology of well '" << getName() << "'.\n" << + "The simulation can proceed like that, but small well elements may cause numerical issues, so it is something to keep an eye on.\n" << + "You can get rid of this message by changing the field " << viewKeyStruct::minElementLengthString() ); + } +} + +void WellGeneratorBase::findPolylineHeadNodeIndex() +{ + // we assume that the first *polyline segment* in the XML input is the well head segment + globalIndex const wellHeadSegId = 0; + + // get the corresponding node indices + globalIndex const ipolyNode_a = m_segmentToPolyNodeMap[wellHeadSegId][0]; + globalIndex const ipolyNode_b = m_segmentToPolyNodeMap[wellHeadSegId][1]; + + // we determine which node is the head node based on depth + // therefore here we throw an error if the well head segment is horizontal + GEOS_THROW_IF( !(m_polyNodeCoords[ipolyNode_a][2] < m_polyNodeCoords[ipolyNode_b][2]) + && !(m_polyNodeCoords[ipolyNode_a][2] > m_polyNodeCoords[ipolyNode_b][2]), + "The head polyline segment cannot be horizontal in well '" << getName() + << "' since we use depth to determine which of its nodes is to head node of the well.\n" + << "If you are trying to set up a horizontal well, please simply add a non-horizontal segment at the top of the well," + << " and this error will go away", + InputError ); + + // detect the top node, assuming z oriented upwards + m_polylineHeadNodeId = + ( m_polyNodeCoords[ipolyNode_a][2] > m_polyNodeCoords[ipolyNode_b][2] ) + ? ipolyNode_a + : ipolyNode_b; + + real64 const headZcoord = m_polyNodeCoords[m_polylineHeadNodeId][2]; + for( globalIndex inode = 0; inode < m_polyNodeCoords.size( 0 ); ++inode ) + { + if( inode == m_polylineHeadNodeId ) + { + continue; + } + real64 const currentZcoord = m_polyNodeCoords[inode][2]; + + GEOS_THROW_IF( !(currentZcoord < headZcoord), + "Error in the topology of well '" << getName() + << "' since we found a well node that is above the head node", + InputError ); + } +} + +void WellGeneratorBase::discretizePolyline() +{ + // initialize well elements and node ids + globalIndex ipolyNodeTop = m_polylineHeadNodeId; + globalIndex iwelemCurrent = 0; + globalIndex isegCurrent = 0; + + // set the location of the first well node and distance from well head + LvArray::tensorOps::copy< 3 >( m_nodeCoords[0], m_polyNodeCoords[ipolyNodeTop] ); + m_nodeDistFromHead[0] = 0.0; + + // note: this part of the code does not support well branching + // TODO: check that there is only one branch + // TODO: read wells with branching (already supported elsewhere in the code) + + // go through the well from top to bottom + for( globalIndex is = 0; is < m_segmentToPolyNodeMap.size( 0 ); ++is ) + { + + GEOS_THROW_IF( isegCurrent == -1, + "Invalid segmentToNode map in well " << getName(), + InputError ); + + globalIndex const ipolyNodeBottom = ( ipolyNodeTop == m_segmentToPolyNodeMap[isegCurrent][0] ) + ? m_segmentToPolyNodeMap[isegCurrent][1] + : m_segmentToPolyNodeMap[isegCurrent][0]; + + real64 vPoly[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_polyNodeCoords[ipolyNodeBottom] ); + LvArray::tensorOps::subtract< 3 >( vPoly, m_polyNodeCoords[ipolyNodeTop] ); + + // add the well elements and well nodes corresponding to this polyline segment + for( localIndex iw = 0; iw < m_numElemsPerSegment; ++iw ) + { + + // 1) set the element location, connectivity + real64 const scaleCenter = (iw + 0.5) / static_cast< real64 >(m_numElemsPerSegment); + LvArray::tensorOps::copy< 3 >( m_elemCenterCoords[iwelemCurrent], vPoly ); + LvArray::tensorOps::scale< 3 >( m_elemCenterCoords[iwelemCurrent], scaleCenter ); + LvArray::tensorOps::add< 3 >( m_elemCenterCoords[iwelemCurrent], m_polyNodeCoords[ipolyNodeTop] ); + + GEOS_THROW_IF( iwelemCurrent >= m_numElems, + "Invalid well topology in well " << getName(), + InputError ); + + globalIndex const iwelemTop = iwelemCurrent - 1; + globalIndex const iwelemBottom = iwelemCurrent + 1; + m_nextElemId[iwelemCurrent] = iwelemTop; + m_prevElemId[iwelemCurrent].resize( 1 ); + m_prevElemId[iwelemCurrent][0] = (iwelemBottom < m_numElems) + ? iwelemBottom + : -1; + + // 2) set the node properties + globalIndex const iwellNodeTop = iwelemCurrent; + globalIndex const iwellNodeBottom = iwelemCurrent+1; + m_elemToNodesMap[iwelemCurrent][LineBlockABC::NodeLocation::TOP] = iwellNodeTop; + m_elemToNodesMap[iwelemCurrent][LineBlockABC::NodeLocation::BOTTOM] = iwellNodeBottom; + + real64 const scaleBottom = (iw + 1.0) / m_numElemsPerSegment; + LvArray::tensorOps::copy< 3 >( m_nodeCoords[iwellNodeBottom], vPoly ); + LvArray::tensorOps::scale< 3 >( m_nodeCoords[iwellNodeBottom], scaleBottom ); + LvArray::tensorOps::add< 3 >( m_nodeCoords[iwellNodeBottom], m_polyNodeCoords[ipolyNodeTop] ); + + // 3) increment the distance from the well head to the bottom of current element + real64 vWellElem[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( m_nodeCoords[iwellNodeBottom] ); + LvArray::tensorOps::subtract< 3 >( vWellElem, m_nodeCoords[iwellNodeTop] ); + m_nodeDistFromHead[iwellNodeBottom] = LvArray::tensorOps::l2Norm< 3 >( vWellElem ); + m_nodeDistFromHead[iwellNodeBottom] += m_nodeDistFromHead[iwellNodeTop]; + + // 4) set element volume + m_elemVolume[iwelemCurrent] = LvArray::tensorOps::l2Norm< 3 >( vWellElem ) * M_PI * m_radius * m_radius; + + // 4) increment the element counter + ++iwelemCurrent; + } + + // then consider the next polyline segment + ipolyNodeTop = ipolyNodeBottom; + isegCurrent = getNextSegmentIndex( isegCurrent, ipolyNodeTop ); + } + + // set the previous index for the bottom segment + m_prevElemId[iwelemCurrent-1].resize( 1 ); + m_prevElemId[iwelemCurrent-1][0] = -1; + +} + + +void WellGeneratorBase::connectPerforationsToWellElements() +{ + + // assign a well element to each perforation + for( globalIndex iperf = 0; iperf < m_numPerforations; ++iperf ) + { + + // get the perforation and its properties + Perforation const & perf = this->getGroup< Perforation >( m_perforationList[iperf] ); + m_perfDistFromHead[iperf] = perf.getDistanceFromWellHead(); + m_perfTransmissibility[iperf] = perf.getWellTransmissibility(); + m_perfSkinFactor[iperf] = perf.getWellSkinFactor(); + + // search in all the elements of this well between head and bottom + globalIndex iwelemTop = 0; + globalIndex iwelemBottom = m_numElems - 1; + + // check the validity of the perforation before starting + real64 const wellLength = m_nodeDistFromHead[m_elemToNodesMap[iwelemBottom][LineBlockABC::NodeLocation::BOTTOM]]; + + GEOS_THROW_IF( m_perfDistFromHead[iperf] > wellLength, + "Distance from perforation " << perf.getName() << " to head is larger than well polyline length for well " << getName() << "\n \n" + << "Here is how the \"distanceFromHead\" keyword is used in the definition of the perforation location: \n" + << "We start from the well head (top of the well) and we measure the linear distance along the well polyline as we go down the well.\n" + << "When we reach the distanceFromHead specified by the user, we place a perforation on the well at this location of the polyline, and connect it to the reservoir element that contains this perforation", + InputError ); + + // start binary search + const globalIndex maxNumSteps = m_numElems + 1; + globalIndex currentNumSteps = 0; + while( iwelemTop < iwelemBottom ) + { + globalIndex iwelemMid = + static_cast< globalIndex >(floor( static_cast< real64 >(iwelemTop + iwelemBottom) / 2.0 )); + real64 const headToBottomDist = + m_nodeDistFromHead[m_elemToNodesMap[iwelemMid][LineBlockABC::NodeLocation::BOTTOM]]; + + if( headToBottomDist < m_perfDistFromHead[iperf] ) + { + iwelemTop = iwelemMid + 1; + } + else + { + iwelemBottom = iwelemMid; + } + + GEOS_THROW_IF( currentNumSteps > maxNumSteps, + "Perforation " << perf.getName() << " cannot be mapped to a well element in well " << getName(), + InputError ); + + currentNumSteps++; + } + + // set the index of the matched element + globalIndex iwelemMatched = iwelemTop; + GEOS_THROW_IF( iwelemMatched >= m_numElems, + "Invalid topology in well " << getName(), + InputError ); + + m_perfElemId[iperf] = iwelemMatched; + + // compute the physical location of the perforation + globalIndex const inodeTop = m_elemToNodesMap[iwelemMatched][LineBlockABC::NodeLocation::TOP]; + globalIndex const inodeBottom = m_elemToNodesMap[iwelemMatched][LineBlockABC::NodeLocation::BOTTOM]; + real64 const elemLength = m_nodeDistFromHead[inodeBottom] - m_nodeDistFromHead[inodeTop]; + real64 const topToPerfDist = m_perfDistFromHead[iperf] - m_nodeDistFromHead[inodeTop]; + + GEOS_THROW_IF( (elemLength <= 0) || (topToPerfDist < 0), + "Invalid topology in well " << getName(), + InputError ); + + LvArray::tensorOps::copy< 3 >( m_perfCoords[iperf], m_nodeCoords[inodeBottom] ); + LvArray::tensorOps::subtract< 3 >( m_perfCoords[iperf], m_nodeCoords[inodeTop] ); + LvArray::tensorOps::scale< 3 >( m_perfCoords[iperf], topToPerfDist / elemLength ); + LvArray::tensorOps::add< 3 >( m_perfCoords[iperf], m_nodeCoords[inodeTop] ); + + } +} + + +globalIndex WellGeneratorBase::getNextSegmentIndex( globalIndex topSegId, + globalIndex currentPolyNodeId ) const +{ + globalIndex nextSegId = -1; + + // get the index of the two segments sharing this node + for( auto is : m_polyNodeToSegmentMap[currentPolyNodeId] ) + { + // if this is not equal to the current segId, we found the next segId + if( is != topSegId ) + { + nextSegId = is; + break; + } + } + + return nextSegId; +} + +void WellGeneratorBase::checkPerforationLocationsValidity() +{ + array1d< array1d< localIndex > > elemToPerfMap; + elemToPerfMap.resize( m_numElems ); + + for( globalIndex iperf = 0; iperf < m_numPerforations; ++iperf ) + { + elemToPerfMap[m_perfElemId[iperf]].emplace_back( iperf ); + } + + // merge perforations to make sure that no well element is shared between two MPI domains + // TODO: instead of merging perforations, split the well elements and do not change the physical location of the + // perforation + int const mpiSize = MpiWrapper::commSize( MPI_COMM_GEOSX ); + if( mpiSize > 1 ) + { + mergePerforations( elemToPerfMap ); + } + + for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) + { + // check that there is always a perforation in the last well element (otherwise, the problem is not well posed) + for( localIndex iwelemPrev = 0; iwelemPrev < m_prevElemId[iwelem].size(); ++iwelemPrev ) + { + GEOS_THROW_IF( m_prevElemId[iwelem][iwelemPrev] == -1 && elemToPerfMap[iwelem].size() == 0, + "The bottom element of well " << getName() << " does not have a perforation. " + << "This is needed to have a well-posed problem. \n\n" + << "Here are the two possible ways to solve this problem: \n\n" + << "1) Adding a perforation located close to the bottom of the well. " + << "To do that, compute the total length of the well polyline (by summing the length of the well segments defined by the keywords \"polylineNodeCoords\" and \"polylineSegmentConn\") " + << "and place a perforation whose \"distanceFromHead\" is slightly smaller than this total length. \n \n" + << "2) Shorten the well polyline. " + << "To do that, reduce the length of the well polyline by shortening the segments defined by the keywords \"polylineNodeCoords\" and \"polylineSegmentConn\", or by removing a segment.", + InputError ); + } + } +} + +void WellGeneratorBase::mergePerforations( array1d< array1d< localIndex > > const & elemToPerfMap ) +{ + + for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) + { + // collect the indices of the elems with more that one perforation + if( elemToPerfMap[iwelem].size() > 1 ) + { + // find the perforation with the largest Peaceman index and keep its location + globalIndex iperfMaxTransmissibility = elemToPerfMap[iwelem][0]; + real64 maxTransmissibility = m_perfTransmissibility[iperfMaxTransmissibility]; + for( localIndex ip = 1; ip < elemToPerfMap[iwelem].size(); ++ip ) + { + if( m_perfTransmissibility[elemToPerfMap[iwelem][ip]] > maxTransmissibility ) + { + iperfMaxTransmissibility = elemToPerfMap[iwelem][ip]; + maxTransmissibility = m_perfTransmissibility[iperfMaxTransmissibility]; + } + } + + // assign the coordinates of the perf with the largest trans to the other perfs on this elem + for( localIndex ip = 0; ip < elemToPerfMap[iwelem].size(); ++ip ) + { + if( elemToPerfMap[iwelem][ip] == iperfMaxTransmissibility ) + { + continue; + } + + GEOS_LOG_RANK_0( "\n \nThe GEOSX wells currently have the following limitation in parallel: \n" + << "We cannot allow an element of the well mesh to have two or more perforations associated with it. \n" + << "So, in the present simulation, perforation #" << elemToPerfMap[iwelem][ip] + << " of well " << getName() + << " is moved from " << m_perfCoords[elemToPerfMap[iwelem][ip]] + << " to " << m_perfCoords[iperfMaxTransmissibility] + << " to make sure that no element of the well mesh has two perforations associated with it. \n" + << "To circumvent this issue, please increase the value of \"numElementsPerSegment\" that controls the number of (uniformly distributed) well mesh elements per segment of the well polyline. \n" + << "Our recommendation is to choose \"numElementsPerSegment\" such that each well mesh element has at most one perforation. \n\n" ); + LvArray::tensorOps::copy< 3 >( m_perfCoords[elemToPerfMap[iwelem][ip]], m_perfCoords[iperfMaxTransmissibility] ); + } + } + } +} + +void WellGeneratorBase::debugWellGeometry() const +{ + if( MpiWrapper::commRank( MPI_COMM_GEOSX ) != 0 ) + { + return; + } + + std::cout << std::endl; + std::cout << "++++++++++++++++++++++++++" << std::endl; + std::cout << "WellGeneratorBase = " << getName() << std::endl; + std::cout << "MPI rank = " << MpiWrapper::commRank( MPI_COMM_GEOSX ) << std::endl << std::endl; + std::cout << "Number of well elements = " << m_numElems << std::endl; + + for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) + { + std::cout << "Well element #" << iwelem << std::endl; + std::cout << "Coordinates of the element center: " << m_elemCenterCoords[iwelem] << std::endl; + if( m_nextElemId[iwelem] < 0 ) + { + std::cout << "No next well element" << std::endl; + } + else + { + std::cout << "Next well element # = " << m_nextElemId[iwelem] << std::endl; + } + if( m_prevElemId[iwelem][0] < 0 ) + { + std::cout << "No previous well element" << std::endl; + } + else + { + std::cout << "Previous well element #" << m_prevElemId[iwelem][0] << std::endl; + } + for( globalIndex inode = 0; inode < m_numNodesPerElem; ++inode ) + { + if( inode == 0 ) + { + std::cout << "First well node: #" << m_elemToNodesMap[iwelem][inode] << std::endl; + } + else + { + std::cout << "Second well node: #" << m_elemToNodesMap[iwelem][inode] << std::endl; + } + } + } + + std::cout << std::endl << "Number of perforations = " << m_numPerforations << std::endl; + + for( globalIndex iperf = 0; iperf < m_numPerforations; ++iperf ) + { + std::cout << "Perforation #" << iperf << std::endl; + std::cout << "Coordinates of the perforation: " << m_perfCoords[iperf] << std::endl; + std::cout << "Is connected to well element #" << m_perfElemId[iperf] << std::endl; + } + std::cout << std::endl; + +} + } diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.hpp b/src/coreComponents/mesh/generators/WellGeneratorBase.hpp index e24cafd8cb7..a4c6cb5f577 100644 --- a/src/coreComponents/mesh/generators/WellGeneratorBase.hpp +++ b/src/coreComponents/mesh/generators/WellGeneratorBase.hpp @@ -20,6 +20,7 @@ #ifndef GEOS_MESH_GENERATORS_WELLGENERATORBASE_HPP_ #define GEOS_MESH_GENERATORS_WELLGENERATORBASE_HPP_ +#include "mesh/generators/WellGeneratorABC.hpp" #include "dataRepository/Group.hpp" #include "codingUtilities/Utilities.hpp" #include "common/DataTypes.hpp" @@ -33,7 +34,7 @@ namespace geos * * This class processes the data of a single well from the XML and generates the well geometry */ -class WellGeneratorBase : public dataRepository::Group +class WellGeneratorBase : public WellGeneratorABC { public: @@ -75,7 +76,8 @@ class WellGeneratorBase : public dataRepository::Group /** * @brief Main function of the class that generates the well geometry */ - virtual void generateWellGeometry( ) = 0; + void generateWellGeometry( ) override; + /** * @name Getters / Setters @@ -88,76 +90,76 @@ class WellGeneratorBase : public dataRepository::Group * @brief Get the global number of well elements. * @return the global number of elements */ - virtual globalIndex numElements() const = 0; + globalIndex numElements() const override { return m_numElems; } /** * @brief Getter to the Segment to PolyNode mapping * @return The Segment to PolyNode mapping as a 2D array */ - virtual const array2d< globalIndex > & getSegmentToPolyNodeMap() const = 0; + const array2d< globalIndex > & getSegmentToPolyNodeMap() const override { return m_segmentToPolyNodeMap; }; /** * @brief Get the number of nodes per well element * @return the number of nodes per well element */ - virtual globalIndex numNodesPerElement() const = 0; + globalIndex numNodesPerElement() const override { return m_numNodesPerElem; } /** * @brief Get the Coordinates of the polyline nodes * @return the Coordinates of the polyline nodes */ - virtual const array2d< real64 > & getPolyNodeCoord() const = 0; + const array2d< real64 > & getPolyNodeCoord() const override { return m_polyNodeCoords; } /** * @return The minimum segment length */ - virtual real64 getMinSegmentLength() const = 0; + real64 getMinSegmentLength() const override { return m_minSegmentLength; } /** * @return The minimum element length */ - virtual real64 getMinElemLength() const = 0; + real64 getMinElemLength() const override { return m_minElemLength; } /** * @return The list of perforation names */ - const string_array & getPerforationList() const { return m_perforationList; } + const string_array & getPerforationList() const override { return m_perforationList; } /** * @brief Get the physical location of the centers of well elements. * @return list of center locations of the well elements */ - virtual arrayView2d< real64 const > getElemCoords() const = 0; + arrayView2d< real64 const > getElemCoords() const override { return m_elemCenterCoords; } /** * @brief Get the global indices mapping an element to the next. * @return list providing the global index of the next element for each element */ - virtual arrayView1d< globalIndex const > getNextElemIndex() const = 0; + arrayView1d< globalIndex const > getNextElemIndex() const override { return m_nextElemId; } /** * @brief Get the global indices mapping an element to the previous ones. * @return list providing the global indices of the previous elements for each element */ - virtual arrayView1d< arrayView1d< globalIndex const > const > getPrevElemIndices() const = 0; + arrayView1d< arrayView1d< globalIndex const > const > getPrevElemIndices() const override { return m_prevElemId.toNestedViewConst(); } /** * @brief Get the global indices of the well nodes nodes connected to each element. * @return list providing the global index of the well nodes for each well element */ - virtual arrayView2d< globalIndex const > getElemToNodesMap() const = 0; + arrayView2d< globalIndex const > getElemToNodesMap() const override { return m_elemToNodesMap; } /** * @brief Get the volume of the well elements. * @return list of volumes of the well elements */ - virtual arrayView1d< real64 const > getElemVolume() const = 0; + arrayView1d< real64 const > getElemVolume() const override { return m_elemVolume; } /** * @brief Get the radius in the well. * @return the radius in the well */ - virtual real64 getElementRadius() const = 0; + real64 getElementRadius() const override { return m_radius; } // getters for node data @@ -165,63 +167,61 @@ class WellGeneratorBase : public dataRepository::Group * @brief Get the global number of well nodes. * @return the global number of nodes */ - virtual globalIndex numNodes() const = 0; + globalIndex numNodes() const override { return m_numNodes; } /** * @brief Get the physical location of the centers of well elements. * @return list of center locations of the well elements */ - virtual arrayView2d< real64 const > getNodeCoords() const = 0; - - + arrayView2d< real64 const > getNodeCoords() const override { return m_nodeCoords; } // getters for perforation data /** * @brief Get the global number of perforations on this well. * @return the global number of elements */ - globalIndex numPerforations() const { return m_numPerforations; } + globalIndex numPerforations() const override { return m_numPerforations; } /** * @brief Get the locations of the perforations. * @return list of locations of all the perforations on the well */ - virtual arrayView2d< real64 const > getPerfCoords() const = 0; + arrayView2d< real64 const > getPerfCoords() const override { return m_perfCoords; } /** * @brief Get the well transmissibility at the perforations. * @return list of well transmissibility at all the perforations on the well */ - virtual arrayView1d< real64 const > getPerfTransmissibility() const = 0; + arrayView1d< real64 const > getPerfTransmissibility() const override { return m_perfTransmissibility; } /** * @brief Get the skin factor at a perforation. * @return the skin factor at a perforation */ - virtual arrayView1d< real64 const > getPerfSkinFactor() const = 0; + arrayView1d< real64 const > getPerfSkinFactor() const override { return m_perfSkinFactor; }; /** * @brief Get the global indices of the well elements connected to each perforation. * @return list providing the global index of the connected well element for each perforation */ - virtual arrayView1d< globalIndex const > getPerfElemIndex() const = 0; + arrayView1d< globalIndex const > getPerfElemIndex() const override { return m_perfElemId; } /** * @returns The number of physical dimensions */ - virtual int getPhysicalDimensionsNumber() const = 0; + int getPhysicalDimensionsNumber() const override { return m_nDims; } /** * Getter for the associated well region name * @return the associated well region name */ - virtual const string getWellRegionName() const = 0; + const string getWellRegionName() const override { return m_wellRegionName; } /** * Getter for the associated well control name * @return the associated well control name */ - virtual const string getWellControlsName() const = 0; + const string getWellControlsName() const override { return m_wellControlsName; } ///@cond DO_NOT_DOCUMENT struct viewKeyStruct @@ -234,17 +234,179 @@ class WellGeneratorBase : public dataRepository::Group constexpr static char const * radiusString() { return "radius"; } constexpr static char const * wellRegionNameString() { return "wellRegionName"; } constexpr static char const * wellControlsNameString() { return "wellControlsName"; } - constexpr static char const * meshNameString() { return "meshName"; } constexpr static char const * perforationString() { return "Perforation"; } }; /// @endcond protected: + + /** + * @brief This function provides capability to post process input values prior to + * any other initialization operations. + */ + void postProcessInput() override; + + /** + * @name Helper functions to construct the geometry of the well + */ + ///@{ + + /** + * @brief Fills the intermediate polyline data structure. + */ + virtual void fillPolylineDataStructure() { }; + + /** + * @brief Map each polyline node to the polyline segment(s) it is connected to. + */ + void constructPolylineNodeToSegmentMap(); + + /** + * @brief Find the head node of the well (i.e., top node of the polyline). + */ + void findPolylineHeadNodeIndex(); + + /** + * @brief Discretize the polyline by placing well elements. + */ + void discretizePolyline(); + + /** + * @brief Map each perforation to a well element. + */ + void connectPerforationsToWellElements(); + + /** + * @brief Make sure that the perforation locations are valid: + * - for partitioning purposes + * - to have a well-posed problem + */ + void checkPerforationLocationsValidity(); + + /** + * @brief Merge perforations on the elements with multiple perforations. + * @param[in] elemToPerfMap Connectivity between well elements and Perforations + */ + void mergePerforations( array1d< array1d< localIndex > > const & elemToPerfMap ); + + /** + * @brief At a given node, find the next segment going in the direction of the bottom of the well. + * @param[in] topSegId index of the top segment + * @param[in] currentNodeId index of the current node + * @return The Id of the next segment + */ + globalIndex getNextSegmentIndex( globalIndex topSegId, + globalIndex currentNodeId ) const; + + ///@} + + /// @cond DO_NOT_DOCUMENT + void debugWellGeometry() const; + /// @endcond + /// Global number of perforations globalIndex m_numPerforations; /// List of perforation names string_array m_perforationList; + + // XML Input + + /// Connectivity between the polyline nodes + array2d< globalIndex > m_segmentToPolyNodeMap; + + /// Number of well elements per polyline interval + int m_numElemsPerSegment; + + /// Min segment length + real64 m_minSegmentLength; + + /// Min well element length + real64 m_minElemLength; + + /// Radius area of the well (assumed to be valid for the entire well) + real64 m_radius; + + /// Name of the corresponding well region + string m_wellRegionName; + + /// Name of the constraints associated with this well + string m_wellControlsName; + + + + // Geometry of the well (later passed to the WellElementSubRegion) + + // well element data + + /// Global number of well elements + globalIndex m_numElems; + + /// Physical location of the center of the well element + array2d< real64 > m_elemCenterCoords; + + /// Global index of the next well element + array1d< globalIndex > m_nextElemId; + + /// Global indices of the prev well elements (maybe need multiple prevs for branching) + array1d< array1d< globalIndex > > m_prevElemId; + + /// Connectivity between elements and nodes + array2d< globalIndex > m_elemToNodesMap; + + /// Volume of well elements + array1d< real64 > m_elemVolume; + + + // well node data + + /// Number of nodes per well element + globalIndex const m_numNodesPerElem; + + /// Global number of well nodes + globalIndex m_numNodes; + + /// Physical location of the nodes + array2d< real64 > m_nodeCoords; + + // perforation data + + /// Absolute physical location of the perforation + array2d< real64 > m_perfCoords; + + /// Well Peaceman index at the perforation + array1d< real64 > m_perfTransmissibility; + + /// Global index of the well element + array1d< globalIndex > m_perfElemId; + + + + // Auxiliary data + + /// Number of physical dimensions + const int m_nDims; + + /// Coordinates of the polyline nodes + array2d< real64 > m_polyNodeCoords; + + /// Map from the polyline nodes to the polyline nodes + array1d< SortedArray< globalIndex > > m_polyNodeToSegmentMap; + + /// Index of the node at the well head + globalIndex m_polylineHeadNodeId; + + /// Physical location of the polyline node wrt to well head + array1d< real64 > m_nodeDistFromHead; + + // Perforation data + + /// Skin Factor at the perforation + array1d< real64 > m_perfSkinFactor; + + /// Physical location of the perforation wrt to well head + array1d< real64 > m_perfDistFromHead; + }; } #endif /* GEOS_MESH_GENERATORS_WELLGENERATORBASE_HPP_ */ diff --git a/src/coreComponents/schema/docs/AcousticElasticSEM.rst b/src/coreComponents/schema/docs/AcousticElasticSEM.rst new file mode 100644 index 00000000000..8a99bd42d0f --- /dev/null +++ b/src/coreComponents/schema/docs/AcousticElasticSEM.rst @@ -0,0 +1,18 @@ + + +========================= ================== ======== ======================================================================================================================================================================================================================================================================================================================== +Name Type Default Description +========================= ================== ======== ======================================================================================================================================================================================================================================================================================================================== +acousticSolverName groupNameRef required Name of the acoustic solver used by the coupled solver +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] +discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +elasticSolverName groupNameRef required Name of the elastic solver used by the coupled solver +initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +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. +LinearSolverParameters node unique :ref:`XML_LinearSolverParameters` +NonlinearSolverParameters node unique :ref:`XML_NonlinearSolverParameters` +========================= ================== ======== ======================================================================================================================================================================================================================================================================================================================== + + diff --git a/src/coreComponents/schema/docs/AcousticElasticSEM_other.rst b/src/coreComponents/schema/docs/AcousticElasticSEM_other.rst new file mode 100644 index 00000000000..3f1b37c3d6c --- /dev/null +++ b/src/coreComponents/schema/docs/AcousticElasticSEM_other.rst @@ -0,0 +1,13 @@ + + +========================= ============================================================================================================================================================== ================================================================ +Name Type Description +========================= ============================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geos_mapBase< std_pair< string, string >, LvArray_Array< string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +SolverStatistics node :ref:`DATASTRUCTURE_SolverStatistics` +========================= ============================================================================================================================================================== ================================================================ + + diff --git a/src/coreComponents/schema/docs/InternalMesh.rst b/src/coreComponents/schema/docs/InternalMesh.rst index 64d2b3e87c0..c8de2fad08c 100644 --- a/src/coreComponents/schema/docs/InternalMesh.rst +++ b/src/coreComponents/schema/docs/InternalMesh.rst @@ -18,6 +18,7 @@ yCoords real64_array required y-coordinates of each mesh block v zBias real64_array {1} Bias of element sizes in the z-direction within each mesh block (dz_left=(1+b)*L/N, dz_right=(1-b)*L/N) zCoords real64_array required z-coordinates of each mesh block vertex InternalWell node :ref:`XML_InternalWell` +VTKWell node :ref:`XML_VTKWell` ================= ================== ======== ======================================================================================================= diff --git a/src/coreComponents/schema/docs/InternalMesh_other.rst b/src/coreComponents/schema/docs/InternalMesh_other.rst index aef746a9df7..8ca4f02bc36 100644 --- a/src/coreComponents/schema/docs/InternalMesh_other.rst +++ b/src/coreComponents/schema/docs/InternalMesh_other.rst @@ -4,6 +4,7 @@ Name Type Description ============ ==== ================================= InternalWell node :ref:`DATASTRUCTURE_InternalWell` +VTKWell node :ref:`DATASTRUCTURE_VTKWell` meshLevels node :ref:`DATASTRUCTURE_meshLevels` ============ ==== ================================= diff --git a/src/coreComponents/schema/docs/InternalWell.rst b/src/coreComponents/schema/docs/InternalWell.rst index 1dd1f482f5d..bb26958fc59 100644 --- a/src/coreComponents/schema/docs/InternalWell.rst +++ b/src/coreComponents/schema/docs/InternalWell.rst @@ -11,8 +11,8 @@ numElementsPerSegment integer required Number of well elements per p polylineNodeCoords real64_array2d required Physical coordinates of the well polyline nodes polylineSegmentConn globalIndex_array2d required Connectivity of the polyline segments radius real64 required Radius of the well [m] -wellControlsName groupNameRef required Name of the set of constraints associated with this well -wellRegionName groupNameRef required Name of the well element region +wellControlsName string required Name of the set of constraints associated with this well +wellRegionName string required Name of the well element region Perforation node :ref:`XML_Perforation` ===================== =================== ======== ==================================================================================================== diff --git a/src/coreComponents/schema/docs/InternalWellbore.rst b/src/coreComponents/schema/docs/InternalWellbore.rst index 07cc08bc83e..d9e3151e2f9 100644 --- a/src/coreComponents/schema/docs/InternalWellbore.rst +++ b/src/coreComponents/schema/docs/InternalWellbore.rst @@ -24,6 +24,7 @@ yBias real64_array {1} Bias of element sizes in zBias real64_array {1} Bias of element sizes in the z-direction within each mesh block (dz_left=(1+b)*L/N, dz_right=(1-b)*L/N) zCoords real64_array required z-coordinates of each mesh block vertex InternalWell node :ref:`XML_InternalWell` +VTKWell node :ref:`XML_VTKWell` =========================== ================== ======== ============================================================================================================================================================================================================================ diff --git a/src/coreComponents/schema/docs/InternalWellbore_other.rst b/src/coreComponents/schema/docs/InternalWellbore_other.rst index 24b84b089e7..14901e94dac 100644 --- a/src/coreComponents/schema/docs/InternalWellbore_other.rst +++ b/src/coreComponents/schema/docs/InternalWellbore_other.rst @@ -8,6 +8,7 @@ ny integer_array Number of elements in the y-direction within each mes xCoords real64_array x-coordinates of each mesh block vertex yCoords real64_array y-coordinates of each mesh block vertex InternalWell node :ref:`DATASTRUCTURE_InternalWell` +VTKWell node :ref:`DATASTRUCTURE_VTKWell` meshLevels node :ref:`DATASTRUCTURE_meshLevels` ============ ============= ============================================================ diff --git a/src/coreComponents/schema/docs/PorousDamageElasticIsotropic.rst b/src/coreComponents/schema/docs/PorousDamageElasticIsotropic.rst new file mode 100644 index 00000000000..7d6d29602c1 --- /dev/null +++ b/src/coreComponents/schema/docs/PorousDamageElasticIsotropic.rst @@ -0,0 +1,13 @@ + + +============================ ============ ======== =========================================== +Name Type Default Description +============================ ============ ======== =========================================== +name groupName required A name is required for any non-unique nodes +permeabilityModelName groupNameRef required Name of the permeability model. +porosityModelName groupNameRef required Name of the porosity model. +solidInternalEnergyModelName groupNameRef Name of the solid internal energy model. +solidModelName groupNameRef required Name of the solid model. +============================ ============ ======== =========================================== + + diff --git a/src/coreComponents/schema/docs/PorousDamageElasticIsotropic_other.rst b/src/coreComponents/schema/docs/PorousDamageElasticIsotropic_other.rst new file mode 100644 index 00000000000..adf1c1b8aec --- /dev/null +++ b/src/coreComponents/schema/docs/PorousDamageElasticIsotropic_other.rst @@ -0,0 +1,9 @@ + + +==== ==== ============================ +Name Type Description +==== ==== ============================ + (no documentation available) +==== ==== ============================ + + diff --git a/src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic.rst b/src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic.rst new file mode 100644 index 00000000000..7d6d29602c1 --- /dev/null +++ b/src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic.rst @@ -0,0 +1,13 @@ + + +============================ ============ ======== =========================================== +Name Type Default Description +============================ ============ ======== =========================================== +name groupName required A name is required for any non-unique nodes +permeabilityModelName groupNameRef required Name of the permeability model. +porosityModelName groupNameRef required Name of the porosity model. +solidInternalEnergyModelName groupNameRef Name of the solid internal energy model. +solidModelName groupNameRef required Name of the solid model. +============================ ============ ======== =========================================== + + diff --git a/src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic_other.rst b/src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic_other.rst new file mode 100644 index 00000000000..adf1c1b8aec --- /dev/null +++ b/src/coreComponents/schema/docs/PorousDamageSpectralElasticIsotropic_other.rst @@ -0,0 +1,9 @@ + + +==== ==== ============================ +Name Type Description +==== ==== ============================ + (no documentation available) +==== ==== ============================ + + diff --git a/src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic.rst b/src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic.rst new file mode 100644 index 00000000000..7d6d29602c1 --- /dev/null +++ b/src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic.rst @@ -0,0 +1,13 @@ + + +============================ ============ ======== =========================================== +Name Type Default Description +============================ ============ ======== =========================================== +name groupName required A name is required for any non-unique nodes +permeabilityModelName groupNameRef required Name of the permeability model. +porosityModelName groupNameRef required Name of the porosity model. +solidInternalEnergyModelName groupNameRef Name of the solid internal energy model. +solidModelName groupNameRef required Name of the solid model. +============================ ============ ======== =========================================== + + diff --git a/src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic_other.rst b/src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic_other.rst new file mode 100644 index 00000000000..adf1c1b8aec --- /dev/null +++ b/src/coreComponents/schema/docs/PorousDamageVolDevElasticIsotropic_other.rst @@ -0,0 +1,9 @@ + + +==== ==== ============================ +Name Type Description +==== ==== ============================ + (no documentation available) +==== ==== ============================ + + diff --git a/src/coreComponents/schema/docs/VTKMesh.rst b/src/coreComponents/schema/docs/VTKMesh.rst index 724901b1747..688b0d98319 100644 --- a/src/coreComponents/schema/docs/VTKMesh.rst +++ b/src/coreComponents/schema/docs/VTKMesh.rst @@ -20,6 +20,7 @@ surfacicFieldsToImport groupNameRef_array {} Surfacic fields to be translate R1Tensor {0,0,0} Translate the coordinates of the vertices by a given vector (prior to scaling) useGlobalIds integer 0 Controls the use of global IDs in the input file for cells and points. If set to 0 (default value), the GlobalId arrays in the input mesh are used if available, and generated otherwise. If set to a negative value, the GlobalId arrays in the input mesh are not used, and generated global Ids are automatically generated. If set to a positive value, the GlobalId arrays in the input mesh are used and required, and the simulation aborts if they are not available InternalWell node :ref:`XML_InternalWell` +VTKWell node :ref:`XML_VTKWell` ====================== ======================== ========= ============================================================================================================================================================================================================================================================================================================================================================================================================================================================================ diff --git a/src/coreComponents/schema/docs/VTKMesh_other.rst b/src/coreComponents/schema/docs/VTKMesh_other.rst index aef746a9df7..8ca4f02bc36 100644 --- a/src/coreComponents/schema/docs/VTKMesh_other.rst +++ b/src/coreComponents/schema/docs/VTKMesh_other.rst @@ -4,6 +4,7 @@ Name Type Description ============ ==== ================================= InternalWell node :ref:`DATASTRUCTURE_InternalWell` +VTKWell node :ref:`DATASTRUCTURE_VTKWell` meshLevels node :ref:`DATASTRUCTURE_meshLevels` ============ ==== ================================= diff --git a/src/coreComponents/schema/docs/VTKWell.rst b/src/coreComponents/schema/docs/VTKWell.rst new file mode 100644 index 00000000000..0df45b1798d --- /dev/null +++ b/src/coreComponents/schema/docs/VTKWell.rst @@ -0,0 +1,17 @@ + + +===================== ========= ======== ==================================================================================================== +Name Type Default Description +===================== ========= ======== ==================================================================================================== +file path required Path to the well file +minElementLength real64 0.001 Minimum length of a well element, computed as (segment length / number of elements per segment ) [m] +minSegmentLength real64 0.01 Minimum length of a well segment [m] +name groupName required A name is required for any non-unique nodes +numElementsPerSegment integer required Number of well elements per polyline segment +radius real64 required Radius of the well [m] +wellControlsName string required Name of the set of constraints associated with this well +wellRegionName string required Name of the well element region +Perforation node :ref:`XML_Perforation` +===================== ========= ======== ==================================================================================================== + + diff --git a/src/coreComponents/schema/docs/VTKWell_other.rst b/src/coreComponents/schema/docs/VTKWell_other.rst new file mode 100644 index 00000000000..638cc779598 --- /dev/null +++ b/src/coreComponents/schema/docs/VTKWell_other.rst @@ -0,0 +1,9 @@ + + +=========== ==== ================================ +Name Type Description +=========== ==== ================================ +Perforation node :ref:`DATASTRUCTURE_Perforation` +=========== ==== ================================ + + diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 53106e73274..866cf51c2c9 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1438,12 +1438,20 @@ stress - traction is applied to the faces as specified by the inner product of i + + + + + + + + @@ -1451,6 +1459,10 @@ stress - traction is applied to the faces as specified by the inner product of i + + + + @@ -1462,6 +1474,12 @@ stress - traction is applied to the faces as specified by the inner product of i + + + + + + @@ -1511,9 +1529,9 @@ stress - traction is applied to the faces as specified by the inner product of i - + - + @@ -1527,6 +1545,27 @@ stress - traction is applied to the faces as specified by the inner product of i + + + + + + + + + + + + + + + + + + + + + @@ -1535,6 +1574,12 @@ stress - traction is applied to the faces as specified by the inner product of i + + + + + + @@ -1597,6 +1642,12 @@ stress - traction is applied to the faces as specified by the inner product of i + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index bfeab0520d8..8670ad568f9 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -380,6 +380,7 @@ + @@ -389,9 +390,15 @@ + + + + + + @@ -411,6 +418,7 @@ + diff --git a/src/coreComponents/unitTests/meshTests/well1.vtk b/src/coreComponents/unitTests/meshTests/well1.vtk new file mode 100644 index 00000000000..3a465cbbcfa --- /dev/null +++ b/src/coreComponents/unitTests/meshTests/well1.vtk @@ -0,0 +1,11 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET POLYDATA +POINTS 2 float +4.5 0 2 4.5 0 0.5 +LINES 2 2 +OFFSETS vtktypeint64 +0 2 +CONNECTIVITY vtktypeint64 +0 1 diff --git a/src/coreComponents/unitTests/meshTests/well2.vtk b/src/coreComponents/unitTests/meshTests/well2.vtk new file mode 100644 index 00000000000..2be27710ea2 --- /dev/null +++ b/src/coreComponents/unitTests/meshTests/well2.vtk @@ -0,0 +1,11 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET POLYDATA +POINTS 2 float +0.5 0 2 0.5 0 0.5 +LINES 2 2 +OFFSETS vtktypeint64 +0 2 +CONNECTIVITY vtktypeint64 +0 1 diff --git a/src/coreComponents/unitTests/wellsTests/testReservoirSinglePhaseMSWells.cpp b/src/coreComponents/unitTests/wellsTests/testReservoirSinglePhaseMSWells.cpp index dc1b2756e10..2948d71f72c 100644 --- a/src/coreComponents/unitTests/wellsTests/testReservoirSinglePhaseMSWells.cpp +++ b/src/coreComponents/unitTests/wellsTests/testReservoirSinglePhaseMSWells.cpp @@ -30,6 +30,8 @@ #include "physicsSolvers/fluidFlow/wells/SinglePhaseWellFields.hpp" #include "physicsSolvers/fluidFlow/wells/WellSolverBaseFields.hpp" +#include "tests/meshDirName.hpp" + using namespace geos; using namespace geos::dataRepository; using namespace geos::constitutive; @@ -37,7 +39,7 @@ using namespace geos::testing; CommandLineOptions g_commandLineOptions; -char const * xmlInput = +char const * PreXmlInput = R"xml( @@ -81,29 +83,10 @@ char const * xmlInput = nx="{3}" ny="{1}" nz="{1}" - cellBlockNames="{cb1}"> - - - - - - + cellBlockNames="{cb1}">)xml"; + +char const * PostXmlInput = + R"xml( @@ -349,6 +332,7 @@ void testNumericalJacobian( SinglePhaseReservoirAndWells< SinglePhaseBase > & so compareLocalMatrices( jacobian.toViewConst(), jacobianFD.toViewConst(), relTol ); } + class SinglePhaseReservoirSolverTest : public ::testing::Test { public: @@ -361,7 +345,6 @@ class SinglePhaseReservoirSolverTest : public ::testing::Test void SetUp() override { - setupProblemFromXML( state.getProblemManager(), xmlInput ); solver = &state.getProblemManager().getPhysicsSolverManager().getGroup< SinglePhaseReservoirAndWells< SinglePhaseBase > >( "reservoirSystem" ); DomainPartition & domain = state.getProblemManager().getDomainPartition(); @@ -372,39 +355,139 @@ class SinglePhaseReservoirSolverTest : public ::testing::Test solver->getSystemRhs(), solver->getSystemSolution() ); - solver->implicitStepSetup( time, dt, domain ); + solver->implicitStepSetup( TIME, DT, domain ); + } + + void TestAssembleCouplingTerms() + { + real64 const perturb = std::sqrt( EPS ); + real64 const tol = 1e-1; // 10% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testNumericalJacobian( *solver, domain, perturb, tol, + [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + solver->assembleCouplingTerms( TIME, DT, domain, solver->getDofManager(), localMatrix, localRhs ); + } ); + } + + void TestAssembleFluxTerms() + { + real64 const perturb = std::sqrt( EPS ); + real64 const tol = 1e-1; // 10% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testNumericalJacobian( *solver, domain, perturb, tol, + [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + solver->wellSolver()->assembleFluxTerms( DT, domain, solver->getDofManager(), localMatrix, localRhs ); + } ); + } + + void TestAssemblePressureRelations() + { + real64 const perturb = std::sqrt( EPS ); + real64 const tol = 1e-1; // 10% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testNumericalJacobian( *solver, domain, perturb, tol, + [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + solver->wellSolver()->assemblePressureRelations( TIME, DT, domain, solver->getDofManager(), localMatrix, localRhs ); + } ); + } + + void TestAssembleAccumulationTerms() + { + real64 const perturb = std::sqrt( EPS ); + real64 const tol = 1e-1; // 10% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testNumericalJacobian( *solver, domain, perturb, tol, + [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + solver->wellSolver()->assembleAccumulationTerms( domain, solver->getDofManager(), localMatrix, localRhs ); + } ); } - static real64 constexpr time = 0.0; - static real64 constexpr dt = 1e4; - static real64 constexpr eps = std::numeric_limits< real64 >::epsilon(); + static real64 constexpr TIME = 0.0; + static real64 constexpr DT = 1e4; + static real64 constexpr EPS = std::numeric_limits< real64 >::epsilon(); GeosxState state; SinglePhaseReservoirAndWells< SinglePhaseBase > * solver; }; -real64 constexpr SinglePhaseReservoirSolverTest::time; -real64 constexpr SinglePhaseReservoirSolverTest::dt; -real64 constexpr SinglePhaseReservoirSolverTest::eps; +real64 constexpr SinglePhaseReservoirSolverTest::TIME; +real64 constexpr SinglePhaseReservoirSolverTest::DT; +real64 constexpr SinglePhaseReservoirSolverTest::EPS; -TEST_F( SinglePhaseReservoirSolverTest, jacobianNumericalCheck_Perforation ) +/** + * @brief Test SinglePhaseReservoirSolver with InternalWell generator + * + */ +class SinglePhaseReservoirSolverInternalWellTest : public SinglePhaseReservoirSolverTest { - real64 const perturb = std::sqrt( eps ); - real64 const tol = 1e-1; // 10% error margin - DomainPartition & domain = state.getProblemManager().getDomainPartition(); +public: - testNumericalJacobian( *solver, domain, perturb, tol, - [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + SinglePhaseReservoirSolverInternalWellTest(): + SinglePhaseReservoirSolverTest() + {} + +protected: + + void SetUp() override { - solver->assembleCouplingTerms( time, dt, domain, solver->getDofManager(), localMatrix, localRhs ); - } ); + string const internalWells = + R"( + + + + + + )"; + + string const xmlInput = PreXmlInput + internalWells + PostXmlInput; + + setupProblemFromXML( state.getProblemManager(), xmlInput.c_str() ); + SinglePhaseReservoirSolverTest::SetUp(); + } +}; + + +TEST_F( SinglePhaseReservoirSolverInternalWellTest, jacobianNumericalCheck_Perforation ) +{ + TestAssembleCouplingTerms(); } -TEST_F( SinglePhaseReservoirSolverTest, jacobianNumericalCheck_Flux ) +TEST_F( SinglePhaseReservoirSolverInternalWellTest, jacobianNumericalCheck_Flux ) { - real64 const perturb = std::sqrt( eps ); + real64 const perturb = std::sqrt( EPS ); real64 const tol = 1e-1; // 10% error margin DomainPartition & domain = state.getProblemManager().getDomainPartition(); @@ -413,40 +496,86 @@ TEST_F( SinglePhaseReservoirSolverTest, jacobianNumericalCheck_Flux ) [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - solver->wellSolver()->assembleFluxTerms( dt, domain, solver->getDofManager(), localMatrix, localRhs ); + solver->wellSolver()->assembleFluxTerms( DT, domain, solver->getDofManager(), localMatrix, localRhs ); } ); } -TEST_F( SinglePhaseReservoirSolverTest, jacobianNumericalCheck_PressureRel ) +TEST_F( SinglePhaseReservoirSolverInternalWellTest, jacobianNumericalCheck_PressureRel ) { - real64 const perturb = std::sqrt( eps ); - real64 const tol = 1e-1; // 10% error margin - - DomainPartition & domain = state.getProblemManager().getDomainPartition(); + TestAssemblePressureRelations(); +} - testNumericalJacobian( *solver, domain, perturb, tol, - [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) - { - solver->wellSolver()->assemblePressureRelations( time, dt, domain, solver->getDofManager(), localMatrix, localRhs ); - } ); +TEST_F( SinglePhaseReservoirSolverInternalWellTest, jacobianNumericalCheck_Accum ) +{ + TestAssembleAccumulationTerms(); } -TEST_F( SinglePhaseReservoirSolverTest, jacobianNumericalCheck_Accum ) + +/** + * @brief Test SinglePhaseReservoirSolver with VTKWell generator + * + */ +class SinglePhaseReservoirSolverVTKWellTest : public SinglePhaseReservoirSolverTest { - real64 const perturb = std::sqrt( eps ); - real64 const tol = 1e-1; // 10% error margin - DomainPartition & domain = state.getProblemManager().getDomainPartition(); +public: - testNumericalJacobian( *solver, domain, perturb, tol, - [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) + SinglePhaseReservoirSolverVTKWellTest(): + SinglePhaseReservoirSolverTest() + {} + +protected: + + void SetUp() override { - solver->wellSolver()->assembleAccumulationTerms( domain, solver->getDofManager(), localMatrix, localRhs ); - } ); + string const vtkWells = GEOS_FMT( R"( + + + + + + )", testMeshDir + "/well1.vtk", + testMeshDir + "/well2.vtk" ); + + string const xmlInput = PreXmlInput + vtkWells + PostXmlInput; + + setupProblemFromXML( state.getProblemManager(), xmlInput.c_str()); + SinglePhaseReservoirSolverTest::SetUp(); + } +}; + + +TEST_F( SinglePhaseReservoirSolverVTKWellTest, jacobianNumericalCheck_Perforation ) +{ + TestAssembleCouplingTerms(); } +TEST_F( SinglePhaseReservoirSolverVTKWellTest, jacobianNumericalCheck_Flux ) +{ + TestAssembleFluxTerms(); +} + +TEST_F( SinglePhaseReservoirSolverVTKWellTest, jacobianNumericalCheck_PressureRel ) +{ + TestAssemblePressureRelations(); +} + +TEST_F( SinglePhaseReservoirSolverVTKWellTest, jacobianNumericalCheck_Accum ) +{ + TestAssembleAccumulationTerms(); +} int main( int argc, char * * argv ) { diff --git a/src/docs/sphinx/CompleteXMLSchema.rst b/src/docs/sphinx/CompleteXMLSchema.rst index 774ada1b92e..751de3551d6 100644 --- a/src/docs/sphinx/CompleteXMLSchema.rst +++ b/src/docs/sphinx/CompleteXMLSchema.rst @@ -1319,6 +1319,13 @@ Element: VTKMesh .. include:: ../../coreComponents/schema/docs/VTKMesh.rst +.. _XML_VTKWell: + +Element: VTKWell +================ +.. include:: ../../coreComponents/schema/docs/VTKWell.rst + + .. _XML_VanGenuchtenBakerRelativePermeability: Element: VanGenuchtenBakerRelativePermeability @@ -2744,6 +2751,13 @@ Datastructure: VTKMesh .. include:: ../../coreComponents/schema/docs/VTKMesh_other.rst +.. _DATASTRUCTURE_VTKWell: + +Datastructure: VTKWell +====================== +.. include:: ../../coreComponents/schema/docs/VTKWell_other.rst + + .. _DATASTRUCTURE_VanGenuchtenBakerRelativePermeability: Datastructure: VanGenuchtenBakerRelativePermeability From dbd1eeb9e0c93e555d500a0526da9e47638fbb93 Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" <50467563+victorapm@users.noreply.github.com> Date: Sun, 28 Jan 2024 12:41:30 -0500 Subject: [PATCH 3/5] Improve `SinglePhasePoromechanicsEmbeddedFractures` and introduce `SinglePhasePoromechanicsConformingFractures` MGR strategy (#2747) --- ...conformingFracture_2d_openingFrac_base.xml | 1 + .../linearAlgebra/CMakeLists.txt | 1 + .../interfaces/hypre/HypreMGR.cpp | 6 + .../interfaces/hypre/HypreSolver.cpp | 2 +- ...ePhasePoromechanicsConformingFractures.hpp | 137 ++++++++++++++++++ ...glePhasePoromechanicsEmbeddedFractures.hpp | 79 +++++----- .../utilities/LinearSolverParameters.hpp | 44 +++--- ...ePhasePoromechanicsConformingFractures.cpp | 5 + 8 files changed, 209 insertions(+), 66 deletions(-) create mode 100644 src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp diff --git a/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml b/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml index 18d744f1300..543423eac0f 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml @@ -15,6 +15,7 @@ newtonMaxIter="10" maxNumConfigurationAttempts="10"/> diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt index 1a1bc7ca108..3ae8298a43d 100644 --- a/src/coreComponents/linearAlgebra/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/CMakeLists.txt @@ -103,6 +103,7 @@ if( ENABLE_HYPRE ) interfaces/hypre/mgrStrategies/LagrangianContactMechanics.hpp interfaces/hypre/mgrStrategies/MultiphasePoromechanics.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp + interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp index b41b7e9577a..37f471729c5 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp @@ -32,6 +32,7 @@ #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseHybridFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanics.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp" +#include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhaseReservoirHybridFVM.hpp" @@ -154,6 +155,11 @@ void hypre::mgr::createMGR( LinearSolverParameters const & params, setStrategy< SinglePhasePoromechanicsEmbeddedFractures >( params.mgr, numComponentsPerField, precond, mgrData ); break; } + case LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFractures: + { + setStrategy< SinglePhasePoromechanicsConformingFractures >( params.mgr, numComponentsPerField, precond, mgrData ); + break; + } case LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsReservoirFVM: { setStrategy< SinglePhasePoromechanicsReservoirFVM >( params.mgr, numComponentsPerField, precond, mgrData ); diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp index 04516475be2..3a54f603c3c 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreSolver.cpp @@ -77,7 +77,7 @@ void createHypreGMRES( LinearSolverParameters const & params, GEOS_LAI_CHECK_ERROR( HYPRE_ParCSRGMRESSetTol( solver.ptr, params.krylov.relTolerance ) ); // Default for now - HYPRE_Int logLevel = (params.logLevel >= 3) ? 1 : 0; + HYPRE_Int logLevel = (params.logLevel >= 3) ? 2 : 0; GEOS_LAI_CHECK_ERROR( HYPRE_ParCSRGMRESSetPrintLevel( solver.ptr, logLevel ) ); // print iteration info GEOS_LAI_CHECK_ERROR( HYPRE_ParCSRGMRESSetLogging( solver.ptr, 1 ) ); /* needed to get run info later */ diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp new file mode 100644 index 00000000000..0dfb0b0ec88 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp @@ -0,0 +1,137 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SinglePhasePoromechanicsConformingFractures.hpp + */ + +#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ +#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ + +#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp" + +#define BRANCH_MGR_FSOLVER + +namespace geos +{ + +namespace hypre +{ + +namespace mgr +{ + +/** + * @brief SinglePhasePoromechanicsConformingFractures strategy. + * + * dofLabel: 0 = displacement, x-component + * dofLabel: 1 = displacement, y-component + * dofLabel: 2 = displacement, z-component + * dofLabel: 3 = pressure (cell elem + fracture elems) + * dofLabel: 4 = face-centered lagrange multiplier (tn) + * dofLabel: 5 = face-centered lagrange multiplier (tt1) + * dofLabel: 6 = face-centered lagrange multiplier (tt2) + + * + * Ingredients: + * 1. Level 1: F-points displacement (4,5,6), C-points pressure (0,1,2,3) + * 2. Level 2: F-points displacement (0,1,2), C-points pressure (3) + * 2. F-points smoother: BoomerAMG, single V-cycle + * 3. C-points coarse-grid/Schur complement solver: BoomerAMG + * 4. Global smoother: none + */ +class SinglePhasePoromechanicsConformingFractures : public MGRStrategyBase< 2 > +{ +public: + + /** + * @brief Constructor. + */ + explicit SinglePhasePoromechanicsConformingFractures( arrayView1d< int const > const & ) + : MGRStrategyBase( 7 ) + { + + // we keep u and p + m_labels[0].push_back( 0 ); + m_labels[0].push_back( 1 ); + m_labels[0].push_back( 2 ); + m_labels[0].push_back( 3 ); + // we keep p + m_labels[1].push_back( 3 ); + + setupLabels(); + + // Level 0 + m_levelFRelaxType[0] = MGRFRelaxationType::none; + m_levelFRelaxIters[0] = 0; + + m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::ilu0; + m_levelGlobalSmootherIters[0] = 1; + + m_levelInterpType[0] = MGRInterpolationType::blockJacobi; + m_levelRestrictType[0] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; + + // Level 1 + m_levelFRelaxType[1] = MGRFRelaxationType::amgVCycle; + m_levelFRelaxIters[1] = 1; + m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::none; + m_levelInterpType[1] = MGRInterpolationType::jacobi; + m_levelRestrictType[1] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; + + } + + /** + * @brief Setup the MGR strategy. + * @param precond preconditioner wrapper + * @param mgrData auxiliary MGR data + */ + void setup( LinearSolverParameters::MGR const &, + HyprePrecWrapper & precond, + HypreMGRData & mgrData ) + { + setReduction( precond, mgrData ); + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetPMaxElmts( precond.ptr, 0 )); + + // Configure the BoomerAMG solver used as F-relaxation for the second level + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); + +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( mgrData.mechSolver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( mgrData.mechSolver.ptr, hypre::getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::chebyshev ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( mgrData.mechSolver.ptr, 1 ) ); +#else + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( mgrData.mechSolver.ptr, 1 ) ); +#endif + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFSolverAtLevel( 1, precond.ptr, mgrData.mechSolver.ptr ) ); + + // Configure the BoomerAMG solver used as mgr coarse solver for the pressure reduced system + setPressureAMG( mgrData.coarseSolver ); + } +}; + +} // namespace mgr + +} // namespace hypre + +} // namespace geos + +#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_*/ diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp index a4ebe344685..36d8051f1ad 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -42,10 +42,11 @@ namespace mgr * dofLabel: 6 = pressure (cell elem + fracture elems) * * Ingredients: - * 1. F-points displacement (0,1,2), C-points pressure (3) - * 2. F-points smoother: AMG, single V-cycle, separate displacement components - * 3. C-points coarse-grid/Schur complement solver: boomer AMG - * 4. Global smoother: none + * 1. Level 1: F-points w (3,4,5), C-points pressure (0,1,2,6) + * 2. Level 2: F-points displacement (0,1,2), C-points pressure (6) + * 3. F-points smoother: BoomerAMG, single V-cycle + * 4. C-points coarse-grid/Schur complement solver: BoomerAMG + * 5. Global smoother: none */ class SinglePhasePoromechanicsEmbeddedFractures : public MGRStrategyBase< 2 > { @@ -57,35 +58,10 @@ class SinglePhasePoromechanicsEmbeddedFractures : public MGRStrategyBase< 2 > explicit SinglePhasePoromechanicsEmbeddedFractures( arrayView1d< int const > const & ) : MGRStrategyBase( 7 ) { - - /// IDEAL strategy but it seg faults - // // we keep u and p - // m_labels[0].push_back( 0 ); - // m_labels[0].push_back( 1 ); - // m_labels[0].push_back( 2 ); - // m_labels[0].push_back( 6 ); - // // we keep p - // m_labels[1].push_back( 6 ); - - // setupLabels(); - - // // Level 0 - // m_levelFRelaxMethod[0] = MGRFRelaxationMethod::singleLevel; - // m_levelFRelaxType[0] = MGRFRelaxationType::gsElimWInverse; // gaussian elimination for the dispJump block - // m_levelInterpType[0] = MGRInterpolationType::blockJacobi; - // m_levelRestrictType[0] = MGRRestrictionType::injection; - // m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; - - // // Level 1 - // m_levelFRelaxMethod[1] = MGRFRelaxationMethod::amgVCycle; - // m_levelInterpType[1] = MGRInterpolationType::jacobi; - // m_levelRestrictType[1] = MGRRestrictionType::injection; - // m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; - - // we keep w and p - m_labels[0].push_back( 3 ); - m_labels[0].push_back( 4 ); - m_labels[0].push_back( 5 ); + // we keep u and p + m_labels[0].push_back( 0 ); + m_labels[0].push_back( 1 ); + m_labels[0].push_back( 2 ); m_labels[0].push_back( 6 ); // we keep p m_labels[1].push_back( 6 ); @@ -93,20 +69,20 @@ class SinglePhasePoromechanicsEmbeddedFractures : public MGRStrategyBase< 2 > setupLabels(); // Level 0 - m_levelFRelaxType[0] = MGRFRelaxationType::amgVCycle; + m_levelFRelaxType[0] = MGRFRelaxationType::gsElimWPivoting; // gaussian elimination for the dispJump block m_levelFRelaxIters[0] = 1; - m_levelInterpType[0] = MGRInterpolationType::jacobi; - m_levelRestrictType[0] = MGRRestrictionType::injection; - m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::nonGalerkin; m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none; + m_levelInterpType[0] = MGRInterpolationType::blockJacobi; + m_levelRestrictType[0] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; // Level 1 - m_levelFRelaxType[1] = MGRFRelaxationType::gsElimWInverse; - m_levelFRelaxIters[1] = 1; - m_levelInterpType[1] = MGRInterpolationType::blockJacobi; - m_levelRestrictType[1] = MGRRestrictionType::injection; - m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::galerkin; + m_levelFRelaxType[1] = MGRFRelaxationType::amgVCycle; m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::none; + m_levelInterpType[1] = MGRInterpolationType::jacobi; + m_levelRestrictType[1] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::nonGalerkin; + } /** @@ -122,8 +98,23 @@ class SinglePhasePoromechanicsEmbeddedFractures : public MGRStrategyBase< 2 > GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetPMaxElmts( precond.ptr, 0 )); - // Configure the BoomerAMG solver used as F-relaxation for the first level - setMechanicsFSolver( precond, mgrData ); + // Configure the BoomerAMG solver used as F-relaxation for the second level + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGCreate( &mgrData.mechSolver.ptr ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetTol( mgrData.mechSolver.ptr, 0.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxIter( mgrData.mechSolver.ptr, 1 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetMaxRowSum( mgrData.mechSolver.ptr, 1.0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetStrongThreshold( mgrData.mechSolver.ptr, 0.6 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetPrintLevel( mgrData.mechSolver.ptr, 0 ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumFunctions( mgrData.mechSolver.ptr, 3 ) ); + +#if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetCoarsenType( mgrData.mechSolver.ptr, hypre::getAMGCoarseningType( LinearSolverParameters::AMG::CoarseningType::PMIS ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxType( mgrData.mechSolver.ptr, hypre::getAMGRelaxationType( LinearSolverParameters::AMG::SmootherType::chebyshev ) ) ); + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetNumSweeps( mgrData.mechSolver.ptr, 1 ) ); +#else + GEOS_LAI_CHECK_ERROR( HYPRE_BoomerAMGSetRelaxOrder( mgrData.mechSolver.ptr, 1 ) ); +#endif + GEOS_LAI_CHECK_ERROR( HYPRE_MGRSetFSolverAtLevel( 1, precond.ptr, mgrData.mechSolver.ptr ) ); // Configure the BoomerAMG solver used as mgr coarse solver for the pressure reduced system setPressureAMG( mgrData.coarseSolver ); diff --git a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp index 46c6280262d..881fcc9984c 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -262,27 +262,28 @@ struct LinearSolverParameters */ enum class StrategyType : integer { - invalid, ///< default value, to ensure solver sets something - singlePhaseReservoirFVM, ///< finite volume single-phase flow with wells - singlePhaseHybridFVM, ///< hybrid finite volume single-phase flow - singlePhaseReservoirHybridFVM, ///< hybrid finite volume single-phase flow with wells - singlePhasePoromechanics, ///< single phase poromechanics with finite volume single phase flow - thermalSinglePhasePoromechanics, ///< thermal single phase poromechanics with finite volume single phase flow - hybridSinglePhasePoromechanics, ///< single phase poromechanics with hybrid finite volume single phase flow - singlePhasePoromechanicsEmbeddedFractures, ///< single phase poromechanics with finite volume single phase flow and embedded fractures - singlePhasePoromechanicsReservoirFVM, ///< single phase poromechanics with finite volume single phase flow with wells - compositionalMultiphaseFVM, ///< finite volume compositional multiphase flow - compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow - compositionalMultiphaseReservoirFVM, ///< finite volume compositional multiphase flow with wells - compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells - reactiveCompositionalMultiphaseOBL, ///< finite volume reactive compositional flow with OBL - thermalCompositionalMultiphaseFVM, ///< finite volume thermal compositional multiphase flow - multiphasePoromechanics, ///< multiphase poromechanics with finite volume compositional multiphase flow - multiphasePoromechanicsReservoirFVM, ///< multiphase poromechanics with finite volume compositional multiphase flow with wells - thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow - hydrofracture, ///< hydrofracture - lagrangianContactMechanics, ///< Lagrangian contact mechanics - solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics + invalid, ///< default value, to ensure solver sets something + singlePhaseReservoirFVM, ///< finite volume single-phase flow with wells + singlePhaseHybridFVM, ///< hybrid finite volume single-phase flow + singlePhaseReservoirHybridFVM, ///< hybrid finite volume single-phase flow with wells + singlePhasePoromechanics, ///< single phase poromechanics with finite volume single phase flow + thermalSinglePhasePoromechanics, ///< thermal single phase poromechanics with finite volume single phase flow + hybridSinglePhasePoromechanics, ///< single phase poromechanics with hybrid finite volume single phase flow + singlePhasePoromechanicsEmbeddedFractures, ///< single phase poromechanics with FV embedded fractures + singlePhasePoromechanicsConformingFractures, ///< single phase poromechanics with FV conforming fractures + singlePhasePoromechanicsReservoirFVM, ///< single phase poromechanics with finite volume single phase flow with wells + compositionalMultiphaseFVM, ///< finite volume compositional multiphase flow + compositionalMultiphaseHybridFVM, ///< hybrid finite volume compositional multiphase flow + compositionalMultiphaseReservoirFVM, ///< finite volume compositional multiphase flow with wells + compositionalMultiphaseReservoirHybridFVM, ///< hybrid finite volume compositional multiphase flow with wells + reactiveCompositionalMultiphaseOBL, ///< finite volume reactive compositional flow with OBL + thermalCompositionalMultiphaseFVM, ///< finite volume thermal compositional multiphase flow + multiphasePoromechanics, ///< multiphase poromechanics with finite volume compositional multiphase flow + multiphasePoromechanicsReservoirFVM, ///< multiphase poromechanics with finite volume compositional multiphase flow with wells + thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow + hydrofracture, ///< hydrofracture + lagrangianContactMechanics, ///< Lagrangian contact mechanics + solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics }; StrategyType strategy = StrategyType::invalid; ///< Predefined MGR solution strategy (solver specific) @@ -361,6 +362,7 @@ ENUM_STRINGS( LinearSolverParameters::MGR::StrategyType, "thermalSinglePhasePoromechanics", "hybridSinglePhasePoromechanics", "singlePhasePoromechanicsEmbeddedFractures", + "singlePhasePoromechanicsConformingFractures", "singlePhasePoromechanicsReservoirFVM", "compositionalMultiphaseFVM", "compositionalMultiphaseHybridFVM", diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp index cce65a5fd42..93a32ddb123 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp @@ -45,6 +45,11 @@ SinglePhasePoromechanicsConformingFractures::SinglePhasePoromechanicsConformingF setApplyDefaultValue( 0 ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "Flag indicating whether the problem is thermal or not. Set isThermal=\"1\" to enable the thermal coupling" ); + + m_linearSolverParameters.get().mgr.strategy = LinearSolverParameters::MGR::StrategyType::singlePhasePoromechanicsConformingFractures; + m_linearSolverParameters.get().mgr.separateComponents = false; + m_linearSolverParameters.get().mgr.displacementFieldName = solidMechanics::totalDisplacement::key(); + m_linearSolverParameters.get().dofsPerNode = 3; } void SinglePhasePoromechanicsConformingFractures::initializePostInitialConditionsPostSubGroups() From 125caabe042c6be489a2e1b5313ca2f08cb55bc6 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 29 Jan 2024 09:18:45 -0600 Subject: [PATCH 4/5] Use output dir (option "-o") for printing out PVT tables (#2946) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Pavel Tomin <“paveltomin@users.noreply.github.com”> --- src/coreComponents/functions/CMakeLists.txt | 2 +- src/coreComponents/functions/TableFunction.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/functions/CMakeLists.txt b/src/coreComponents/functions/CMakeLists.txt index f28de10f61f..01f8d6fed86 100644 --- a/src/coreComponents/functions/CMakeLists.txt +++ b/src/coreComponents/functions/CMakeLists.txt @@ -27,7 +27,7 @@ if( ENABLE_MATHPRESSO ) endif() -set( dependencyList ${parallelDeps} codingUtilities dataRepository ) +set( dependencyList ${parallelDeps} codingUtilities dataRepository fileIO ) if( ENABLE_MATHPRESSO ) list( APPEND dependencyList mathpresso ) diff --git a/src/coreComponents/functions/TableFunction.cpp b/src/coreComponents/functions/TableFunction.cpp index 892b2289627..0a8751c28ab 100644 --- a/src/coreComponents/functions/TableFunction.cpp +++ b/src/coreComponents/functions/TableFunction.cpp @@ -19,6 +19,7 @@ #include "TableFunction.hpp" #include "codingUtilities/Parsing.hpp" #include "common/DataTypes.hpp" +#include "fileIO/Outputs/OutputBase.hpp" #include @@ -182,7 +183,7 @@ void TableFunction::checkCoord( real64 const coord, localIndex const dim ) const void TableFunction::print( std::string const & filename ) const { - std::ofstream os( filename + ".csv" ); + std::ofstream os( joinPath( OutputBase::getOutputDirectory(), filename + ".csv" ) ); integer const numDimensions = LvArray::integerConversion< integer >( m_coordinates.size() ); From 3cdbde87c3ca7d356e4d70ad47faa3e383ed872f Mon Sep 17 00:00:00 2001 From: TotoGaz <49004943+TotoGaz@users.noreply.github.com> Date: Mon, 29 Jan 2024 10:47:43 -0800 Subject: [PATCH 5/5] Update LvArray (#2956) Add new pretty printers. --- src/coreComponents/LvArray | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/LvArray b/src/coreComponents/LvArray index 6057f598b00..6cb244ecf76 160000 --- a/src/coreComponents/LvArray +++ b/src/coreComponents/LvArray @@ -1 +1 @@ -Subproject commit 6057f598b005db6efd132f60c298ad9c827fe3cd +Subproject commit 6cb244ecf76810fe738dca0a9201ca539533a343