diff --git a/src/coreComponents/fileIO/CMakeLists.txt b/src/coreComponents/fileIO/CMakeLists.txt index eedca99b65c..8e448471b8f 100644 --- a/src/coreComponents/fileIO/CMakeLists.txt +++ b/src/coreComponents/fileIO/CMakeLists.txt @@ -25,6 +25,7 @@ Contains: set( fileIO_headers LogLevelsInfo.hpp Outputs/BlueprintOutput.hpp + Outputs/CommonMeshOutputUtilities.hpp Outputs/MemoryStatsOutput.hpp Outputs/OutputBase.hpp Outputs/OutputManager.hpp @@ -32,6 +33,7 @@ set( fileIO_headers Outputs/PythonOutput.hpp Outputs/RestartOutput.hpp Outputs/TimeHistoryOutput.hpp + Outputs/XDMFOutput.hpp timeHistory/HDFFile.hpp timeHistory/HistoryCollectionBase.hpp timeHistory/BufferedHistoryIO.hpp @@ -44,6 +46,7 @@ set( fileIO_headers # set( fileIO_sources Outputs/BlueprintOutput.cpp + Outputs/CommonMeshOutputUtilities.cpp Outputs/MemoryStatsOutput.cpp Outputs/OutputBase.cpp Outputs/OutputManager.cpp @@ -51,6 +54,7 @@ set( fileIO_sources Outputs/PythonOutput.cpp Outputs/RestartOutput.cpp Outputs/TimeHistoryOutput.cpp + Outputs/XDMFOutput.cpp timeHistory/HDFFile.cpp timeHistory/HistoryCollectionBase.cpp timeHistory/PackCollection.cpp @@ -128,4 +132,3 @@ if( GEOS_ENABLE_TESTS ) add_subdirectory( timeHistory/unitTests ) add_subdirectory( Outputs/unitTests ) endif() - diff --git a/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp b/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp index a3ca597d4a8..dd6124a5f79 100644 --- a/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp +++ b/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp @@ -17,101 +17,20 @@ * @file BlueprintOutput.cpp */ -/// Source includes #include "BlueprintOutput.hpp" +#include "CommonMeshOutputUtilities.hpp" #include "common/TimingMacros.hpp" +#include "common/MpiWrapper.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/MeshLevel.hpp" -#include "dataRepository/ConduitRestart.hpp" -// TPL includes #include #include -#include namespace geos { -namespace internal -{ - -/** - * @brief @return The Blueprint shape from the GEOSX element type string. - * @param elementType the elementType to look up. - */ -string toBlueprintShape( ElementType const elementType ) -{ - switch( elementType ) - { - case ElementType::Tetrahedron: return "tet"; - case ElementType::Hexahedron: return "hex"; - default: - { - GEOS_ERROR( "No Blueprint type for element type: " << elementType ); - return {}; - } - } -} - -static stdVector< int > getBlueprintNodeOrdering( ElementType const elementType ) -{ - // Same as VTK, but kept separate for flexibility - switch( elementType ) - { - case ElementType::Vertex: return { 0 }; - case ElementType::Line: return { 0, 1 }; - case ElementType::Triangle: return { 0, 1, 2 }; - case ElementType::Quadrilateral: return { 0, 1, 2, 3 }; // TODO check - case ElementType::Polygon: return { }; // TODO - case ElementType::Tetrahedron: return { 1, 0, 2, 3 }; - case ElementType::Pyramid: return { 0, 3, 2, 1, 4, 0, 0, 0 }; - case ElementType::Wedge: return { 0, 4, 2, 1, 5, 3, 0, 0 }; - case ElementType::Hexahedron: return { 0, 1, 3, 2, 4, 5, 7, 6 }; - case ElementType::Prism5: return { }; // TODO - case ElementType::Prism6: return { }; // TODO - case ElementType::Prism7: return { }; // TODO - case ElementType::Prism8: return { }; // TODO - case ElementType::Prism9: return { }; // TODO - case ElementType::Prism10: return { }; // TODO - case ElementType::Prism11: return { }; // TODO - case ElementType::Polyhedron: return { }; // TODO - } - return {}; -} - -/** - * @brief Outputs the element to node map of @p subRegion to @p connectivity in VTK order. - * @param subRegion The sub-region to output. - * @param connectivity The Conduit Node to output to. - */ -void reorderElementToNodeMap( CellElementSubRegion const & subRegion, conduit::Node & connectivity ) -{ - GEOS_MARK_FUNCTION; - - arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemToNodeMap = subRegion.nodeList(); - localIndex const numElems = elemToNodeMap.size( 0 ); - localIndex const numNodesPerElem = elemToNodeMap.size( 1 ); - stdVector< int > const vtkOrdering = getBlueprintNodeOrdering( subRegion.getElementType() ); - GEOS_ERROR_IF_NE( localIndex( vtkOrdering.size() ), numNodesPerElem ); - - constexpr int conduitTypeID = dataRepository::conduitTypeInfo< localIndex >::id; - conduit::DataType const dtype( conduitTypeID, elemToNodeMap.size() ); - connectivity.set( dtype ); - - localIndex * const reorderedConnectivity = connectivity.value(); - forAll< serialPolicy >( numElems, [reorderedConnectivity, numNodesPerElem, elemToNodeMap, &vtkOrdering] ( localIndex const i ) - { - for( localIndex j = 0; j < numNodesPerElem; ++j ) - { - reorderedConnectivity[ i * numNodesPerElem + j ] = elemToNodeMap( i, vtkOrdering[ j ] ); - } - } ); -} - -} /// namespace internal; - -/////////////////////////////////////////////////////////////////////////////////////////////////// BlueprintOutput::BlueprintOutput( string const & name, dataRepository::Group * const parent ): OutputBase( name, parent ) @@ -126,9 +45,28 @@ BlueprintOutput::BlueprintOutput( string const & name, setApplyDefaultValue( false ). setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "If true writes out data associated with every quadrature point." ); + + registerWrapper( "writeParallelFiles", &m_writeParallelFiles ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If non-zero, write one heavy-data file per rank." ); + + registerWrapper( "writeGhostObjects", &m_writeGhostObjects ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If zero, skip ghost-owned cell objects." ); + + registerWrapper( "writeFieldData", &m_writeFieldData ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If zero, skip all field data." ); + + registerWrapper( "writeGhostFieldData", &m_writeGhostFieldData ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If zero, skip fields that contain ghost metadata." ); } -/////////////////////////////////////////////////////////////////////////////////////////////////// bool BlueprintOutput::execute( real64 const time_n, real64 const GEOS_UNUSED_PARAM( dt ), integer const cycleNumber, @@ -144,173 +82,45 @@ bool BlueprintOutput::execute( real64 const time_n, MeshLevel const & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); conduit::Node meshRoot; - conduit::Node & mesh = meshRoot[ "mesh" ]; - conduit::Node & coordset = mesh[ "coordsets/nodes" ]; - conduit::Node & topologies = mesh[ "topologies" ]; - - mesh[ "state/time" ] = time_n; - mesh[ "state/cycle" ] = cycleNumber; - - addNodalData( meshLevel.getNodeManager(), coordset, topologies, mesh[ "fields" ] ); - - dataRepository::Group averagedElementData( "averagedElementData", this ); - addElementData( meshLevel.getElemManager(), coordset, topologies, mesh[ "fields" ], averagedElementData ); - - /// The Blueprint will complain if the fields node is present but empty. - if( mesh[ "fields" ].number_of_children() == 0 ) - { - mesh.remove( "fields" ); - } + outputUtilities::buildCommonMeshTree( meshLevel, + time_n, + cycleNumber, + outputUtilities::CommonMeshOptions{ + m_plotLevel, + m_outputFullQuadratureData, + m_writeGhostObjects, + m_writeFieldData, + m_writeGhostFieldData + }, + *this, + meshRoot ); + + conduit::Node const & mesh = meshRoot.fetch_existing( "mesh" ); - /// Verify that the mesh conforms to the Blueprint. conduit::Node info; GEOS_ASSERT_MSG( conduit::blueprint::verify( "mesh", meshRoot, info ), info.to_json() ); - /// Generate the Blueprint index. conduit::Node fileRoot; conduit::Node & index = fileRoot[ "blueprint_index/mesh" ]; - conduit::blueprint::mesh::generate_index( mesh, "mesh", MpiWrapper::commSize(), index ); + conduit::index_t const numberOfDomains = m_writeParallelFiles != 0 ? MpiWrapper::commSize() : 1; + conduit::blueprint::mesh::generate_index( mesh, "mesh", numberOfDomains, index ); - /// Verify that the index conforms to the Blueprint. info.reset(); GEOS_ASSERT_MSG( conduit::blueprint::mesh::index::verify( index, info ), info.to_json() ); - /// Write out the root index file, then write out the mesh. - string const completePath = GEOS_FMT( "{}/blueprintFiles/cycle_{:07}", getOutputDirectory(), cycleNumber ); - string const filePathForRank = dataRepository::writeRootFile( fileRoot, completePath ); - conduit::relay::io::save( meshRoot, filePathForRank, "hdf5" ); + outputUtilities::writeCommonMeshHDF5( meshRoot, + getOutputDirectory(), + cycleNumber, + outputUtilities::CommonHDFWriteOptions{ + "blueprintFiles", + m_writeParallelFiles + }, + &fileRoot ); } return false; } -/////////////////////////////////////////////////////////////////////////////////////////////////// -void BlueprintOutput::addNodalData( NodeManager const & nodeManager, - conduit::Node & coordset, - conduit::Node & topologies, - conduit::Node & fields ) -{ - GEOS_MARK_FUNCTION; - - /// Populate the coordset group - coordset[ "type" ] = "explicit"; - dataRepository::wrapperHelpers::populateMCArray( nodeManager.referencePosition(), - coordset[ "values" ], - { "x", "y", "z" } ); - - /// Create the points topology - string const coordsetName = coordset.name(); - conduit::Node & nodeTopology = topologies[ coordsetName ]; - nodeTopology[ "coordset" ] = coordsetName; - - /// TODO: Once VisIT supports the implicit "points" topology we can just do the following. - /// See https://github.com/visit-dav/visit/issues/4593 - // nodeTopology[ "type" ] = "points"; - - nodeTopology[ "type" ] = "unstructured"; - nodeTopology[ "elements/shape" ] = "point"; - conduit::Node & connectivity = nodeTopology[ "elements/connectivity" ]; - - localIndex const numNodes = nodeManager.size(); - constexpr int conduitTypeID = dataRepository::conduitTypeInfo< localIndex >::id; - conduit::DataType const dtype( conduitTypeID, numNodes ); - connectivity.set( dtype ); - - localIndex * const nodeIDs = connectivity.value(); - forAll< serialPolicy >( numNodes, [nodeIDs] ( localIndex const i ) - { - nodeIDs[ i ] = i; - } ); - - /// Write out the fields. - writeOutWrappersAsFields( nodeManager, fields, coordsetName ); -} - -/////////////////////////////////////////////////////////////////////////////////////////////////// -void BlueprintOutput::addElementData( ElementRegionManager const & elemRegionManager, - conduit::Node & coordset, - conduit::Node & topologies, - conduit::Node & fields, - dataRepository::Group & averagedElementData ) -{ - GEOS_MARK_FUNCTION; - - elemRegionManager.forElementSubRegionsComplete< CellElementSubRegion >( - [&] ( localIndex, localIndex, ElementRegionBase const & region, CellElementSubRegion const & subRegion ) - { - string const topologyName = region.getName() + "-" + subRegion.getName(); - - /// Create the topology representing the sub-region. - conduit::Node & topology = topologies[ topologyName ]; - topology[ "coordset" ] = coordset.name(); - topology[ "type" ] = "unstructured"; - topology[ "elements/shape" ] = internal::toBlueprintShape( subRegion.getElementType() ); - internal::reorderElementToNodeMap( subRegion, topology[ "elements/connectivity" ] ); - - /// Write out the fields. - writeOutWrappersAsFields( subRegion, fields, topologyName ); - - /// Write out the quadrature averaged constitutive data and the full data if requested. - Group & averagedSubRegionData = averagedElementData.registerGroup( topologyName ); - subRegion.getConstitutiveModels().forSubGroups( [&]( dataRepository::Group const & constitutiveModel ) - { - writeOutConstitutiveData( constitutiveModel, fields, topologyName, averagedSubRegionData ); - - if( m_outputFullQuadratureData ) - { - writeOutWrappersAsFields( constitutiveModel, fields, topologyName, constitutiveModel.getName() ); - } - } ); - } ); -} - -/////////////////////////////////////////////////////////////////////////////////////////////////// -void BlueprintOutput::writeOutWrappersAsFields( Group const & group, - conduit::Node & fields, - string const & topology, - string const & prefix ) -{ - GEOS_MARK_FUNCTION; - - group.forWrappers( [&] ( dataRepository::WrapperBase const & wrapper ) - { - if( wrapper.getPlotLevel() <= m_plotLevel && wrapper.sizedFromParent() ) - { - string const name = prefix.empty() ? wrapper.getName() : prefix + "-" + wrapper.getName(); - - // conduit::Node & field = fields[ name ]; - // field[ "association" ] = "element"; - // field[ "volume_dependent" ] = "false"; - // field[ "topology" ] = topology; - // wrapper.populateMCArray( field[ "values" ] ); - - /// TODO: Replace with code above once https://github.com/visit-dav/visit/issues/4637 is fixed and released. - wrapper.addBlueprintField( fields, name, topology ); - } - } ); -} - -/////////////////////////////////////////////////////////////////////////////////////////////////// -void BlueprintOutput::writeOutConstitutiveData( dataRepository::Group const & constitutiveModel, - conduit::Node & fields, - string const & topology, - dataRepository::Group & averagedSubRegionData ) -{ - GEOS_MARK_FUNCTION; - - Group & averagedConstitutiveData = averagedSubRegionData.registerGroup( constitutiveModel.getName() ); - - constitutiveModel.forWrappers( [&] ( dataRepository::WrapperBase const & wrapper ) - { - if( wrapper.getPlotLevel() <= m_plotLevel && wrapper.sizedFromParent() ) - { - string const fieldName = constitutiveModel.getName() + "-quadrature-averaged-" + wrapper.getName(); - averagedConstitutiveData.registerWrapper( wrapper.averageOverSecondDim( fieldName, averagedConstitutiveData ) ). - addBlueprintField( fields, fieldName, topology ); - } - } ); -} - REGISTER_CATALOG_ENTRY( OutputBase, BlueprintOutput, string const &, dataRepository::Group * const ) } /* namespace geos */ diff --git a/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp b/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp index 6c6adb34bbe..4aeef3720d4 100644 --- a/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp +++ b/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp @@ -25,11 +25,6 @@ namespace geos { -// Forward declarations -class MeshLevel; -class NodeManager; -class ElementRegionManager; - /** * @class BlueprintOutput * @brief A class for creating Conduit blueprint-based outputs. @@ -83,53 +78,23 @@ class BlueprintOutput : public OutputBase private: - /** - * @brief Create the Blueprint coordinate set, the nodal topology and register nodal fields. - * @param nodeManager The NodeManager to write out. - * @param coordset The Blueprint coordinate set to populate. - * @param topologies The Node containing all the topologies. - * @param fields The Node containing all the fields. - */ - void addNodalData( NodeManager const & nodeManager, - conduit::Node & coordset, - conduit::Node & topologies, - conduit::Node & fields ); - - /** - * @brief Create the blueprint topologies and register element fields. - * @param elemRegionManager The ElementRegionManager to write out. - * @param coordset The Blueprint coordinate set to associate the topologies with. - * @param topologies The Node containing all the topologies. - * @param fields The Node containing all the fields. - */ - void addElementData( ElementRegionManager const & elemRegionManager, - conduit::Node & coordset, - conduit::Node & topologies, - conduit::Node & fields, - dataRepository::Group & averagedElementData ); - - /** - * @brief Write out all the children with the appropriate plot level of @p group. - * @param group The group that contains the children to write out. - * @param fields The Blueprint "fields" Node. - * @param topology The name of the Blueprint "topology" that @p group is associated with. - * @param prefix The string to prepend to the name of the field. If not specified no prefix is used. - */ - void writeOutWrappersAsFields( Group const & group, - conduit::Node & fields, - string const & topology, - string const & prefix="" ); - - void writeOutConstitutiveData( dataRepository::Group const & constitutiveModel, - conduit::Node & fields, - string const & topology, - dataRepository::Group & averagedElementData ); - // Used to determine which fields to write out. dataRepository::PlotLevel m_plotLevel = dataRepository::PlotLevel::LEVEL_1; // If true will write out the full quadrature data, otherwise it is averaged over. int m_outputFullQuadratureData = 0; + + // If non-zero, write one heavy-data file per rank. + int m_writeParallelFiles = 1; + + // If zero, skip ghost-owned cell objects. + int m_writeGhostObjects = 1; + + // If zero, skip all field data. + int m_writeFieldData = 1; + + // If zero, skip fields that contain ghost metadata (e.g., ghostRank). + int m_writeGhostFieldData = 1; }; diff --git a/src/coreComponents/fileIO/Outputs/CommonMeshOutputUtilities.cpp b/src/coreComponents/fileIO/Outputs/CommonMeshOutputUtilities.cpp new file mode 100644 index 00000000000..8ac9784e5c7 --- /dev/null +++ b/src/coreComponents/fileIO/Outputs/CommonMeshOutputUtilities.cpp @@ -0,0 +1,461 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CommonMeshOutputUtilities.cpp + */ + +#include "CommonMeshOutputUtilities.hpp" + +#include "common/MpiWrapper.hpp" +#include "common/Path.hpp" +#include "common/TimingMacros.hpp" +#include "common/logger/Logger.hpp" +#include "dataRepository/ConduitRestart.hpp" +#include "dataRepository/wrapperHelpers.hpp" +#include "mesh/ElementRegionManager.hpp" +#include "mesh/ElementSubRegionBase.hpp" +#include "mesh/MeshLevel.hpp" +#include "mesh/NodeManager.hpp" +#include "mesh/ObjectManagerBase.hpp" + +#include +#include + +namespace geos +{ +namespace outputUtilities +{ +namespace internal +{ + +string toBlueprintShape( ElementType const elementType ) +{ + switch( elementType ) + { + case ElementType::Tetrahedron: return "tet"; + case ElementType::Hexahedron: return "hex"; + default: + { + GEOS_ERROR( "No Blueprint type for element type: " << elementType ); + return {}; + } + } +} + +static stdVector< int > getBlueprintNodeOrdering( ElementType const elementType ) +{ + switch( elementType ) + { + case ElementType::Vertex: return { 0 }; + case ElementType::Line: return { 0, 1 }; + case ElementType::Triangle: return { 0, 1, 2 }; + case ElementType::Quadrilateral: return { 0, 1, 2, 3 }; + case ElementType::Polygon: return { }; + case ElementType::Tetrahedron: return { 1, 0, 2, 3 }; + case ElementType::Pyramid: return { 0, 3, 2, 1, 4, 0, 0, 0 }; + case ElementType::Wedge: return { 0, 4, 2, 1, 5, 3, 0, 0 }; + case ElementType::Hexahedron: return { 0, 1, 3, 2, 4, 5, 7, 6 }; + case ElementType::Prism5: return { }; + case ElementType::Prism6: return { }; + case ElementType::Prism7: return { }; + case ElementType::Prism8: return { }; + case ElementType::Prism9: return { }; + case ElementType::Prism10: return { }; + case ElementType::Prism11: return { }; + case ElementType::Polyhedron: return { }; + } + return {}; +} + +bool isGhostField( string const & fieldName ) +{ + return fieldName == ObjectManagerBase::viewKeyStruct::ghostRankString(); +} + +void collectElementIndicesToWrite( CellElementSubRegion const & subRegion, + int const writeGhostObjects, + stdVector< localIndex > & elementIndices, + bool & contiguousPrefix ) +{ + localIndex const numElems = subRegion.size(); + elementIndices.clear(); + elementIndices.reserve( numElems ); + + if( writeGhostObjects != 0 ) + { + for( localIndex ei = 0; ei < numElems; ++ei ) + { + elementIndices.emplace_back( ei ); + } + contiguousPrefix = true; + return; + } + + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + for( localIndex ei = 0; ei < numElems; ++ei ) + { + if( ghostRank[ei] < 0 ) + { + elementIndices.emplace_back( ei ); + } + } + + contiguousPrefix = true; + for( localIndex i = 0; i < LvArray::integerConversion< localIndex >( elementIndices.size() ); ++i ) + { + if( elementIndices[i] != i ) + { + contiguousPrefix = false; + break; + } + } +} + +void reorderElementToNodeMap( CellElementSubRegion const & subRegion, + stdVector< localIndex > const & elementIndices, + conduit::Node & connectivity ) +{ + GEOS_MARK_FUNCTION; + + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemToNodeMap = subRegion.nodeList(); + localIndex const numElems = LvArray::integerConversion< localIndex >( elementIndices.size() ); + localIndex const numNodesPerElem = elemToNodeMap.size( 1 ); + + stdVector< int > const vtkOrdering = getBlueprintNodeOrdering( subRegion.getElementType() ); + GEOS_ERROR_IF_NE( localIndex( vtkOrdering.size() ), numNodesPerElem ); + + constexpr int conduitTypeID = dataRepository::conduitTypeInfo< localIndex >::id; + conduit::DataType const dtype( conduitTypeID, numElems * numNodesPerElem ); + connectivity.set( dtype ); + + localIndex * const reorderedConnectivity = connectivity.value(); + forAll< serialPolicy >( numElems, [reorderedConnectivity, numNodesPerElem, elemToNodeMap, &vtkOrdering, &elementIndices] ( localIndex const i ) + { + localIndex const elemIndex = elementIndices[i]; + for( localIndex j = 0; j < numNodesPerElem; ++j ) + { + reorderedConnectivity[ i * numNodesPerElem + j ] = elemToNodeMap( elemIndex, vtkOrdering[ j ] ); + } + } ); +} + +void trimFieldValuesNode( conduit::Node & valuesNode, localIndex const count ) +{ + conduit::DataType dtype = valuesNode.dtype(); + conduit::index_t const originalCount = dtype.number_of_elements(); + GEOS_ERROR_IF_LT( originalCount, conduit::index_t( count ) ); + if( originalCount == conduit::index_t( count ) ) + { + return; + } + + void * dataPtr = const_cast< void * >( valuesNode.data_ptr() ); + dtype.set_number_of_elements( count ); + valuesNode.set_external( dtype, dataPtr ); +} + +void trimNewFieldsToCount( conduit::Node & fields, + stdVector< string > const & newFieldNames, + localIndex const count, + bool const contiguousPrefix ) +{ + if( count < 0 ) + { + return; + } + + if( !contiguousPrefix ) + { + for( string const & fieldName : newFieldNames ) + { + fields.remove( fieldName ); + } + + GEOS_LOG_RANK_0( "Warning: dropped field output for a non-contiguous non-ghost element layout " + "(ghost filtering requested)." ); + return; + } + + for( string const & fieldName : newFieldNames ) + { + conduit::Node & field = fields.fetch_existing( fieldName ); + conduit::Node & values = field.fetch_existing( "values" ); + trimFieldValuesNode( values, count ); + } +} + +void writeOutWrappersAsFields( dataRepository::Group const & group, + conduit::Node & fields, + string const & topology, + CommonMeshOptions const & options, + localIndex const countForFieldTrim, + bool const contiguousPrefix, + string const & prefix = "" ) +{ + GEOS_MARK_FUNCTION; + + if( options.writeFieldData == 0 ) + { + return; + } + + group.forWrappers( [&] ( dataRepository::WrapperBase const & wrapper ) + { + if( wrapper.getPlotLevel() > options.plotLevel || !wrapper.sizedFromParent() ) + { + return; + } + + if( options.writeGhostFieldData == 0 && isGhostField( wrapper.getName() ) ) + { + return; + } + + string const name = prefix.empty() ? wrapper.getName() : prefix + "-" + wrapper.getName(); + + conduit::index_t const numChildrenBefore = fields.number_of_children(); + wrapper.addBlueprintField( fields, name, topology ); + + stdVector< string > newFieldNames; + newFieldNames.reserve( fields.number_of_children() - numChildrenBefore ); + for( conduit::index_t i = numChildrenBefore; i < fields.number_of_children(); ++i ) + { + newFieldNames.emplace_back( fields.child( i ).name() ); + } + trimNewFieldsToCount( fields, newFieldNames, countForFieldTrim, contiguousPrefix ); + } ); +} + +void writeOutConstitutiveData( dataRepository::Group const & constitutiveModel, + conduit::Node & fields, + string const & topology, + CommonMeshOptions const & options, + localIndex const countForFieldTrim, + bool const contiguousPrefix, + dataRepository::Group & averagedElementData ) +{ + GEOS_MARK_FUNCTION; + + if( options.writeFieldData == 0 ) + { + return; + } + + dataRepository::Group & averagedConstitutiveData = averagedElementData.registerGroup( constitutiveModel.getName() ); + + constitutiveModel.forWrappers( [&] ( dataRepository::WrapperBase const & wrapper ) + { + if( wrapper.getPlotLevel() > options.plotLevel || !wrapper.sizedFromParent() ) + { + return; + } + + if( options.writeGhostFieldData == 0 && isGhostField( wrapper.getName() ) ) + { + return; + } + + string const fieldName = constitutiveModel.getName() + "-quadrature-averaged-" + wrapper.getName(); + conduit::index_t const numChildrenBefore = fields.number_of_children(); + averagedConstitutiveData.registerWrapper( wrapper.averageOverSecondDim( fieldName, averagedConstitutiveData ) ). + addBlueprintField( fields, fieldName, topology ); + + stdVector< string > newFieldNames; + newFieldNames.reserve( fields.number_of_children() - numChildrenBefore ); + for( conduit::index_t i = numChildrenBefore; i < fields.number_of_children(); ++i ) + { + newFieldNames.emplace_back( fields.child( i ).name() ); + } + trimNewFieldsToCount( fields, newFieldNames, countForFieldTrim, contiguousPrefix ); + } ); +} + +void addNodalData( NodeManager const & nodeManager, + conduit::Node & coordset, + conduit::Node & topologies, + conduit::Node & fields, + CommonMeshOptions const & options ) +{ + GEOS_MARK_FUNCTION; + + coordset[ "type" ] = "explicit"; + dataRepository::wrapperHelpers::populateMCArray( nodeManager.referencePosition(), + coordset[ "values" ], + { "x", "y", "z" } ); + + string const coordsetName = coordset.name(); + conduit::Node & nodeTopology = topologies[ coordsetName ]; + nodeTopology[ "coordset" ] = coordsetName; + nodeTopology[ "type" ] = "unstructured"; + nodeTopology[ "elements/shape" ] = "point"; + conduit::Node & connectivity = nodeTopology[ "elements/connectivity" ]; + + localIndex const numNodes = nodeManager.size(); + constexpr int conduitTypeID = dataRepository::conduitTypeInfo< localIndex >::id; + conduit::DataType const dtype( conduitTypeID, numNodes ); + connectivity.set( dtype ); + + localIndex * const nodeIDs = connectivity.value(); + forAll< serialPolicy >( numNodes, [nodeIDs] ( localIndex const i ) + { + nodeIDs[i] = i; + } ); + + writeOutWrappersAsFields( nodeManager, fields, coordsetName, options, -1, true ); +} + +void addElementData( ElementRegionManager const & elemRegionManager, + conduit::Node & coordset, + conduit::Node & topologies, + conduit::Node & fields, + CommonMeshOptions const & options, + dataRepository::Group & averagedElementData ) +{ + GEOS_MARK_FUNCTION; + + elemRegionManager.forElementSubRegionsComplete< CellElementSubRegion >( + [&] ( localIndex, localIndex, ElementRegionBase const & region, CellElementSubRegion const & subRegion ) + { + stdVector< localIndex > elementIndices; + bool contiguousPrefix = true; + collectElementIndicesToWrite( subRegion, options.writeGhostObjects, elementIndices, contiguousPrefix ); + + if( elementIndices.empty() ) + { + return; + } + + string const topologyName = region.getName() + "-" + subRegion.getName(); + + conduit::Node & topology = topologies[ topologyName ]; + topology[ "coordset" ] = coordset.name(); + topology[ "type" ] = "unstructured"; + topology[ "elements/shape" ] = toBlueprintShape( subRegion.getElementType() ); + reorderElementToNodeMap( subRegion, elementIndices, topology[ "elements/connectivity" ] ); + + localIndex const numElemsToWrite = LvArray::integerConversion< localIndex >( elementIndices.size() ); + writeOutWrappersAsFields( subRegion, fields, topologyName, options, numElemsToWrite, contiguousPrefix ); + + dataRepository::Group & averagedSubRegionData = averagedElementData.registerGroup( topologyName ); + subRegion.getConstitutiveModels().forSubGroups( [&]( dataRepository::Group const & constitutiveModel ) + { + writeOutConstitutiveData( constitutiveModel, + fields, + topologyName, + options, + numElemsToWrite, + contiguousPrefix, + averagedSubRegionData ); + + if( options.outputFullQuadratureData ) + { + writeOutWrappersAsFields( constitutiveModel, + fields, + topologyName, + options, + numElemsToWrite, + contiguousPrefix, + constitutiveModel.getName() ); + } + } ); + } ); +} + +} // namespace internal + +void buildCommonMeshTree( MeshLevel const & meshLevel, + real64 const time_n, + integer const cycleNumber, + CommonMeshOptions const & options, + dataRepository::Group & parentGroup, + conduit::Node & meshRoot ) +{ + GEOS_MARK_FUNCTION; + + conduit::Node & mesh = meshRoot[ "mesh" ]; + conduit::Node & coordset = mesh[ "coordsets/nodes" ]; + conduit::Node & topologies = mesh[ "topologies" ]; + + mesh[ "state/time" ] = time_n; + mesh[ "state/cycle" ] = cycleNumber; + + internal::addNodalData( meshLevel.getNodeManager(), coordset, topologies, mesh[ "fields" ], options ); + + dataRepository::Group averagedElementData( "averagedElementData", &parentGroup ); + internal::addElementData( meshLevel.getElemManager(), + coordset, + topologies, + mesh[ "fields" ], + options, + averagedElementData ); + + if( mesh[ "fields" ].number_of_children() == 0 ) + { + mesh.remove( "fields" ); + } +} + +CommonHDFWriteResult writeCommonMeshHDF5( conduit::Node & meshRoot, + string const & outputDirectory, + integer const cycleNumber, + CommonHDFWriteOptions const & options, + conduit::Node * rootNode ) +{ + GEOS_MARK_FUNCTION; + + GEOS_ERROR_IF( options.outputSubdirectory.empty(), "Output subdirectory must not be empty." ); + + CommonHDFWriteResult result; + result.numberOfFiles = options.writeParallelFiles != 0 ? MpiWrapper::commSize() : 1; + result.fileRank = options.writeParallelFiles != 0 ? MpiWrapper::commRank() : 0; + result.wroteFileOnThisRank = options.writeParallelFiles != 0 || MpiWrapper::commRank() == 0; + result.cyclePath = GEOS_FMT( "{}/{}/cycle_{:07}", outputDirectory, options.outputSubdirectory, cycleNumber ); + result.filePathForRank = GEOS_FMT( "{}/rank_{:07}.hdf5", result.cyclePath, result.fileRank ); + result.hdfFileName = splitPath( result.filePathForRank ).second; + + conduit::Node fallbackRoot; + conduit::Node & root = rootNode != nullptr ? *rootNode : fallbackRoot; + + if( options.writeParallelFiles != 0 ) + { + result.filePathForRank = dataRepository::writeRootFile( root, result.cyclePath ); + result.hdfFileName = splitPath( result.filePathForRank ).second; + conduit::relay::io::save( meshRoot, result.filePathForRank, "hdf5" ); + return result; + } + + if( MpiWrapper::commRank() == 0 ) + { + makeDirsForPath( result.cyclePath ); + string const rootFileName = splitPath( result.cyclePath ).second; + + root[ "protocol/name" ] = "hdf5"; + root[ "protocol/version" ] = CONDUIT_VERSION; + root[ "number_of_files" ] = 1; + root[ "file_pattern" ] = rootFileName + "/rank_%07d.hdf5"; + root[ "number_of_trees" ] = 1; + root[ "tree_pattern" ] = "/"; + + conduit::relay::io::save( root, result.cyclePath + ".root", "hdf5" ); + conduit::relay::io::save( meshRoot, result.filePathForRank, "hdf5" ); + } + + MpiWrapper::barrier( MPI_COMM_GEOS ); + return result; +} + +} /* namespace outputUtilities */ +} /* namespace geos */ diff --git a/src/coreComponents/fileIO/Outputs/CommonMeshOutputUtilities.hpp b/src/coreComponents/fileIO/Outputs/CommonMeshOutputUtilities.hpp new file mode 100644 index 00000000000..95b41216525 --- /dev/null +++ b/src/coreComponents/fileIO/Outputs/CommonMeshOutputUtilities.hpp @@ -0,0 +1,113 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CommonMeshOutputUtilities.hpp + */ + +#ifndef GEOS_FILEIO_OUTPUTS_COMMONMESHOUTPUTUTILITIES_HPP_ +#define GEOS_FILEIO_OUTPUTS_COMMONMESHOUTPUTUTILITIES_HPP_ + +#include "common/DataTypes.hpp" +#include "dataRepository/RestartFlags.hpp" + +namespace conduit +{ +class Node; +} + +namespace geos +{ + +class MeshLevel; + +namespace dataRepository +{ +class Group; +} + +namespace outputUtilities +{ + +/** + * @brief Common options for building a Conduit mesh tree for post-processing metadata. + */ +struct CommonMeshOptions +{ + dataRepository::PlotLevel plotLevel = dataRepository::PlotLevel::LEVEL_1; + int outputFullQuadratureData = 0; + int writeGhostObjects = 1; + int writeFieldData = 1; + int writeGhostFieldData = 1; +}; + +/** + * @brief Common options for writing the shared HDF5 heavy data. + */ +struct CommonHDFWriteOptions +{ + string outputSubdirectory; + int writeParallelFiles = 1; +}; + +/** + * @brief Output details returned by the shared HDF5 writer. + */ +struct CommonHDFWriteResult +{ + string cyclePath; + string filePathForRank; + string hdfFileName; + integer numberOfFiles = 0; + integer fileRank = 0; + int wroteFileOnThisRank = 0; +}; + +/** + * @brief Build a Conduit mesh tree with common GEOS mesh/field data. + * @param meshLevel The mesh level to serialize. + * @param time_n Current simulation time. + * @param cycleNumber Current cycle number. + * @param options Shared mesh output options. + * @param parentGroup Group used for temporary averaged constitutive fields. + * @param meshRoot Root node to populate. Expected layout: /mesh/... + */ +void buildCommonMeshTree( MeshLevel const & meshLevel, + real64 const time_n, + integer const cycleNumber, + CommonMeshOptions const & options, + dataRepository::Group & parentGroup, + conduit::Node & meshRoot ); + +/** + * @brief Write shared HDF5 heavy data and an associated root file. + * @param meshRoot Root Conduit node to write. + * @param outputDirectory Base output directory. + * @param cycleNumber Current cycle number. + * @param options Shared HDF write options. + * @param rootNode Optional pre-populated root node to write (.root file payload). + * @return Information about files written and file naming. + */ +CommonHDFWriteResult writeCommonMeshHDF5( conduit::Node & meshRoot, + string const & outputDirectory, + integer const cycleNumber, + CommonHDFWriteOptions const & options, + conduit::Node * rootNode = nullptr ); + +} /* namespace outputUtilities */ + +} /* namespace geos */ + +#endif /* GEOS_FILEIO_OUTPUTS_COMMONMESHOUTPUTUTILITIES_HPP_ */ diff --git a/src/coreComponents/fileIO/Outputs/XDMFOutput.cpp b/src/coreComponents/fileIO/Outputs/XDMFOutput.cpp new file mode 100644 index 00000000000..280b6dbffc2 --- /dev/null +++ b/src/coreComponents/fileIO/Outputs/XDMFOutput.cpp @@ -0,0 +1,250 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file XDMFOutput.cpp + */ + +#include "XDMFOutput.hpp" + +#include "CommonMeshOutputUtilities.hpp" +#include "common/Path.hpp" +#include "common/TimingMacros.hpp" +#include "mesh/DomainPartition.hpp" +#include "mesh/MeshLevel.hpp" + +#include + +#include +#include +#include + +namespace geos +{ +namespace internal +{ + +struct TopologyInfo +{ + string xdmfType; + localIndex nodesPerElement; +}; + +TopologyInfo toXdmfTopology( string const & blueprintShape ) +{ + if( blueprintShape == "hex" ) + { + return { "Hexahedron", 8 }; + } + if( blueprintShape == "tet" ) + { + return { "Tetrahedron", 4 }; + } + + GEOS_ERROR( "Unsupported Blueprint shape for XDMF: " << blueprintShape ); + return {}; +} + +} // namespace internal + +XDMFOutput::XDMFOutput( string const & name, + dataRepository::Group * const parent ): + OutputBase( name, parent ) +{ + registerWrapper( "plotLevel", &m_plotLevel ). + setApplyDefaultValue( dataRepository::PlotLevel::LEVEL_1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setRTTypeName( rtTypes::CustomTypes::plotLevel ). + setDescription( "Determines which fields to write." ); + + registerWrapper( "outputFullQuadratureData", &m_outputFullQuadratureData ). + setApplyDefaultValue( false ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If true writes out data associated with every quadrature point." ); + + registerWrapper( "writeParallelFiles", &m_writeParallelFiles ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If non-zero, write one heavy-data file per rank." ); + + registerWrapper( "writeGhostObjects", &m_writeGhostObjects ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If zero, skip ghost-owned cell objects." ); + + registerWrapper( "writeFieldData", &m_writeFieldData ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If zero, skip all field data." ); + + registerWrapper( "writeGhostFieldData", &m_writeGhostFieldData ). + setApplyDefaultValue( 1 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setDescription( "If zero, skip fields that contain ghost metadata." ); +} + +bool XDMFOutput::execute( real64 const time_n, + real64 const GEOS_UNUSED_PARAM( dt ), + integer const cycleNumber, + integer const GEOS_UNUSED_PARAM( eventCounter ), + real64 const GEOS_UNUSED_PARAM( eventProgress ), + DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + + { + Timer timer( m_outputTimer ); + + MeshLevel const & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); + + conduit::Node meshRoot; + outputUtilities::buildCommonMeshTree( meshLevel, + time_n, + cycleNumber, + outputUtilities::CommonMeshOptions{ + m_plotLevel, + m_outputFullQuadratureData, + m_writeGhostObjects, + m_writeFieldData, + m_writeGhostFieldData + }, + *this, + meshRoot ); + + outputUtilities::CommonHDFWriteResult const writeResult = + outputUtilities::writeCommonMeshHDF5( meshRoot, + getOutputDirectory(), + cycleNumber, + outputUtilities::CommonHDFWriteOptions{ + "xdmfFiles", + m_writeParallelFiles + } ); + + if( writeResult.wroteFileOnThisRank == 0 ) + { + return false; + } + + string const xdmfPath = GEOS_FMT( "{}/rank_{:07}.xmf", writeResult.cyclePath, writeResult.fileRank ); + std::ofstream xdmfFile( xdmfPath ); + GEOS_ERROR_IF( !xdmfFile.is_open(), GEOS_FMT( "Failed to open XDMF file for writing: {}", xdmfPath ) ); + xdmfFile << buildXdmfDocument( meshRoot.fetch_existing( "mesh" ), writeResult.hdfFileName, time_n ); + } + + return false; +} + +string XDMFOutput::buildXdmfDocument( conduit::Node const & mesh, + string const & hdfFileName, + real64 const time ) +{ + conduit::Node const & coordset = mesh.fetch_existing( "coordsets/nodes" ); + conduit::Node const & topologies = mesh.fetch_existing( "topologies" ); + conduit::Node const * fields = mesh.has_path( "fields" ) ? &mesh.fetch_existing( "fields" ) : nullptr; + + string const coordsetName = "nodes"; + string const xPath = hdfFileName + ":/mesh/coordsets/nodes/values/x"; + string const yPath = hdfFileName + ":/mesh/coordsets/nodes/values/y"; + string const zPath = hdfFileName + ":/mesh/coordsets/nodes/values/z"; + localIndex const numberOfNodes = coordset.fetch_existing( "values/x" ).dtype().number_of_elements(); + + std::unordered_map< string, stdVector< conduit::Node const * > > elementFieldsByTopology; + stdVector< conduit::Node const * > nodalFields; + if( fields != nullptr ) + { + for( conduit::index_t i = 0; i < fields->number_of_children(); ++i ) + { + conduit::Node const & field = fields->child( i ); + string const topology = field.fetch_existing( "topology" ).as_string(); + if( topology == coordsetName ) + { + nodalFields.emplace_back( &field ); + } + else + { + elementFieldsByTopology[topology].emplace_back( &field ); + } + } + } + + std::ostringstream os; + os << "\n"; + os << "\n"; + os << "\n"; + os << " \n"; + + for( conduit::index_t i = 0; i < topologies.number_of_children(); ++i ) + { + conduit::Node const & topology = topologies.child( i ); + string const topologyName = topology.name(); + string const shape = topology.fetch_existing( "elements/shape" ).as_string(); + if( shape == "point" ) + { + continue; + } + + internal::TopologyInfo const topoInfo = internal::toXdmfTopology( shape ); + conduit::Node const & connectivity = topology.fetch_existing( "elements/connectivity" ); + localIndex const totalConnectivitySize = connectivity.dtype().number_of_elements(); + localIndex const numberOfElements = totalConnectivitySize / topoInfo.nodesPerElement; + string const connectivityPath = hdfFileName + ":/mesh/topologies/" + topologyName + "/elements/connectivity"; + + os << " \n"; + os << " \n"; + } + + os << " \n"; + os << "\n"; + return os.str(); +} + +REGISTER_CATALOG_ENTRY( OutputBase, XDMFOutput, string const &, dataRepository::Group * const ) + +} /* namespace geos */ diff --git a/src/coreComponents/fileIO/Outputs/XDMFOutput.hpp b/src/coreComponents/fileIO/Outputs/XDMFOutput.hpp new file mode 100644 index 00000000000..b75f90c8472 --- /dev/null +++ b/src/coreComponents/fileIO/Outputs/XDMFOutput.hpp @@ -0,0 +1,84 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file XDMFOutput.hpp + */ + +#ifndef GEOS_FILEIO_OUTPUTS_XDMFOUTPUT_HPP_ +#define GEOS_FILEIO_OUTPUTS_XDMFOUTPUT_HPP_ + +#include "fileIO/Outputs/OutputBase.hpp" + +namespace conduit +{ +class Node; +} + +namespace geos +{ + +/** + * @class XDMFOutput + * @brief A class for creating XDMF outputs with heavy data stored in HDF5. + */ +class XDMFOutput : public OutputBase +{ +public: + + XDMFOutput( string const & name, + Group * const parent ); + + virtual ~XDMFOutput() override + {} + + static string catalogName() { return "XDMF"; } + + virtual bool execute( real64 const time_n, + real64 const dt, + integer const cycleNumber, + integer const eventCounter, + real64 const eventProgress, + DomainPartition & domain ) override; + + virtual void cleanup( real64 const time_n, + integer const cycleNumber, + integer const eventCounter, + real64 const eventProgress, + DomainPartition & domain ) override + { execute( time_n, 0, cycleNumber, eventCounter, eventProgress, domain ); } + +private: + + static string buildXdmfDocument( conduit::Node const & mesh, + string const & hdfFileName, + real64 time ); + + dataRepository::PlotLevel m_plotLevel = dataRepository::PlotLevel::LEVEL_1; + + int m_outputFullQuadratureData = 0; + + int m_writeParallelFiles = 1; + + int m_writeGhostObjects = 1; + + int m_writeFieldData = 1; + + int m_writeGhostFieldData = 1; +}; + +} // namespace geos + +#endif // GEOS_FILEIO_OUTPUTS_XDMFOUTPUT_HPP_