diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 5076ede8bd3..8e8f20b5ff1 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr4040-16993-1393f80 + baseline: integratedTests/baseline_integratedTests-pr4077-17009-64cf777 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 024e6787153..9de6e3cc27c 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -5,6 +5,9 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #4077 (2026-06-19) +Add input mesh parameter to output Pvd results + PR #4040 (2026-06-16) Move relperm driver to use new constitutive driver framework diff --git a/src/coreComponents/mesh/CellElementSubRegion.cpp b/src/coreComponents/mesh/CellElementSubRegion.cpp index 64e82d1da65..7bf60d54a40 100644 --- a/src/coreComponents/mesh/CellElementSubRegion.cpp +++ b/src/coreComponents/mesh/CellElementSubRegion.cpp @@ -91,7 +91,9 @@ void CellElementSubRegion::copyFromCellBlock( CellBlockABC const & cellBlock ) { using ArrayType = camp::first< decltype( tupleOfTypes ) >; auto const src = Wrapper< ArrayType >::cast( wrapper ).reference().toViewConst(); - ArrayType & dst = this->registerWrapper( wrapper.getName(), std::make_unique< ArrayType >() ).reference(); + ArrayType & dst = this->registerWrapper( wrapper.getName(), std::make_unique< ArrayType >() ) + .setPlotLevel( wrapper.getPlotLevel() ) + .reference(); // This is a hack since Array's copy ctor does not accept ArrayView source dst.resize( ArrayType::NDIM, src.dims() ); std::copy( src.data(), src.data() + src.size(), dst.data() ); diff --git a/src/coreComponents/mesh/ElementSubRegionBase.cpp b/src/coreComponents/mesh/ElementSubRegionBase.cpp index 6234dbaca94..70662168d50 100644 --- a/src/coreComponents/mesh/ElementSubRegionBase.cpp +++ b/src/coreComponents/mesh/ElementSubRegionBase.cpp @@ -47,6 +47,7 @@ ElementSubRegionBase::ElementSubRegionBase( string const & name, Group * const p registerWrapper( viewKeyStruct::elementVolumeString(), &m_elementVolume ). setPlotLevel( PlotLevel::LEVEL_1 ); + } ElementSubRegionBase::~ElementSubRegionBase() diff --git a/src/coreComponents/mesh/FaceElementSubRegion.cpp b/src/coreComponents/mesh/FaceElementSubRegion.cpp index 877a559633b..09990e0f55f 100644 --- a/src/coreComponents/mesh/FaceElementSubRegion.cpp +++ b/src/coreComponents/mesh/FaceElementSubRegion.cpp @@ -15,6 +15,7 @@ #include "FaceElementSubRegion.hpp" #include "common/GEOS_RAJA_Interface.hpp" +#include "common/TypeDispatch.hpp" #include "NodeManager.hpp" #include "MeshLevel.hpp" @@ -214,7 +215,20 @@ void FaceElementSubRegion::copyFromCellBlock( FaceBlockABC const & faceBlock ) m_newFaceElements.insert( i ); } - // TODO We still need to be able to import fields on the FaceElementSubRegion. + faceBlock.forExternalProperties( [&]( WrapperBase const & wrapper ) + { + types::dispatch( types::ListofTypeList< types::StandardArrays >{}, [&]( auto tupleOfTypes ) + { + using ArrayType = camp::first< decltype( tupleOfTypes ) >; + auto const src = Wrapper< ArrayType >::cast( wrapper ).reference().toViewConst(); + ArrayType & dst = this->registerWrapper( wrapper.getName(), std::make_unique< ArrayType >() ) + .setPlotLevel( wrapper.getPlotLevel() ) + .reference(); + // This is a hack since Array's copy ctor does not accept ArrayView source + dst.resize( ArrayType::NDIM, src.dims() ); + std::copy( src.data(), src.data() + src.size(), dst.data() ); + }, wrapper ); + } ); } void FaceElementSubRegion::setupRelatedObjectsInRelations( MeshLevel const & mesh ) diff --git a/src/coreComponents/mesh/generators/FaceBlock.hpp b/src/coreComponents/mesh/generators/FaceBlock.hpp index 55b0d5840b8..31060e8e03b 100644 --- a/src/coreComponents/mesh/generators/FaceBlock.hpp +++ b/src/coreComponents/mesh/generators/FaceBlock.hpp @@ -120,11 +120,37 @@ class FaceBlock : public FaceBlockABC */ void set2dElemsToCollocatedNodesBuckets( ArrayOfArrays< array1d< globalIndex > > && collocatedNodesBuckets ); + /** + * @brief Add a property to the FaceBlock. + * @tparam T type of the property + * @param[in] propertyName the name of the property + * @return a non-const reference to the property + */ + template< typename T > + T & addProperty( string const & propertyName ) + { + m_externalPropertyNames.emplace_back( propertyName ); + return this->registerWrapper< T >( propertyName ).reference(); + } + private: + std::list< dataRepository::WrapperBase const * > getExternalProperties() const override + { + std::list< dataRepository::WrapperBase const * > result; + for( string const & externalPropertyName : m_externalPropertyNames ) + { + result.push_back( &this->getWrapperBase( externalPropertyName ) ); + } + return result; + } + localIndex m_num2dElements; localIndex m_num2dFaces; + /// Names of the properties registered from an external mesh + string_array m_externalPropertyNames; + ArrayOfArrays< localIndex > m_2dElemToNodes; ArrayOfArrays< localIndex > m_2dElemToEdges; ArrayOfArrays< localIndex > m_2dElemToFaces; diff --git a/src/coreComponents/mesh/generators/FaceBlockABC.hpp b/src/coreComponents/mesh/generators/FaceBlockABC.hpp index 2df82db9913..3708622349c 100644 --- a/src/coreComponents/mesh/generators/FaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/FaceBlockABC.hpp @@ -20,6 +20,8 @@ #include "dataRepository/Group.hpp" #include "common/DataTypes.hpp" +#include + namespace geos { @@ -170,6 +172,31 @@ class FaceBlockABC : public dataRepository::Group * @return The mapping relationship as an array. */ virtual array1d< globalIndex > localToGlobalMap() const = 0; + + /** + * @brief Apply a lambda to each external (passthrough) property registered on this face block. + * @tparam LAMBDA type of the lambda + * @param lambda the lambda to apply; receives a @p dataRepository::WrapperBase const & + */ + template< typename LAMBDA > + void forExternalProperties( LAMBDA && lambda ) const + { + for( auto * wrapperBase : this->getExternalProperties() ) + { + lambda( *wrapperBase ); + } + } + +private: + + /** + * @brief Returns the list of external (passthrough) property wrappers. + * @return List of pointers to external property wrappers. + */ + virtual std::list< dataRepository::WrapperBase const * > getExternalProperties() const + { + return {}; + } }; } diff --git a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp index 2dd962a156f..c13c66408a6 100644 --- a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp @@ -18,6 +18,7 @@ #include "mesh/generators/VTKUtilities.hpp" #include "mesh/generators/CollocatedNodes.hpp" +#include "common/MpiWrapper.hpp" #include "dataRepository/Group.hpp" #include @@ -624,10 +625,116 @@ array1d< globalIndex > buildLocalToGlobal( vtkIdTypeArray const * faceMeshCellGl } +/** + * @brief Register a passthrough field on a FaceBlock, ensuring MPI consistency. + * + * On ranks where the VTK array is absent (e.g. empty fracture partition), an empty + * wrapper of the correct type is still registered so all ranks have the same wrappers + * during ghost communication. + */ +static void writeFacePassthroughField( vtkDataSet & faceMesh, + localIndex const numElems, + string const & fieldName, + FaceBlock & faceBlock ) +{ + vtkDataArray * srcArray = faceMesh.GetCellData()->GetArray( fieldName.c_str() ); + + // Collectively determine type: isReal (1/0) and numComp, or -1 if absent on this rank. + int localIsReal = ( srcArray != nullptr ) + ? ( ( srcArray->GetDataType() == VTK_FLOAT || srcArray->GetDataType() == VTK_DOUBLE ) ? 1 : 0 ) + : -1; + int localNumComp = ( srcArray != nullptr ) ? srcArray->GetNumberOfComponents() : -1; + + int const globalIsReal = MpiWrapper::allReduce( localIsReal, MpiWrapper::Reduction::Max, MPI_COMM_GEOS ); + int const globalNumComp = MpiWrapper::allReduce( localNumComp, MpiWrapper::Reduction::Max, MPI_COMM_GEOS ); + + if( globalNumComp == -1 ) + { + // Field absent on every rank — not a passthrough field of this mesh. + GEOS_WARNING( GEOS_FMT( "Passthrough field '{}' not found in any fracture VTK CellData; skipping.", fieldName ) ); + return; + } + + if( srcArray == nullptr ) + { + // This rank has no data for this field (empty partition). Register an empty + // wrapper so ghost communication finds the same wrapper on all ranks. + if( globalIsReal == 1 ) + { + if( globalNumComp == 1 ) + { + faceBlock.addProperty< array1d< real64 > >( fieldName ).resize( numElems ); + } + else + { + faceBlock.addProperty< array2d< real64 > >( fieldName ).resize( numElems, globalNumComp ); + } + } + else + { + if( globalNumComp == 1 ) + { + faceBlock.addProperty< array1d< integer > >( fieldName ).resize( numElems ); + } + else + { + faceBlock.addProperty< array2d< integer > >( fieldName ).resize( numElems, globalNumComp ); + } + } + faceBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + return; + } + + // This rank has the data — fill it. + bool const isReal = ( globalIsReal == 1 ); + if( isReal ) + { + if( globalNumComp == 1 ) + { + array1d< real64 > & dst = faceBlock.addProperty< array1d< real64 > >( fieldName ); + dst.resize( numElems ); + faceBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + for( localIndex i = 0; i < numElems; ++i ) + dst[i] = static_cast< real64 >( srcArray->GetTuple1( i ) ); + } + else + { + array2d< real64 > & dst = faceBlock.addProperty< array2d< real64 > >( fieldName ); + dst.resize( numElems, globalNumComp ); + faceBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + for( localIndex i = 0; i < numElems; ++i ) + for( int c = 0; c < globalNumComp; ++c ) + dst( i, c ) = static_cast< real64 >( srcArray->GetComponent( i, c ) ); + } + } + else + { + if( globalNumComp == 1 ) + { + array1d< integer > & dst = faceBlock.addProperty< array1d< integer > >( fieldName ); + dst.resize( numElems ); + faceBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + for( localIndex i = 0; i < numElems; ++i ) + dst[i] = static_cast< integer >( srcArray->GetTuple1( i ) ); + } + else + { + array2d< integer > & dst = faceBlock.addProperty< array2d< integer > >( fieldName ); + dst.resize( numElems, globalNumComp ); + faceBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + for( localIndex i = 0; i < numElems; ++i ) + for( int c = 0; c < globalNumComp; ++c ) + dst( i, c ) = static_cast< integer >( srcArray->GetComponent( i, c ) ); + } + } +} + + void importFractureNetwork( string const & faceBlockName, vtkSmartPointer< vtkDataSet > faceMesh, vtkSmartPointer< vtkDataSet > mesh, - CellBlockManager & cellBlockManager ) + CellBlockManager & cellBlockManager, + Span< string const > passthroughFieldNames ) { ArrayOfArrays< localIndex > const faceToNodes = cellBlockManager.getFaceToNodes(); vtk::internal::ElementToFace const elemToFaces( cellBlockManager.getCellBlocks() ); @@ -686,6 +793,16 @@ void importFractureNetwork( string const & faceBlockName, buildCollocatedNodesBucketsOf2dElemsMap( build2dElemTo2dNodes( faceMesh ), buildCollocatedNodesMap( collocatedNodes ) ) ); + + // Import passthrough fields from the face mesh CellData. + // writeFacePassthroughField performs an MPI_Allreduce to determine the field type + // collectively, ensuring all ranks register the same wrappers even when some ranks + // have an empty fracture partition (needed for ghost communication consistency). + localIndex const numElems = LvArray::integerConversion< localIndex >( num2dElements ); + for( string const & fieldName : passthroughFieldNames ) + { + writeFacePassthroughField( *faceMesh, numElems, fieldName, faceBlock ); + } } } // end of namespace diff --git a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.hpp b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.hpp index f812071440d..b4efcb8d647 100644 --- a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.hpp @@ -19,6 +19,7 @@ #include "CellBlockManager.hpp" #include "common/DataTypes.hpp" +#include "common/Span.hpp" #include #include @@ -32,11 +33,13 @@ namespace geos::vtk * @param faceMesh[in] The vtk mesh for the face block. * @param mesh[in] The vtk volumic mesh. * @param cellBlockManager[inout] The cell block manager that will receive the face block information. + * @param passthroughFieldNames[in] Names of VTK CellData arrays to import and pass through to output as-is. */ void importFractureNetwork( string const & faceBlockName, vtkSmartPointer< vtkDataSet > faceMesh, vtkSmartPointer< vtkDataSet > mesh, - CellBlockManager & cellBlockManager ); + CellBlockManager & cellBlockManager, + Span< string const > passthroughFieldNames = {} ); } #endif // include guard diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp index 4a881414782..c9d0d708173 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp @@ -20,6 +20,7 @@ #include "VTKMeshGenerator.hpp" #include "common/DataTypes.hpp" +#include "common/Span.hpp" #include "mesh/ExternalDataSourceManager.hpp" #include "mesh/LogLevelsInfo.hpp" #include "mesh/generators/VTKFaceBlockUtilities.hpp" @@ -98,6 +99,11 @@ VTKMeshGenerator::VTKMeshGenerator( string const & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Name of the VTK data source" ); + registerWrapper( viewKeyStruct::passthroughFieldsString(), &m_passthroughFieldNames ). + setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Names of VTK CellData arrays to import and pass through to output as-is, without mapping to a GEOS field." ); + addLogLevel< logInfo::VTKSteps >(); } @@ -234,7 +240,8 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager writeNodes( getLogLevel(), *m_vtkMesh, m_nodesetNames, cellBlockManager, m_translate, m_scale ); GEOS_LOG_LEVEL_RANK_0( logInfo::VTKSteps, GEOS_FMT( "{} '{}': writing cells...", catalogName(), getName() ) ); - writeCells( getLogLevel(), *m_vtkMesh, m_cellMap, m_structuredIndexAttributeName, cellBlockManager ); + writeCells( getLogLevel(), *m_vtkMesh, m_cellMap, m_structuredIndexAttributeName, cellBlockManager, + Span< string const >( m_passthroughFieldNames.begin(), m_passthroughFieldNames.end() ) ); GEOS_LOG_LEVEL_RANK_0( logInfo::VTKSteps, GEOS_FMT( "{} '{}': writing surfaces...", catalogName(), getName() ) ); writeSurfaces( getLogLevel(), *m_vtkMesh, m_cellMap, cellBlockManager ); @@ -248,7 +255,8 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager for( auto const & [name, mesh]: m_faceBlockMeshes ) { - vtk::importFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager ); + vtk::importFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager, + Span< string const >( m_passthroughFieldNames.begin(), m_passthroughFieldNames.end() ) ); } GEOS_LOG_LEVEL_RANK_0( logInfo::VTKSteps, GEOS_FMT( "{} '{}': done!", catalogName(), getName() ) ); diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp index e31ada4bb1e..5fc838121ed 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp @@ -117,6 +117,7 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase constexpr static char const * partitionFractureWeightString() { return "partitionFractureWeight"; } constexpr static char const * useGlobalIdsString() { return "useGlobalIds"; } constexpr static char const * dataSourceString() { return "dataSourceName"; } + constexpr static char const * passthroughFieldsString() { return "passThroughFields"; } }; struct groupKeyStruct @@ -180,6 +181,9 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase /// Repository of VTK objects VTKHierarchicalDataSource * m_dataSource; + /// Names of VTK CellData arrays to import and pass through to output as-is + string_array m_passthroughFieldNames; + }; } // namespace geos diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 2bf0d4f02e5..baa80e40166 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -3561,12 +3561,89 @@ void writeStructuredIndex( vtkDataSet & mesh, } ); } +static void writePassthroughField( vtkDataSet & mesh, + Span< vtkIdType const > const cellIds, + string const & fieldName, + CellBlock & cellBlock ) +{ + vtkDataArray * srcArray = mesh.GetCellData()->GetArray( fieldName.c_str() ); + if( srcArray == nullptr ) + { + GEOS_WARNING( GEOS_FMT( "Passthrough field '{}' not found in VTK CellData; skipping.", fieldName ) ); + return; + } + + localIndex const numElems = LvArray::integerConversion< localIndex >( cellIds.size() ); + integer const numComp = srcArray->GetNumberOfComponents(); + + // Map VTK float types to real64; everything else (integers, id types) to integer. + bool const isRealArray = ( srcArray->GetDataType() == VTK_FLOAT || srcArray->GetDataType() == VTK_DOUBLE ); + + if( isRealArray ) + { + if( numComp == 1 ) + { + array1d< real64 > & dst = cellBlock.addProperty< array1d< real64 > >( fieldName ); + dst.resize( numElems ); + cellBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + forAll< parallelHostPolicy >( numElems, [&dst, srcArray, cellIds]( localIndex const i ) + { + dst[i] = static_cast< real64 >( srcArray->GetTuple1( cellIds[i] ) ); + } ); + } + else + { + array2d< real64 > & dst = cellBlock.addProperty< array2d< real64 > >( fieldName ); + dst.resize( numElems, numComp ); + cellBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + forAll< parallelHostPolicy >( numElems, [&dst, srcArray, cellIds, numComp]( localIndex const i ) + { + for( integer c = 0; c < numComp; ++c ) + dst( i, c ) = static_cast< real64 >( srcArray->GetComponent( cellIds[i], c ) ); + } ); + } + } + else + { + if( numComp == 1 ) + { + array1d< integer > & dst = cellBlock.addProperty< array1d< integer > >( fieldName ); + dst.resize( numElems ); + cellBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + forAll< parallelHostPolicy >( numElems, [&dst, srcArray, cellIds]( localIndex const i ) + { + dst[i] = static_cast< integer >( srcArray->GetTuple1( cellIds[i] ) ); + } ); + } + else + { + array2d< integer > & dst = cellBlock.addProperty< array2d< integer > >( fieldName ); + dst.resize( numElems, numComp ); + cellBlock.getWrapperBase( fieldName ).setPlotLevel( dataRepository::PlotLevel::LEVEL_1 ); + forAll< parallelHostPolicy >( numElems, [&dst, srcArray, cellIds, numComp]( localIndex const i ) + { + for( integer c = 0; c < numComp; ++c ) + dst( i, c ) = static_cast< integer >( srcArray->GetComponent( cellIds[i], c ) ); + } ); + } + } +} + void writeCells( integer const logLevel, vtkDataSet & mesh, vtk::CellMapType const & cellMap, string const & structuredIndexAttributeName, - CellBlockManager & cellBlockManager ) + CellBlockManager & cellBlockManager, + Span< string const > passthroughFieldNames ) { + if( !passthroughFieldNames.empty() ) + { + string fieldList; + for( string const & n : passthroughFieldNames ) + fieldList += "'" + n + "' "; + GEOS_LOG_RANK_0( GEOS_FMT( "VTKMesh: propagating passthrough fields: {}", fieldList ) ); + } + // Creates a new cell block for each region and for each type of cell. for( auto const & typeRegions : cellMap ) { @@ -3595,6 +3672,11 @@ void writeCells( integer const logLevel, { writeStructuredIndex( mesh, structuredIndexAttributeName, cellIds, cellBlock ); } + + for( string const & fieldName : passthroughFieldNames ) + { + writePassthroughField( mesh, cellIds, fieldName, cellBlock ); + } } } } diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp index a0181fbfd8b..af919f43e93 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp @@ -288,12 +288,14 @@ void writeNodes( integer const logLevel, * @param[in] cellMap Map from the surfaces index to the list of cells in this surface in this rank. * @param[in] structuredIndexAttributeName name of the VTK cell array to use as "structured index" attribute (if non-empty) * @param[in] cellBlockManager The instance that stores the cell blocks. + * @param[in] passthroughFieldNames Names of VTK CellData arrays to import and pass through to output as-is. */ void writeCells( integer const logLevel, vtkDataSet & mesh, vtk::CellMapType const & cellMap, string const & structuredIndexAttributeName, - CellBlockManager & cellBlockManager ); + CellBlockManager & cellBlockManager, + Span< string const > passthroughFieldNames = {} ); /** * @brief Build the "surface" node sets from the surface information.