From 53aa3b638f4f95da9a8eeab225f4674da09a064e Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 26 Jun 2025 20:45:12 -0500 Subject: [PATCH 01/30] Draft --- src/cmake/GeosxOptions.cmake | 2 + .../physicsSolvers/CMakeLists.txt | 7 +- .../physicsSolvers/gravity/CMakeLists.txt | 48 +++ .../physicsSolvers/gravity/GravityFE.cpp | 364 ++++++++++++++++++ .../physicsSolvers/gravity/GravityFE.hpp | 136 +++++++ .../gravity/GravityFEKernel.hpp | 165 ++++++++ .../gravity/GravitySolverBase.cpp | 215 +++++++++++ .../gravity/GravitySolverBase.hpp | 137 +++++++ 8 files changed, 1073 insertions(+), 1 deletion(-) create mode 100644 src/coreComponents/physicsSolvers/gravity/CMakeLists.txt create mode 100644 src/coreComponents/physicsSolvers/gravity/GravityFE.cpp create mode 100644 src/coreComponents/physicsSolvers/gravity/GravityFE.hpp create mode 100644 src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp create mode 100644 src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp create mode 100644 src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp diff --git a/src/cmake/GeosxOptions.cmake b/src/cmake/GeosxOptions.cmake index c9072b081a4..265a6cd23bd 100644 --- a/src/cmake/GeosxOptions.cmake +++ b/src/cmake/GeosxOptions.cmake @@ -126,6 +126,8 @@ option( GEOS_ENABLE_SIMPLEPDE "Enables simple PDE physics package" ON ) option( GEOS_ENABLE_SOLIDMECHANICS "Enables solid mechanics physics package" ON ) option( GEOS_ENABLE_SURFACEGENERATION "Enables surface generation physics package" ON ) option( GEOS_ENABLE_WAVEPROPAGATION "Enables wave propagation physics package" ON ) +option( GEOS_ENABLE_GRAVITY "Enables gravity physics package" ON ) + #set(CMAKE_POSITION_INDEPENDENT_CODE ON CACHE BOOL "" FORCE) #blt_append_custom_compiler_flag(FLAGS_VAR CMAKE_CXX_FLAGS DEFAULT -rdynamic) diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index 246fc46be5f..1fcbd15f371 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -113,6 +113,11 @@ if( GEOS_ENABLE_WAVEPROPAGATION ) list( APPEND dependencyList wavePropagationSolvers ) endif() +if ( GEOS_ENABLE_GRAVITY ) + add_subdirectory( gravity ) + list( APPEND dependencyList gravitySolvers ) +endif() + if( GEOS_ENABLE_MULTIPHYSICS ) add_subdirectory( multiphysics ) list( APPEND dependencyList multiPhysicsSolvers ) @@ -134,4 +139,4 @@ if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU") target_link_options(physicsSolvers PUBLIC -Wl,--no-as-needed) endif() -install( TARGETS physicsSolvers LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/lib ) \ No newline at end of file +install( TARGETS physicsSolvers LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/lib ) diff --git a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt new file mode 100644 index 00000000000..a05504c0211 --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt @@ -0,0 +1,48 @@ +# 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. +# +#-------------------------------------------------------------------------------------------------- + +#[[ +Package: gravitySolvers +Contains: +#]] + +# Specify solver headers +set( gravitySolvers_headers + GravitySolverBase.hpp + GravityFE.hpp + GravityFEKernel.hpp ) + +# Specify solver sources +set( gravitySolvers_sources + GravitySolverBase.cpp + GravityFE.cpp ) + +set( dependencyList ${parallelDeps} physicsSolversBase ) + +geos_decorate_link_dependencies( LIST decoratedDependencies DEPENDENCIES ${dependencyList} ) + +blt_add_library( NAME gravitySolvers + SOURCES ${gravitySolvers_sources} + HEADERS ${gravitySolvers_headers} + DEPENDS_ON ${decoratedDependencies} ${externalComponentDeps} + OBJECT ${GEOS_BUILD_OBJ_LIBS} + SHARED ${GEOS_BUILD_SHARED_LIBS} + ) + +target_include_directories( gravitySolvers PUBLIC ${CMAKE_SOURCE_DIR}/coreComponents ) + +install( TARGETS gravitySolvers LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/lib ) + +if( externalComponentDeps ) + target_include_directories( gravitySolvers PUBLIC ${CMAKE_SOURCE_DIR}/externalComponents ) +endif() diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp new file mode 100644 index 00000000000..045df4cfb5f --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -0,0 +1,364 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 GravityFE.cpp + */ +#include + +#include "GravityFE.hpp" +#include "GravityFEKernel.hpp" + +#include "discretizationMethods/NumericalMethodsManager.hpp" +#include "finiteElement/FiniteElementDiscretization.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/ElementType.hpp" +#include "mesh/DomainPartition.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" + + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; + +constexpr real64 GRAVITATIONAL_CONSTANT = 6.67430e-11; // in m3.kg-1.s-2 // Older value: 6.6738480e-11 +constexpr localIndex MAX_SUPPORT_POINTS = 8; + + + +GravityFE::GravityFE( const std::string & name, + Group * const parent ): + GravitySolverBase( name, parent ) +{} + +GravityFE::~GravityFE() +{} + + +void GravityFE::initializePreSubGroups() +{ + GravitySolverBase::initializePreSubGroups(); + + DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteElementDiscretizationManager const & + feDiscretizationManager = numericalMethodManager.getFiniteElementDiscretizationManager(); + + FiniteElementDiscretization const * const + feDiscretization = feDiscretizationManager.getGroupPointer< FiniteElementDiscretization >( m_discretizationName ); + GEOS_THROW_IF( feDiscretization == nullptr, + getName() << ": FE discretization not found: " << m_discretizationName, + InputError ); +} + + +void GravityFE::registerDataOnMesh( Group & meshBodies ) +{ + + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + string_array const & ) + { + NodeManager & nodeManager = mesh.getNodeManager(); + + + nodeManager.registerField< fields::VolumeIntegral >( this->getName() ); + + ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion ) + { + subRegion.registerField< fields::MediumDensity >( this->getName() ); + subRegion.registerField< fields::Adjoint >( this->getName() ); + subRegion.registerField< fields::VolumeIntegral2d >( this->getName() ); + + // Assume the maximum number of points per cell is MAX_SUPPORT_POINTS... + subRegion.getField< fields::VolumeIntegral2d >().resizeDimension< 1 >( MAX_SUPPORT_POINTS ); + } ); + } ); +} + + +void GravityFE::postInputInitialization() +{ + GravitySolverBase::postInputInitialization(); +} + + +void GravityFE::initializePostInitialConditionsPreSubGroups() +{ + GravitySolverBase::initializePostInitialConditionsPreSubGroups(); +} + + +real64 GravityFE::explicitStepModeling( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + + int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); + + array1d< real64 > localGzAtStations( m_stationCoordinates.size( 0 ) ); + localGzAtStations.setValues< parallelHostPolicy >( 0. ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), + [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + // Step 1: Precompute volumeIntegral once and for all. + NodeManager & nodeManager = mesh.getNodeManager(); + arrayView1d< integer const > const & nodeGhostRank = nodeManager.ghostRank(); + + // VolumeIntegral matrix to be computed in this function. + arrayView1d< real64 > const volumeIntegral = nodeManager.getField< fields::VolumeIntegral >(); + volumeIntegral.setValues< parallelHostPolicy >( 0. ); + + +#if 0 + // Print size and type of "CellElements" partitions. + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + GEOS_LOG_RANK_0( "GravityFE: Rank 0: Partition[ " << elementSubRegion.getElementType() << " ] = "<< elementSubRegion.size()); + } ); +#endif + + // Loop over all sub-regions in regions of type "CellElements". + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()); + if( this->getLogLevel()>3 ) + { + for( localIndex i=0; i const gzAtStations = m_gzAtStations.toView(); + for( localIndex istation=0; istation gzAtStations = m_gzAtStations.toView(); + gzAtStations = localGzAtStations; +#endif + + + for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) + { + auto const & coords = m_stationCoordinates[iStation]; + std::ostringstream logStream; + logStream << std::fixed << std::setprecision( 2 ); + logStream << "GravityFE: station[" << std::setw( 5 ) << iStation << "] " + << std::setw( 15 ) << coords[0] << " " + << std::setw( 15 ) << coords[1] << " " + << std::setw( 10 ) << coords[2] << " " + << std::scientific << std::setprecision( 6 ) + << std::setw( 14 ) << gzAtStations[iStation]; + GEOS_LOG_RANK_0( logStream.str() ); + } + + + + // Dump result to disk... + if( (rank==0) && ( this->m_outputGz == 1 )) + { + GravitySolverBase::saveGz( time_n, cycleNumber, this->m_outputGzBasename, gzAtStations ); + } + + return dt; +} + + +real64 GravityFE::explicitStepAdjoint( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( time_n, cycleNumber ); + + //int const rank = MpiWrapper::commRank( MPI_COMM_GEOSX ); + + array1d< real64 > localGzAtStations( m_stationCoordinates.size( 0 ) ); + localGzAtStations.setValues< parallelHostPolicy >( 0. ); + + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + + auto const & nodePosition = nodeManager.referencePosition(); + + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + + arrayView1d< real64 const > const residue = m_residue.toViewConst(); + + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< real64 > volumeIntegral2d = elementSubRegion.getReference< array2d< real64 > >( fields::VolumeIntegral2d::key()); + arrayView1d< real64 > adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()); + adjoint.zero(); + + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition().toViewConst(); + + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + + localIndex const numSupportPoints = fe.getNumSupportPoints(); + GEOS_ERROR_IF_GT_MSG( numSupportPoints, MAX_SUPPORT_POINTS, GEOS_FMT( "Maximum number of SupportPoints is {}", MAX_SUPPORT_POINTS ) ); + + + finiteElement::dispatchlowOrder3D( fe, [&] ( auto const finiteElement ) + { + using FE_TYPE = TYPEOFREF( finiteElement ); + + gravityFEKernel:: + VolumeIntegralKernel_uni2< FE_TYPE > kernel( finiteElement ); + kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > + ( elementSubRegion.size(), + X, + elemsToNodes, + volumeIntegral2d ); + } ); + + + for( localIndex iStation=0; iStation( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + for( localIndex iLoc=0; iLoc < numSupportPoints; ++iLoc ) + { + localIndex a = elemsToNodes[k][iLoc]; + + real64 dx = nodePosition[a][0] - coords[0]; + real64 dy = nodePosition[a][1] - coords[1]; + real64 dz = nodePosition[a][2] - coords[2]; + real64 r2 = dx*dx + dy*dy + dz*dz; + real64 r3 = sqrt( r2 ) * r2; + adjoint[k] += GRAVITATIONAL_CONSTANT * volumeIntegral2d[k][iLoc] * res * dz / r3; + + } + + } ); // Loop elem + } //Loop station + + adjoint.move( LvArray::MemorySpace::host, true ); + + if( this->getLogLevel() > 2 ) + { + for( localIndex i=0; i; + using ATOMIC_POLICY = parallelDeviceAtomic; + + + GravityFE() = delete; + + GravityFE( const std::string & name, + Group * const parent ); + + virtual ~GravityFE() override; + + + GravityFE( GravityFE const & ) = delete; + GravityFE( GravityFE && ) = default; + + GravityFE & operator=( GravityFE const & ) = delete; + GravityFE & operator=( GravityFE && ) = delete; + + + static string catalogName() { return "GravityFE"; } + string getCatalogName() const override { return catalogName(); } + + + virtual void initializePreSubGroups() override; + + virtual void registerDataOnMesh( Group & meshBodies ) override final; + + + /** + * @defgroup Solver Interface Functions + * + * These functions provide the primary interface that is required for derived classes + */ + /**@{*/ + virtual + real64 explicitStepModeling( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override; + virtual + real64 explicitStepAdjoint( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override; + /**@}*/ + + +protected: + + virtual void postInputInitialization() override final; + + virtual void initializePostInitialConditionsPreSubGroups() override final; + + +private: + + +}; + + + +namespace fields +{ +DECLARE_FIELD( MediumDensity, + "mediumDensity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Medium density of the cell" ); + +DECLARE_FIELD( VolumeIntegral, + "volumeIntegral", + array1d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "VolumeIntegral of the cell." ); + +DECLARE_FIELD( VolumeIntegral2d, + "volumeIntegral2d", + array2d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "VolumeIntegral for adjoint computation." ); + +DECLARE_FIELD( Adjoint, + "adjoint", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Adjoint field." ); +} + + + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp new file mode 100644 index 00000000000..3c41586f7de --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp @@ -0,0 +1,165 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 GravityFEKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ +#define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ + +//#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" + +namespace geos +{ + +/// Namespace to contain the gravity kernels. +namespace gravityFEKernel +{ + +template< typename FE_TYPE > +struct DensityVolumeIntegralKernel +{ + explicit DensityVolumeIntegralKernel( FE_TYPE const & finiteElement ) + : m_finiteElement( finiteElement ) + {} + + /** + * @brief Launches the precomputation of the DensityVolumeIntegral matrix. + * @tparam EXEC_POLICY the execution policy + * @tparam ATOMIC_POLICY the atomic policy + * @param[in] size Number of elements in the subregion + * @param[in] X Node coordinates + * @param[in] elemsToNodes Element-to-node connectivity + * @param[in] density Cell-wise density + * @param[out] volumeIntegral Output mass matrix diagonal + */ + template< typename EXEC_POLICY, typename ATOMIC_POLICY > + void launch( localIndex const size, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes, + arrayView1d< real64 const > const density, + arrayView1d< real64 > const volumeIntegral ) const + { + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + constexpr localIndex numNodesPerElem = FE_TYPE::numNodes; + constexpr localIndex numQuadraturePoints = FE_TYPE::numQuadraturePoints; + + auto const rho = density[k]; + + // Gather nodal coordinates for this element + real64 xLocal[numNodesPerElem][3]; + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + localIndex const nodeIdx = elemsToNodes( k, a ); + for( localIndex i = 0; i < 3; ++i ) + { + xLocal[a][i] = X( nodeIdx, i ); + } + } + + real64 N[numNodesPerElem]; + real64 gradN[numNodesPerElem][3]; + + for( localIndex q = 0; q < numQuadraturePoints; ++q ) + { + FE_TYPE::calcN( q, N ); + real64 const detJ = m_finiteElement.template getGradN< FE_TYPE >( k, q, xLocal, gradN ); + real64 const weightedJacobian = rho * detJ; + + for( localIndex a = 0; a < numNodesPerElem; ++a ) + { + localIndex const nodeIdx = elemsToNodes( k, a ); + real64 const localIncrement = weightedJacobian * N[a]; + RAJA::atomicAdd< ATOMIC_POLICY >( &volumeIntegral[nodeIdx], localIncrement ); + } + } + } ); + } + +private: + FE_TYPE const & m_finiteElement; +}; + + + +template< typename FE_TYPE > +struct VolumeIntegralKernel_uni2 +{ + explicit VolumeIntegralKernel_uni2( FE_TYPE const & finiteElement ) + : m_finiteElement( finiteElement ) + {} + + /** + * @brief Launches the precomputation of the volume integral matrix. + * @tparam EXEC_POLICY the execution policy + * @tparam ATOMIC_POLICY the atomic policy + * @param[in] size Number of elements in the subregion + * @param[in] X Node coordinates + * @param[in] elemsToNodes Element-to-node connectivity + * @param[out] volumeIntegral Output volume integral matrix + */ + template< typename EXEC_POLICY, typename ATOMIC_POLICY > + void launch( localIndex const size, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X, + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes, + arrayView2d< real64 > const volumeIntegral ) const + { + forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + constexpr localIndex numNodes = FE_TYPE::numNodes; + constexpr localIndex numQuadPoints = FE_TYPE::numQuadraturePoints; + + // Gather nodal coordinates for this element + real64 xLocal[numNodes][3]; + for( localIndex a = 0; a < numNodes; ++a ) + { + localIndex const nodeIdx = elemsToNodes( k, a ); + for( localIndex i = 0; i < 3; ++i ) + { + xLocal[a][i] = X( nodeIdx, i ); + } + } + + real64 N[numNodes]; + real64 gradN[numNodes][3]; + + for( localIndex q = 0; q < numQuadPoints; ++q ) + { + FE_TYPE::calcN( q, N ); + real64 const detJ = m_finiteElement.template getGradN< FE_TYPE >( k, q, xLocal, gradN ); + + for( localIndex a = 0; a < numNodes; ++a ) + { + real64 const localIncrement = detJ * N[a]; + RAJA::atomicAdd< ATOMIC_POLICY >( &volumeIntegral( k, a ), localIncrement ); + } + } + } ); + } + +private: + FE_TYPE const & m_finiteElement; +}; + + + +} // namespace gravityFEKernel + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp new file mode 100644 index 00000000000..a2867c5c9ab --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp @@ -0,0 +1,215 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 GravitySolverBase.cpp + */ + +#include "GravitySolverBase.hpp" + +#include "dataRepository/KeyNames.hpp" +#include "finiteElement/FiniteElementDiscretization.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" +#include + +namespace geos +{ + +using namespace dataRepository; + + + +const std::unordered_map< std::string, GravitySolverBase::GravityMode > GravitySolverBase::modeMap = { + { "modeling", GravitySolverBase::GravityMode::Modeling }, + { "adjoint", GravitySolverBase::GravityMode::Adjoint } +}; + + +GravitySolverBase::GravitySolverBase( const std::string & name, Group * const parent ) + : PhysicsSolverBase( name, parent ) +{ + registerWrapper( viewKeyStruct::modeString(), &m_modeString ) + .setInputFlag( InputFlags::OPTIONAL ) + .setApplyDefaultValue( "modeling" ) + .setDescription( "Mode: 'modeling' (default) or 'adjoint'" ); + + registerWrapper( viewKeyStruct::stationCoordinatesString(), &m_stationCoordinates ) + .setInputFlag( InputFlags::REQUIRED ) + .setSizedFromParent( 0 ) + .setDescription( "Coordinates (x,y,z) of the gravimeter stations" ); + + registerWrapper( viewKeyStruct::outputGzString(), &m_outputGz ) + .setInputFlag( InputFlags::OPTIONAL ) + .setApplyDefaultValue( 1 ) + .setDescription( "Flag to dump to disk field recorded at gravimeters" ); + + registerWrapper( viewKeyStruct::outputGzBasenameString(), &m_outputGzBasename ) + .setInputFlag( InputFlags::OPTIONAL ) + .setApplyDefaultValue( "gz" ) + .setDescription( "Basename to dump to disk field recorded at gravimeters" ); + + registerWrapper( viewKeyStruct::residueString(), &m_residue ) + .setInputFlag( InputFlags::OPTIONAL ) + .setSizedFromParent( 0 ) + .setDescription( "Residue at each station" ); + + registerWrapper( viewKeyStruct::gzAtStationsString(), &m_gzAtStations ) + .setInputFlag( InputFlags::FALSE ) + .setSizedFromParent( 0 ) + .setDescription( "Gz value at each station" ); +} + +GravitySolverBase::~GravitySolverBase() = default; + +void GravitySolverBase::reinit() +{ + initializePostInitialConditionsPreSubGroups(); +} + +void GravitySolverBase::initializePreSubGroups() +{ + PhysicsSolverBase::initializePreSubGroups(); +} + +void GravitySolverBase::initializePostInitialConditionsPreSubGroups() +{ + PhysicsSolverBase::initializePostInitialConditionsPreSubGroups(); + + + try + { + m_mode = parseMode( m_modeString ); + } + catch( const std::invalid_argument & e ) + { + GEOS_THROW( e.what(), InputError ); + } + + if( m_mode == GravityMode::Adjoint ) + { + GEOS_THROW_IF( m_stationCoordinates.size( 0 ) != m_residue.size( 0 ), + "GravitySolverBase: Residue size does not match the number of stations", + InputError ); + } + + + +#if 0 + if( m_modeString == "modeling" ) + { + m_mode = GravityMode_Modeling; + } + else if( m_modeString == "adjoint" ) + { + m_mode = GravityMode_Adjoint; + GEOS_THROW_IF( m_stationCoordinates.size( 0 ) != m_residue.size( 0 ), + "GravitySolverBase: Residue size does not match the number of stations", + InputError ); + } + else + { + GEOS_THROW( "Invalid mode string: " + m_modeString, InputError ); + } +#endif + constexpr auto yesno = [] (int flag) noexcept->const char * { return flag ? "yes" : "no"; }; + + LogPart gravitySolverLog( "Gravity Solver: ", MpiWrapper::commRank() == 0 ); + gravitySolverLog.begin(); + GEOS_LOG_RANK_0( "Name: " << getName()); + GEOS_LOG_RANK_0( "Mode: " << m_modeString ); + GEOS_LOG_RANK_0( "Number of stations: " << m_stationCoordinates.size( 0 )); + GEOS_LOG_RANK_0( "Output Gz to file: " << yesno( m_outputGz )); + GEOS_LOG_RANK_0( "Output Gz basename: " << m_outputGzBasename ); + GEOS_LOG_RANK_0( "Log level: " << getLogLevel()); + GEOS_LOG_RANK_0( " Output Gz to logs: " << yesno( getLogLevel() > 1 )); + GEOS_LOG_RANK_0( " Output Adjoint to logs: " << yesno( getLogLevel() > 2 )); + GEOS_LOG_RANK_0( " Output Properties to logs: " << yesno( getLogLevel() > 3 )); + gravitySolverLog.end(); +} + +void GravitySolverBase::postInputInitialization() +{ + PhysicsSolverBase::postInputInitialization(); + + GEOS_THROW_IF( m_stationCoordinates.size( 1 ) != 3, + "GravitySolverBase: Invalid number of physical coordinates for the gravimeter stations", + InputError ); + + m_gzAtStations.resize( m_stationCoordinates.size( 0 )); +} + +real64 GravitySolverBase::solverStep( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) +{ + GEOS_LOG_RANK_0( "GravitySolverBase: time_n: " << time_n << " dt: " << dt ); + return explicitStep( time_n, dt, cycleNumber, domain ); +} + +real64 GravitySolverBase::explicitStep( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) +{ + switch( m_mode ) + { + case GravityMode::Modeling: + return explicitStepModeling( time_n, dt, cycleNumber, domain ); + case GravityMode::Adjoint: + return explicitStepAdjoint( time_n, dt, cycleNumber, domain ); + default: + GEOS_THROW( "GravitySolverBase::explicitStep: Mode not supported", InputError ); + } +} + +void GravitySolverBase::saveGz( real64 const & time_n, + integer const cycleNumber, + string const basename, + arrayView1d< real64 > const & gzAtStations ) +{ + // Convert to float32 + std::vector< float > tmp32( gzAtStations.size()); + std::transform( gzAtStations.begin(), gzAtStations.end(), tmp32.begin(), + [] ( real64 val )-> float { return static_cast< float >(val); } ); + + // Binary data + std::filesystem::path dataFilename = basename + GEOS_FMT( "_{:015}.H@", time_n ); + std::ofstream fout( dataFilename, std::ios::binary ); + if( !fout ) + throw std::ios_base::failure( "Failed to open data file" ); + fout.write( reinterpret_cast< const char * >(tmp32.data()), tmp32.size() * sizeof(float)); + + // Header file + const std::string headerFilename = basename + GEOS_FMT( "_{:015}.H", time_n ); + std::ofstream fheader( headerFilename ); + if( !fheader ) + throw std::ios_base::failure( "Failed to open header file" ); + + fheader << "n1=" << gzAtStations.size() << '\n' + << "d1=1.\n" + << "o1=0.\n" + << "label1=STATION\n" + << "esize=4\n" + << "data_format=native_float\n" + << "data_label=GZ\n" + << "in=" << dataFilename << '\n' + << GEOS_FMT( "time={:015}", time_n ) << '\n' + << "cycle=" << cycleNumber << '\n'; +} + +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp new file mode 100644 index 00000000000..0eaee0e6065 --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp @@ -0,0 +1,137 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 GravitySolverBase.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYSOLVERBASE_HPP_ +#define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYSOLVERBASE_HPP_ + +#include "mesh/MeshFields.hpp" +#include "physicsSolvers/PhysicsSolverBase.hpp" + + +namespace geos +{ + +class GravitySolverBase : public PhysicsSolverBase +{ +public: + + /// The default nullary constructor is disabled to avoid compiler auto-generation: + GravitySolverBase() = delete; + + /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) + GravitySolverBase( const std::string & name, + Group * const parent ); + + virtual ~GravitySolverBase() override; + + GravitySolverBase( GravitySolverBase const & ) = delete; + GravitySolverBase( GravitySolverBase && ) = default; + + GravitySolverBase & operator=( GravitySolverBase const & ) = delete; + GravitySolverBase & operator=( GravitySolverBase && ) = delete; + + virtual void postInputInitialization() override; + + virtual void initializePreSubGroups() override; + virtual void initializePostInitialConditionsPreSubGroups() override; + + + virtual real64 solverStep( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override; + + + virtual real64 explicitStep( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override; + + + /// Bind data between input XML file and source code + struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct + { + static constexpr char const * modeString() { return "mode"; } + static constexpr char const * stationCoordinatesString() { return "stationCoordinates"; } + static constexpr char const * outputPropertyLogString() { return "outputPropertyLog"; } + static constexpr char const * outputGzLogString() { return "outputGzLog"; } + static constexpr char const * outputGzString() { return "outputGz"; } + static constexpr char const * outputGzBasenameString() { return "outputGzBasename"; } + static constexpr char const * outputAdjointLogString() { return "outputAdjointLog"; } + static constexpr char const * residueString() { return "residue"; } + static constexpr char const * gzAtStationsString() { return "gzAtStations"; } + }; + + void reinit() override final; + + +protected: + + void saveGz( real64 const & time_n, + integer const cycleNumber, + string const basename, + arrayView1d< real64 > const & gzAtStations ); + + virtual real64 explicitStepModeling( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) = 0; + + virtual real64 explicitStepAdjoint( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) = 0; + + + enum class GravityMode { Modeling, Adjoint }; + static const std::unordered_map< std::string, GravityMode > modeMap; + + + GravityMode parseMode( const std::string & modeStr ) + { + auto it = modeMap.find( modeStr ); + if( it != modeMap.end()) + return it->second; + throw std::invalid_argument( "Invalid mode string: " + modeStr ); + } + + + + std::string m_modeString; + + GravityMode m_mode; + + + /// Coordinates of the gravimeter stations + array2d< real64 > m_stationCoordinates; + + /// Dump vertical component Gz to disk + localIndex m_outputGz; + string m_outputGzBasename; + + /// Residue observed at station (only for adjoint computation) + array1d< real64 > m_residue; + + /// Gz component recorded at stations + array1d< real64 > m_gzAtStations; +}; + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYSOLVERBASE_HPP_ From aea99ebdf62570a332038a0e4c01b9fd493fb360 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Fri, 27 Jun 2025 12:51:19 -0500 Subject: [PATCH 02/30] Clean up. --- .../physicsSolvers/gravity/GravityFE.cpp | 45 ++------------ .../physicsSolvers/gravity/GravityFE.hpp | 9 +-- .../gravity/GravityFEKernel.hpp | 12 ++-- .../gravity/GravitySolverBase.cpp | 61 ++++++++----------- .../gravity/GravitySolverBase.hpp | 43 +++++-------- 5 files changed, 52 insertions(+), 118 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 045df4cfb5f..6348eedac1f 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -17,18 +17,12 @@ /** * @file GravityFE.cpp */ -#include #include "GravityFE.hpp" #include "GravityFEKernel.hpp" -#include "discretizationMethods/NumericalMethodsManager.hpp" #include "finiteElement/FiniteElementDiscretization.hpp" -#include "fieldSpecification/FieldSpecificationManager.hpp" -#include "mainInterface/ProblemManager.hpp" -#include "mesh/ElementType.hpp" #include "mesh/DomainPartition.hpp" -#include "mesh/mpiCommunications/CommunicationTools.hpp" namespace geos @@ -70,7 +64,6 @@ void GravityFE::initializePreSubGroups() void GravityFE::registerDataOnMesh( Group & meshBodies ) { - forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, MeshLevel & mesh, string_array const & ) @@ -131,16 +124,6 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, arrayView1d< real64 > const volumeIntegral = nodeManager.getField< fields::VolumeIntegral >(); volumeIntegral.setValues< parallelHostPolicy >( 0. ); - -#if 0 - // Print size and type of "CellElements" partitions. - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & elementSubRegion ) - { - GEOS_LOG_RANK_0( "GravityFE: Rank 0: Partition[ " << elementSubRegion.getElementType() << " ] = "<< elementSubRegion.size()); - } ); -#endif - // Loop over all sub-regions in regions of type "CellElements". mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & elementSubRegion ) @@ -166,7 +149,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, using FE_TYPE = TYPEOFREF( finiteElement ); gravityFEKernel:: - DensityVolumeIntegralKernel< FE_TYPE > kernel( finiteElement ); + ForwardVolumeIntegralKernel< FE_TYPE > kernel( finiteElement ); kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > ( elementSubRegion.size(), X, @@ -183,7 +166,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, volumeIntegral[a]=0.; } } ); - } ); // loop on cellElements + } ); // loop cellElements // Step #2: Compute contribution to all stations. @@ -218,17 +201,8 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, // In place: Source==Destination. MpiWrapper::allReduce( localGzAtStations, localGzAtStations, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); -#if 0 - arrayView1d< real64 > const gzAtStations = m_gzAtStations.toView(); - for( localIndex istation=0; istation gzAtStations = m_gzAtStations.toView(); gzAtStations = localGzAtStations; -#endif - for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) { @@ -244,8 +218,6 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, GEOS_LOG_RANK_0( logStream.str() ); } - - // Dump result to disk... if( (rank==0) && ( this->m_outputGz == 1 )) { @@ -264,8 +236,6 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( time_n, cycleNumber ); - //int const rank = MpiWrapper::commRank( MPI_COMM_GEOSX ); - array1d< real64 > localGzAtStations( m_stationCoordinates.size( 0 ) ); localGzAtStations.setValues< parallelHostPolicy >( 0. ); @@ -304,7 +274,7 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, using FE_TYPE = TYPEOFREF( finiteElement ); gravityFEKernel:: - VolumeIntegralKernel_uni2< FE_TYPE > kernel( finiteElement ); + AdjointVolumeIntegralKernel< FE_TYPE > kernel( finiteElement ); kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > ( elementSubRegion.size(), X, @@ -317,11 +287,6 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, { // Deal with one station. auto const & coords = m_stationCoordinates[iStation]; -#if 0 - real64 const coords[3] = { m_stationCoordinates[iStation][0], - m_stationCoordinates[iStation][1], - m_stationCoordinates[iStation][2] }; -#endif real64 const res=residue[iStation]; forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) @@ -336,9 +301,7 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, real64 r2 = dx*dx + dy*dy + dz*dz; real64 r3 = sqrt( r2 ) * r2; adjoint[k] += GRAVITATIONAL_CONSTANT * volumeIntegral2d[k][iLoc] * res * dz / r3; - } - } ); // Loop elem } //Loop station @@ -352,7 +315,7 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, } } - } );// Loop subregion + } ); // Loop subregion } ); // Loop mesh return dt; diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp index d4759a85272..2cfcf2db090 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp @@ -21,9 +21,8 @@ #ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ #define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ -#include "physicsSolvers/PhysicsSolverBase.hpp" #include "GravitySolverBase.hpp" - +#include "mesh/MeshFields.hpp" namespace geos @@ -86,14 +85,9 @@ class GravityFE : public GravitySolverBase virtual void initializePostInitialConditionsPreSubGroups() override final; - -private: - - }; - namespace fields { DECLARE_FIELD( MediumDensity, @@ -130,7 +124,6 @@ DECLARE_FIELD( Adjoint, } - } /* namespace geos */ #endif /* GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp index 3c41586f7de..4438e3a9ead 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp @@ -21,7 +21,6 @@ #ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ #define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ -//#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" namespace geos { @@ -31,9 +30,9 @@ namespace gravityFEKernel { template< typename FE_TYPE > -struct DensityVolumeIntegralKernel +struct ForwardVolumeIntegralKernel { - explicit DensityVolumeIntegralKernel( FE_TYPE const & finiteElement ) + explicit ForwardVolumeIntegralKernel( FE_TYPE const & finiteElement ) : m_finiteElement( finiteElement ) {} @@ -98,9 +97,9 @@ struct DensityVolumeIntegralKernel template< typename FE_TYPE > -struct VolumeIntegralKernel_uni2 +struct AdjointVolumeIntegralKernel { - explicit VolumeIntegralKernel_uni2( FE_TYPE const & finiteElement ) + explicit AdjointVolumeIntegralKernel( FE_TYPE const & finiteElement ) : m_finiteElement( finiteElement ) {} @@ -156,10 +155,7 @@ struct VolumeIntegralKernel_uni2 FE_TYPE const & m_finiteElement; }; - - } // namespace gravityFEKernel - } // namespace geos #endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp index a2867c5c9ab..cbd2f8e644e 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp @@ -18,27 +18,14 @@ */ #include "GravitySolverBase.hpp" - -#include "dataRepository/KeyNames.hpp" -#include "finiteElement/FiniteElementDiscretization.hpp" -#include "fieldSpecification/FieldSpecificationManager.hpp" -#include "mainInterface/ProblemManager.hpp" -#include "mesh/mpiCommunications/CommunicationTools.hpp" #include + namespace geos { using namespace dataRepository; - - -const std::unordered_map< std::string, GravitySolverBase::GravityMode > GravitySolverBase::modeMap = { - { "modeling", GravitySolverBase::GravityMode::Modeling }, - { "adjoint", GravitySolverBase::GravityMode::Adjoint } -}; - - GravitySolverBase::GravitySolverBase( const std::string & name, Group * const parent ) : PhysicsSolverBase( name, parent ) { @@ -75,6 +62,31 @@ GravitySolverBase::GravitySolverBase( const std::string & name, Group * const pa GravitySolverBase::~GravitySolverBase() = default; + +GravitySolverBase::GravityMode GravitySolverBase::parseMode( const std::string & modeStr ) +{ + std::string lowerMode = modeStr; + std::transform( lowerMode.begin(), lowerMode.end(), lowerMode.begin(), + [] ( unsigned char c ) -> unsigned char { return std::tolower( c ); } ); + + auto it = modeMap.find( lowerMode ); + if( it != modeMap.end()) + return it->second; + + std::ostringstream oss; + oss << "Invalid gravity mode string: '" << modeStr << "'. Valid options are: "; + bool first = true; + for( const auto & pair : modeMap ) + { + if( !first ) + oss << ", "; + oss << pair.first; + first = false; + } + throw std::invalid_argument( oss.str()); +} + + void GravitySolverBase::reinit() { initializePostInitialConditionsPreSubGroups(); @@ -89,7 +101,6 @@ void GravitySolverBase::initializePostInitialConditionsPreSubGroups() { PhysicsSolverBase::initializePostInitialConditionsPreSubGroups(); - try { m_mode = parseMode( m_modeString ); @@ -107,24 +118,6 @@ void GravitySolverBase::initializePostInitialConditionsPreSubGroups() } - -#if 0 - if( m_modeString == "modeling" ) - { - m_mode = GravityMode_Modeling; - } - else if( m_modeString == "adjoint" ) - { - m_mode = GravityMode_Adjoint; - GEOS_THROW_IF( m_stationCoordinates.size( 0 ) != m_residue.size( 0 ), - "GravitySolverBase: Residue size does not match the number of stations", - InputError ); - } - else - { - GEOS_THROW( "Invalid mode string: " + m_modeString, InputError ); - } -#endif constexpr auto yesno = [] (int flag) noexcept->const char * { return flag ? "yes" : "no"; }; LogPart gravitySolverLog( "Gravity Solver: ", MpiWrapper::commRank() == 0 ); @@ -195,7 +188,7 @@ void GravitySolverBase::saveGz( real64 const & time_n, fout.write( reinterpret_cast< const char * >(tmp32.data()), tmp32.size() * sizeof(float)); // Header file - const std::string headerFilename = basename + GEOS_FMT( "_{:015}.H", time_n ); + std::filesystem::path headerFilename = basename + GEOS_FMT( "_{:015}.H", time_n ); std::ofstream fheader( headerFilename ); if( !fheader ) throw std::ios_base::failure( "Failed to open header file" ); diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp index 0eaee0e6065..478e94e95ab 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp @@ -20,7 +20,6 @@ #ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYSOLVERBASE_HPP_ #define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYSOLVERBASE_HPP_ -#include "mesh/MeshFields.hpp" #include "physicsSolvers/PhysicsSolverBase.hpp" @@ -69,11 +68,8 @@ class GravitySolverBase : public PhysicsSolverBase { static constexpr char const * modeString() { return "mode"; } static constexpr char const * stationCoordinatesString() { return "stationCoordinates"; } - static constexpr char const * outputPropertyLogString() { return "outputPropertyLog"; } - static constexpr char const * outputGzLogString() { return "outputGzLog"; } static constexpr char const * outputGzString() { return "outputGz"; } static constexpr char const * outputGzBasenameString() { return "outputGzBasename"; } - static constexpr char const * outputAdjointLogString() { return "outputAdjointLog"; } static constexpr char const * residueString() { return "residue"; } static constexpr char const * gzAtStationsString() { return "gzAtStations"; } }; @@ -83,10 +79,18 @@ class GravitySolverBase : public PhysicsSolverBase protected: - void saveGz( real64 const & time_n, - integer const cycleNumber, - string const basename, - arrayView1d< real64 > const & gzAtStations ); + enum class GravityMode { Modeling, Adjoint }; + + inline static const std::unordered_map< std::string, GravityMode > modeMap = { + {"modeling", GravityMode::Modeling}, + {"adjoint", GravityMode::Adjoint} + }; + + static GravityMode parseMode( const std::string & modeStr ); + + std::string m_modeString; + GravityMode m_mode; + virtual real64 explicitStepModeling( real64 const & time_n, real64 const & dt, @@ -98,25 +102,10 @@ class GravitySolverBase : public PhysicsSolverBase integer const cycleNumber, DomainPartition & domain ) = 0; - - enum class GravityMode { Modeling, Adjoint }; - static const std::unordered_map< std::string, GravityMode > modeMap; - - - GravityMode parseMode( const std::string & modeStr ) - { - auto it = modeMap.find( modeStr ); - if( it != modeMap.end()) - return it->second; - throw std::invalid_argument( "Invalid mode string: " + modeStr ); - } - - - - std::string m_modeString; - - GravityMode m_mode; - + void saveGz( real64 const & time_n, + integer const cycleNumber, + string const basename, + arrayView1d< real64 > const & gzAtStations ); /// Coordinates of the gravimeter stations array2d< real64 > m_stationCoordinates; From d42d7bfa6813589a19b683ac096d5b1a90accf7c Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Mon, 30 Jun 2025 12:47:22 -0500 Subject: [PATCH 03/30] Fixed station output. --- .../physicsSolvers/gravity/GravityFE.cpp | 30 +++++++++++-------- 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 6348eedac1f..0800e67d3d4 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -202,20 +202,26 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, MpiWrapper::allReduce( localGzAtStations, localGzAtStations, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); arrayView1d< real64 > gzAtStations = m_gzAtStations.toView(); - gzAtStations = localGzAtStations; + for( localIndex i = 0; i < gzAtStations.size(); ++i ) + { + gzAtStations[i] = localGzAtStations[i]; + } - for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) + if( this->getLogLevel()>1 ) { - auto const & coords = m_stationCoordinates[iStation]; - std::ostringstream logStream; - logStream << std::fixed << std::setprecision( 2 ); - logStream << "GravityFE: station[" << std::setw( 5 ) << iStation << "] " - << std::setw( 15 ) << coords[0] << " " - << std::setw( 15 ) << coords[1] << " " - << std::setw( 10 ) << coords[2] << " " - << std::scientific << std::setprecision( 6 ) - << std::setw( 14 ) << gzAtStations[iStation]; - GEOS_LOG_RANK_0( logStream.str() ); + for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) + { + auto const & coords = m_stationCoordinates[iStation]; + std::ostringstream logStream; + logStream << std::fixed << std::setprecision( 2 ); + logStream << "GravityFE: station[" << std::setw( 5 ) << iStation << "] " + << std::setw( 15 ) << coords[0] << " " + << std::setw( 15 ) << coords[1] << " " + << std::setw( 10 ) << coords[2] << " " + << std::scientific << std::setprecision( 6 ) + << std::setw( 14 ) << gzAtStations[iStation]; + GEOS_LOG_RANK_0( logStream.str() ); + } } // Dump result to disk... From 2ab1cafe638ceb0d526f277045afef72ed6c068b Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Tue, 1 Jul 2025 11:37:12 -0500 Subject: [PATCH 04/30] Draft. --- .../dataRepository/python/PyWrapper.cpp | 85 +++++++++++-- .../dataRepository/python/PyWrapper.hpp | 16 +++ .../dataRepository/unitTests/CMakeLists.txt | 5 + .../unitTests/testPyWrapper_setValue.cpp | 116 ++++++++++++++++++ 4 files changed, 211 insertions(+), 11 deletions(-) create mode 100644 src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp diff --git a/src/coreComponents/dataRepository/python/PyWrapper.cpp b/src/coreComponents/dataRepository/python/PyWrapper.cpp index f90c4e01715..f87d342bffe 100644 --- a/src/coreComponents/dataRepository/python/PyWrapper.cpp +++ b/src/coreComponents/dataRepository/python/PyWrapper.cpp @@ -19,7 +19,7 @@ // Source includes #include "PyWrapper.hpp" - +#include "dataRepository/Wrapper.hpp" #define VERIFY_NON_NULL_SELF( self ) \ PYTHON_ERROR_IF( self == nullptr, PyExc_RuntimeError, "Passed a nullptr as self.", nullptr ) @@ -32,16 +32,6 @@ namespace geos namespace python { -struct PyWrapper -{ - PyObject_HEAD - - static constexpr char const * docString = - "A Python interface to geos::dataRepository::WrapperBase."; - - dataRepository::WrapperBase * wrapper; -}; - /** * */ @@ -100,10 +90,83 @@ static PyObject * PyWrapper_value( PyWrapper * const self, PyObject * const args return ret; } + +template< typename T > +bool trySetValue( dataRepository::WrapperBase * base, T const & value ) +{ + auto * wrapper = dynamic_cast< dataRepository::Wrapper< T > * >( base ); + if( wrapper ) + { + wrapper->reference() = value; + return true; + } + return false; +} + + +PyObject * PyWrapper_setValue( PyWrapper * const self, PyObject * const args ) +{ + VERIFY_NON_NULL_SELF( self ); + VERIFY_INITIALIZED( self ); + + PyObject * pyValue = nullptr; + if( !PyArg_ParseTuple( args, "O", &pyValue ) ) + { + return nullptr; + } + + try + { + bool success = false; + + if( PyUnicode_Check( pyValue ) ) + { + const char * strVal = PyUnicode_AsUTF8( pyValue ); + success = trySetValue( self->wrapper, std::string( strVal ) ); + if( !success ) + PyErr_SetString( PyExc_TypeError, "Wrapper is not a std::string." ); + } + else if( PyFloat_Check( pyValue ) ) + { + double val = PyFloat_AsDouble( pyValue ); + success = trySetValue( self->wrapper, val ); + if( !success ) + PyErr_SetString( PyExc_TypeError, "Wrapper is not a double." ); + } + else if( PyLong_Check( pyValue ) ) + { + long val = PyLong_AsLong( pyValue ); + success = trySetValue( self->wrapper, static_cast< int >( val ) ) || + trySetValue( self->wrapper, val ); + if( !success ) + PyErr_SetString( PyExc_TypeError, "Wrapper is not an int or long." ); + } + else + { + PyErr_SetString( PyExc_TypeError, "Unsupported Python type for assignment." ); + return nullptr; + } + + if( !success ) + { + return nullptr; + } + + Py_RETURN_NONE; + } + catch( std::exception const & e ) + { + PyErr_SetString( PyExc_RuntimeError, e.what() ); + return nullptr; + } +} + + BEGIN_ALLOW_DESIGNATED_INITIALIZERS static PyMethodDef PyWrapperMethods[] = { { "value", (PyCFunction) PyWrapper_value, METH_VARARGS, PyWrapper_valueDocString }, + { "set_value", (PyCFunction) PyWrapper_setValue, METH_VARARGS, "Set the value of a scalar wrapper (string, int, float)." }, { nullptr, nullptr, 0, nullptr } // Sentinel }; diff --git a/src/coreComponents/dataRepository/python/PyWrapper.hpp b/src/coreComponents/dataRepository/python/PyWrapper.hpp index c059459e397..aad6414dbb1 100644 --- a/src/coreComponents/dataRepository/python/PyWrapper.hpp +++ b/src/coreComponents/dataRepository/python/PyWrapper.hpp @@ -24,6 +24,18 @@ namespace geos namespace python { + +struct PyWrapper +{ + PyObject_HEAD + + static constexpr char const * docString = + "A Python interface to geos::dataRepository::WrapperBase."; + + dataRepository::WrapperBase * wrapper; +}; + + /** * */ @@ -34,6 +46,10 @@ PyObject * createNewPyWrapper( dataRepository::WrapperBase & wrapper ); */ PyTypeObject * getPyWrapperType(); + +PyObject * PyWrapper_setValue( PyWrapper * self, PyObject * args ); + + } // namespace python } // namespace geos diff --git a/src/coreComponents/dataRepository/unitTests/CMakeLists.txt b/src/coreComponents/dataRepository/unitTests/CMakeLists.txt index eb42386660e..34fb3f41b84 100644 --- a/src/coreComponents/dataRepository/unitTests/CMakeLists.txt +++ b/src/coreComponents/dataRepository/unitTests/CMakeLists.txt @@ -7,6 +7,11 @@ set( dataRepository_tests testXmlWrapper.cpp testBufferOps.cpp ) + +if( ENABLE_PYGEOSX ) + list( APPEND dataRepository_tests testPyWrapper_setValue.cpp ) +endif() + set( dependencyList ${parallelDeps} gtest dataRepository ) if( ENABLE_CUDA AND ENABLE_CUDA_NVTOOLSEXT ) diff --git a/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp b/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp new file mode 100644 index 00000000000..aec76e978a7 --- /dev/null +++ b/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp @@ -0,0 +1,116 @@ +#include +#include +#include "dataRepository/Group.hpp" +#include "dataRepository/Wrapper.hpp" + +#include "dataRepository/python/PyWrapper.hpp" + + + +using namespace geos::python; +using namespace geos::dataRepository; + + +class PythonEnvironment : public ::testing::Environment +{ +public: + void SetUp() override { Py_Initialize(); } + void TearDown() override { Py_Finalize(); } +}; + +::testing::Environment * const pythonEnv = ::testing::AddGlobalTestEnvironment( new PythonEnvironment ); + + +template< typename T > +void testSetValue( PyObject * arg, const T & expected, const std::string & typeName ) +{ + conduit::Node node; + Group rootGroup( "root", node ); + Group dummyGroup( "dummy", &rootGroup ); + T * value = new T( expected ); + auto * wrapper = new Wrapper< T >( "test", dummyGroup, value ); + + PyWrapper pyWrapper; + pyWrapper.wrapper = wrapper; + + PyObject * result = PyWrapper_setValue( &pyWrapper, arg ); + ASSERT_NE( result, nullptr ) << "Failed to set value for type: " << typeName; + + EXPECT_EQ( wrapper->reference(), expected ); + + delete wrapper; + Py_DECREF( arg ); + Py_XDECREF( result ); +} + + +TEST( PyWrapperSetValue, SetString ) +{ + PyObject * arg = Py_BuildValue( "(s)", "hello" ); + testSetValue< std::string >( arg, "hello", "string" ); +} + +TEST( PyWrapperSetValue, SetDouble ) +{ + PyObject * arg = Py_BuildValue( "(d)", 3.14 ); + testSetValue< double >( arg, 3.14, "double" ); +} + +TEST( PyWrapperSetValue, SetInt ) +{ + PyObject * arg = Py_BuildValue( "(i)", 42 ); + testSetValue< int >( arg, 42, "int" ); +} + +TEST( PyWrapperSetValue, SetLong ) +{ + PyObject * arg = Py_BuildValue( "(l)", 123456789L ); + testSetValue< long >( arg, 123456789L, "long" ); +} + + +TEST( PyWrapperSetValue, TypeMismatch ) +{ + conduit::Node node; + Group rootGroup( "root", node ); + Group dummyGroup( "dummy", &rootGroup ); + int * value = new int(0); + auto * wrapper = new Wrapper< int >( "test", dummyGroup, value ); + + PyWrapper pyWrapper; + pyWrapper.wrapper = wrapper; + + PyObject * arg = Py_BuildValue( "(s)", "not an int" ); + PyObject * result = PyWrapper_setValue( &pyWrapper, arg ); + + EXPECT_EQ( result, nullptr ); + EXPECT_TRUE( PyErr_Occurred() ); + PyErr_Clear(); + + delete wrapper; + Py_DECREF( arg ); +} + + +TEST( PyWrapperSetValue, UnsupportedType ) +{ + conduit::Node node; + Group rootGroup( "root", node ); + Group dummyGroup( "dummy", &rootGroup ); + auto * wrapper = new Wrapper< std::vector< int > >( "test", dummyGroup, new std::vector< int >() ); + + PyWrapper pyWrapper; + pyWrapper.wrapper = wrapper; + + PyObject * dict = PyDict_New(); + PyObject * arg = Py_BuildValue( "(O)", dict ); + PyObject * result = PyWrapper_setValue( &pyWrapper, arg ); + + EXPECT_EQ( result, nullptr ); + EXPECT_TRUE( PyErr_Occurred()); + PyErr_Clear(); + + delete wrapper; + Py_DECREF( dict ); + Py_DECREF( arg ); +} From be4753050bc2d4a2f07e1255847af859f942f2c8 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Tue, 1 Jul 2025 12:17:32 -0500 Subject: [PATCH 05/30] Added doc. --- .../dataRepository/python/PyWrapper.hpp | 26 ++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/dataRepository/python/PyWrapper.hpp b/src/coreComponents/dataRepository/python/PyWrapper.hpp index aad6414dbb1..d9fcf73fdcc 100644 --- a/src/coreComponents/dataRepository/python/PyWrapper.hpp +++ b/src/coreComponents/dataRepository/python/PyWrapper.hpp @@ -25,13 +25,21 @@ namespace python { +/** + * @brief Python wrapper object for geos::dataRepository::WrapperBase. + * + * This struct defines the Python-facing interface for accessing and modifying + * wrapped C++ values from Python. + */ struct PyWrapper { PyObject_HEAD + /// Documentation string for the PyWrapper Python type. static constexpr char const * docString = "A Python interface to geos::dataRepository::WrapperBase."; + /// Pointer to the underlying C++ WrapperBase instance. dataRepository::WrapperBase * wrapper; }; @@ -46,7 +54,23 @@ PyObject * createNewPyWrapper( dataRepository::WrapperBase & wrapper ); */ PyTypeObject * getPyWrapperType(); - +/** + * @brief Set the value of a scalar wrapper from Python. + * + * This function allows assigning a scalar value from Python into the underlying + * C++ Wrapper instance. Supported types include: + * - std::string + * - int + * - long + * - double + * + * If the type of the Python object does not match the wrapped C++ type, + * a Python TypeError is raised. Same for unsupported types. + * + * @param self Pointer to the PyWrapper instance. + * @param args A Python tuple containing a single value to assign. + * @return Py_None on success, or nullptr on failure with a Python exception set. + */ PyObject * PyWrapper_setValue( PyWrapper * self, PyObject * args ); From 2cc94b8a12bfbac7b57436bfa03625a3f611de52 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Tue, 1 Jul 2025 12:32:21 -0500 Subject: [PATCH 06/30] Added copyrights. --- .../unitTests/testPyWrapper_setValue.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp b/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp index aec76e978a7..aef65a795ea 100644 --- a/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp +++ b/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp @@ -1,3 +1,18 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + #include #include #include "dataRepository/Group.hpp" From 76b439aea06be2a4caedca7ffbfc2917127e0538 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 10 Jul 2025 11:39:12 -0500 Subject: [PATCH 07/30] Added unit tests. --- .../physicsSolvers/gravity/GravityFE.cpp | 1 - .../physicsSolvers/gravity/GravityFE.hpp | 3 + .../gravity/GravityFEKernel.hpp | 2 + .../gravity/GravitySolverBase.cpp | 45 ++-- src/coreComponents/unitTests/CMakeLists.txt | 4 + .../unitTests/gravityTests/CMakeLists.txt | 39 ++++ .../gravityTests/testGravitySolver.cpp | 197 +++++++++++++++++ .../gravityTests/testGravitySolverAdjoint.cpp | 202 ++++++++++++++++++ 8 files changed, 475 insertions(+), 18 deletions(-) create mode 100644 src/coreComponents/unitTests/gravityTests/CMakeLists.txt create mode 100644 src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp create mode 100644 src/coreComponents/unitTests/gravityTests/testGravitySolverAdjoint.cpp diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 0800e67d3d4..881bf006849 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -31,7 +31,6 @@ namespace geos using namespace constitutive; using namespace dataRepository; -constexpr real64 GRAVITATIONAL_CONSTANT = 6.67430e-11; // in m3.kg-1.s-2 // Older value: 6.6738480e-11 constexpr localIndex MAX_SUPPORT_POINTS = 8; diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp index 2cfcf2db090..85d9bdf1833 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp @@ -28,6 +28,9 @@ namespace geos { +constexpr real64 GRAVITATIONAL_CONSTANT = 6.67430e-11; // in m³·kg⁻¹·s⁻² + + class GravityFE : public GravitySolverBase { public: diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp index 4438e3a9ead..bf2db96bd07 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp @@ -21,6 +21,8 @@ #ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ #define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFEKERNEL_HPP_ +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" namespace geos { diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp index cbd2f8e644e..5989ef81c31 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp @@ -112,26 +112,37 @@ void GravitySolverBase::initializePostInitialConditionsPreSubGroups() if( m_mode == GravityMode::Adjoint ) { - GEOS_THROW_IF( m_stationCoordinates.size( 0 ) != m_residue.size( 0 ), - "GravitySolverBase: Residue size does not match the number of stations", - InputError ); + + + const auto stationCount = m_stationCoordinates.size( 0 ); + const auto residueCount = m_residue.size( 0 ); + + GEOS_THROW_IF( + stationCount != residueCount, + "GravitySolverBase: Residue size (" + std::to_string( residueCount ) + + ") does not match the number of stations (" + std::to_string( stationCount ) + ")", + InputError + ); } - constexpr auto yesno = [] (int flag) noexcept->const char * { return flag ? "yes" : "no"; }; - - LogPart gravitySolverLog( "Gravity Solver: ", MpiWrapper::commRank() == 0 ); - gravitySolverLog.begin(); - GEOS_LOG_RANK_0( "Name: " << getName()); - GEOS_LOG_RANK_0( "Mode: " << m_modeString ); - GEOS_LOG_RANK_0( "Number of stations: " << m_stationCoordinates.size( 0 )); - GEOS_LOG_RANK_0( "Output Gz to file: " << yesno( m_outputGz )); - GEOS_LOG_RANK_0( "Output Gz basename: " << m_outputGzBasename ); - GEOS_LOG_RANK_0( "Log level: " << getLogLevel()); - GEOS_LOG_RANK_0( " Output Gz to logs: " << yesno( getLogLevel() > 1 )); - GEOS_LOG_RANK_0( " Output Adjoint to logs: " << yesno( getLogLevel() > 2 )); - GEOS_LOG_RANK_0( " Output Properties to logs: " << yesno( getLogLevel() > 3 )); - gravitySolverLog.end(); + if( this->getLogLevel()>0 ) + { + constexpr auto yesno = [] (int flag) noexcept->const char * { return flag ? "yes" : "no"; }; + + LogPart gravitySolverLog( "Gravity Solver: ", MpiWrapper::commRank() == 0 ); + gravitySolverLog.begin(); + GEOS_LOG_RANK_0( "Name: " << getName()); + GEOS_LOG_RANK_0( "Mode: " << m_modeString ); + GEOS_LOG_RANK_0( "Number of stations: " << m_stationCoordinates.size( 0 )); + GEOS_LOG_RANK_0( "Output Gz to file: " << yesno( m_outputGz )); + GEOS_LOG_RANK_0( "Output Gz basename: " << m_outputGzBasename ); + GEOS_LOG_RANK_0( "Log level: " << getLogLevel()); + GEOS_LOG_RANK_0( " Output Gz to logs: " << yesno( getLogLevel() > 1 )); + GEOS_LOG_RANK_0( " Output Adjoint to logs: " << yesno( getLogLevel() > 2 )); + GEOS_LOG_RANK_0( " Output Properties to logs: " << yesno( getLogLevel() > 3 )); + gravitySolverLog.end(); + } } void GravitySolverBase::postInputInitialization() diff --git a/src/coreComponents/unitTests/CMakeLists.txt b/src/coreComponents/unitTests/CMakeLists.txt index 8b760067a42..ef38b36a9b5 100644 --- a/src/coreComponents/unitTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/CMakeLists.txt @@ -34,3 +34,7 @@ if( GEOS_ENABLE_FLUIDFLOW ) endif() add_subdirectory( wellsTests ) add_subdirectory( wavePropagationTests ) +if( GEOS_ENABLE_GRAVITY ) + add_subdirectory( gravityTests ) +endif() + diff --git a/src/coreComponents/unitTests/gravityTests/CMakeLists.txt b/src/coreComponents/unitTests/gravityTests/CMakeLists.txt new file mode 100644 index 00000000000..629d5307b08 --- /dev/null +++ b/src/coreComponents/unitTests/gravityTests/CMakeLists.txt @@ -0,0 +1,39 @@ +# Specify list of tests +set( gtest_geosx_tests + testGravitySolver.cpp + testGravitySolverAdjoint.cpp + ) + +set( tplDependencyList ${parallelDeps} gtest ) + +set( dependencyList mainInterface ) + +geos_decorate_link_dependencies( LIST decoratedDependencies + DEPENDENCIES ${dependencyList} ) + +# Add gtest C++ based tests +foreach(test ${gtest_geosx_tests}) + get_filename_component( test_name ${test} NAME_WE ) + blt_add_executable( NAME ${test_name} + SOURCES ${test} + OUTPUT_DIR ${TEST_OUTPUT_DIRECTORY} + DEPENDS_ON ${decoratedDependencies} ${tplDependencyList} ) + + # Guard to prevent GCC (version 8) from giving warnings due + # to some sort of possible conversion from int to long unsigned. + # See also discussion here: https://github.com/GEOS-DEV/LvArray/pull/250 + if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 8.3.1) + target_compile_options(${test_name} PRIVATE "-Wno-stringop-overflow") + endif() + endif() + + geos_add_test( NAME ${test_name} + COMMAND ${test_name} ) + +endforeach() + +# For some reason, BLT is not setting CUDA language for these source files +if ( ENABLE_CUDA ) + set_source_files_properties( ${gtest_geosx_tests} PROPERTIES LANGUAGE CUDA ) +endif() diff --git a/src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp b/src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp new file mode 100644 index 00000000000..01b00f051f0 --- /dev/null +++ b/src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp @@ -0,0 +1,197 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include + +#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" // For setupProblemFromXML +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/gravity/GravityFE.hpp" + +using namespace geos; +using namespace geos::dataRepository; +using namespace geos::testing; + +CommandLineOptions g_commandLineOptions; + + +// Shared constants +constexpr std::array< real64, 6 > PRISM_BOUNDS = { 10000, 14000, 8000, 11000, 2000, 2200 }; +constexpr real64 DX = 500.0; +constexpr real64 DY = 500.0; +constexpr real64 Z_STATION = 10.0; +constexpr localIndex NX = 41; +constexpr localIndex NY = 31; +constexpr real64 DENSITY_CONTRAST = 900.0; + + +/** + * @brief Analytic solution for gravity anomaly (gz) due to a rectangular prism. + * + * @param station Coordinates of the observation point (x, y, z) + * @param prism_bounds Bounds of the prism: {x0, x1, y0, y1, z0, z1} + * @param density_contrast Density contrast of the prism (kg/m^3) + * @return double Gravity anomaly in m/s^2 + */ +double analyticPrismGz( const std::array< double, 3 > & station, + const std::array< double, 6 > & prism_bounds, + double density_contrast ) +{ + const double x = station[0]; + const double y = station[1]; + const double z = station[2]; + + const double x0 = prism_bounds[0]; + const double x1 = prism_bounds[1]; + const double y0 = prism_bounds[2]; + const double y1 = prism_bounds[3]; + const double z0 = prism_bounds[4]; + const double z1 = prism_bounds[5]; + + const double epsilon = 1e-16; + double val = 0.0; + + for( int i = 0; i < 2; ++i ) + { + for( int j = 0; j < 2; ++j ) + { + for( int k = 0; k < 2; ++k ) + { + double xi = (i == 0) ? x0 : x1; + double yj = (j == 0) ? y0 : y1; + double zk = (k == 0) ? z0 : z1; + + double dx = xi - x; + double dy = yj - y; + double dz = zk - z; + + double r = std::sqrt( dx * dx + dy * dy + dz * dz ); + + if( r > epsilon ) + { + double sign = ((i + j + k) % 2 == 0) ? 1.0 : -1.0; + double arg1 = std::atan2( dx * dy, dz * r + epsilon ); + double arg2 = std::log( r + dy + epsilon ); + double arg3 = std::log( r + dx + epsilon ); + + val += sign * (dz * arg1 - dx * arg2 - dy * arg3); + } + } + } + } + + return -GRAVITATIONAL_CONSTANT * density_contrast * val; +} + + + +class GravitySolverTest : public ::testing::Test +{ +public: + GravitySolverTest() + : state( std::make_unique< CommandLineOptions >( g_commandLineOptions )) + {} + +protected: + void SetUp() override + { + stations.clear(); + std::ostringstream coords; + coords << "{ "; + for( localIndex j = 0; j < NY; ++j ) + { + for( localIndex i = 0; i < NX; ++i ) + { + std::array< real64, 3 > station = { i * DX, j * DY, Z_STATION }; + stations.push_back( station ); + if( i != 0 || j != 0 ) coords << ", "; + coords << "{ " << station[0] << ", " << station[1] << ", " << station[2] << " }"; + } + } + coords << " }"; + + std::ostringstream xml; + xml << "\n" + << "\n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << " \n" + << "\n"; + + setupProblemFromXML( state.getProblemManager(), xml.str().c_str()); + } + + GeosxState state; + std::vector< std::array< real64, 3 > > stations; + GravityFE * gravitySolver; +}; + +TEST_F( GravitySolverTest, RectangularPrismAnomaly ) +{ + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + gravitySolver = &state.getProblemManager().getPhysicsSolverManager().getGroup< GravityFE >( "gravity" ); + + gravitySolver->explicitStepModeling( 0.0, 1.0e6, 0, domain ); + + auto const & gz_computed = gravitySolver->getReference< array1d< real64 > >( "gzAtStations" ); + + EXPECT_EQ( gz_computed.size(), stations.size()); + + for( localIndex i = 0; i < gz_computed.size(); ++i ) + { + real64 gz_expected = analyticPrismGz( stations[i], PRISM_BOUNDS, DENSITY_CONTRAST ); + EXPECT_NEAR( gz_computed[i], gz_expected, 1e-6 ) + << "Station " << i << " mismatch: computed=" << gz_computed[i] + << ", expected=" << gz_expected; + } +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; +} diff --git a/src/coreComponents/unitTests/gravityTests/testGravitySolverAdjoint.cpp b/src/coreComponents/unitTests/gravityTests/testGravitySolverAdjoint.cpp new file mode 100644 index 00000000000..dd7c4c4c7e6 --- /dev/null +++ b/src/coreComponents/unitTests/gravityTests/testGravitySolverAdjoint.cpp @@ -0,0 +1,202 @@ + +#include +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/gravity/GravityFE.hpp" +#include "mesh/DomainPartition.hpp" +#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" +#include "physicsSolvers/PhysicsSolverBase.hpp" +#include "finiteElement/FiniteElementDiscretization.hpp" + + +#include "common/DataTypes.hpp" +#include +#include +#include +#include + +using namespace geos; +using namespace geos::testing; +using namespace geos::dataRepository; + +using namespace geos::finiteElement; + +CommandLineOptions g_commandLineOptions; + +// Reusable adjoint test function +template< typename ForwardOp, typename AdjointOp > +std::pair< bool, double > adjointTest( + ForwardOp A, + AdjointOp AT, + std::vector< double > const & x, + std::vector< double > const & y, + double tol = 1e-13, + bool verbose = true ) +{ + auto dot = []( std::vector< double > const & a, std::vector< double > const & b ) -> double { + return std::inner_product( a.begin(), a.end(), b.begin(), 0.0 ); + }; + + auto norm = [&]( std::vector< double > const & v ) -> double { + return std::sqrt( dot( v, v )); + }; + + std::vector< double > Ax = A( x ); + std::vector< double > ATy = AT( y ); + + double IP1 = dot( Ax, y ); + double IP2 = dot( x, ATy ); + + double error = std::abs( IP1 - IP2 ) / std::max( norm( Ax ) * norm( y ), norm( ATy ) * norm( x )); + bool passed = error < tol; + + if( verbose ) + { + std::cout << "\n=== Adjoint Test Summary ===\n"; + std::cout << " = " << IP1 << "\n"; + std::cout << " = " << IP2 << "\n"; + std::cout << "Passed = " << std::boolalpha << passed << "\n"; + std::cout << "Error = " << error << "\n"; + } + + return { passed, error }; +} + +TEST( GravityFEKernelTest, AdjointConsistencyWithGravityFE ) +{ + GeosxState state( std::make_unique< CommandLineOptions >( g_commandLineOptions )); + + const char * xmlInput = + "" + "" + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + ""; + + setupProblemFromXML( state.getProblemManager(), xmlInput ); + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + GravityFE & gravityFE = state.getProblemManager().getPhysicsSolverManager().getGroup< GravityFE >( "gravity" ); + + // Set up synthetic station and residue + array2d< real64 > & stations = gravityFE.getReference< array2d< real64 > >( "stationCoordinates" ); + stations.resize( 3, 3 ); + stations( 0, 0 ) = 12557.71; stations( 0, 1 ) = 8075.03; stations( 0, 2 ) = 2527.50; + stations( 1, 0 ) = 10892.84; stations( 1, 1 ) = 10209.41; stations( 1, 2 ) = 2567.67; + stations( 2, 0 ) = 13568.72; stations( 2, 1 ) = 8260.82; stations( 2, 2 ) = 2542.19; + + + + array1d< real64 > & residueField = gravityFE.getReference< array1d< real64 > >( "residue" ); + + residueField.resize( 3 ); + residueField[0] = 1.0; + residueField[1] = 1.0; + residueField[2] = 1.0; + + + // Forward operator + auto A = [&]( std::vector< double > const & xVec ) -> std::vector< double > + { + gravityFE.forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( std::string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()); + for( localIndex i = 0; i < density.size(); ++i ) + { + const_cast< real64 & >(density[i]) = xVec[i]; + } + } ); + } ); + + + + gravityFE.explicitStepModeling( 0.0, 1.0, 0, domain ); + + arrayView1d< real64 const > gz = gravityFE.getReference< array1d< real64 > >( "gzAtStations" ).toViewConst(); + return std::vector< double >( gz.begin(), gz.end()); + }; + + + +// Adjoint operator + auto AT = [&]( std::vector< double > const & yVec ) -> std::vector< double > + { + // Set the residue field from yVec + array1d< real64 > & residueInner = gravityFE.getReference< array1d< real64 > >( "residue" ); + for( localIndex i = 0; i < residueInner.size(); ++i ) + residueInner[i] = yVec[i ]; + + // Run the adjoint step + gravityFE.explicitStepAdjoint( 0.0, 1.0, 0, domain ); + + // Extract adjoint field using the same loop structure + std::vector< double > adjointValues; + gravityFE.forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( std::string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + arrayView1d< real64 > const adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()); + for( localIndex i = 0; i < adjoint.size(); ++i ) + { + adjointValues.push_back( adjoint[i] ); + } + } ); + } ); + + return adjointValues; + }; + + + + std::vector< double > x( 12*12*12, 1.0 ); // synthetic density + + std::vector< double > y( 3, 2.0 ); // synthetic residue + + + + auto [passed, error] = adjointTest( A, AT, x, y ); + EXPECT_TRUE( passed ); + EXPECT_LT( error, 1e-13 ); + +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; +} From ea5766caa9872dfe0e0eb1733aef15fb26c44626 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Fri, 22 Aug 2025 16:25:01 -0500 Subject: [PATCH 08/30] Added first draft --- src/coreComponents/fileIO/CMakeLists.txt | 7 +- .../ScatterDataHistoryCollection.cpp | 235 ++++++++++++++++++ .../ScatterDataHistoryCollection.hpp | 154 ++++++++++++ .../timeHistory/ScatterDataProvider.hpp | 86 +++++++ 4 files changed, 480 insertions(+), 2 deletions(-) create mode 100644 src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp create mode 100644 src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp create mode 100644 src/coreComponents/fileIO/timeHistory/ScatterDataProvider.hpp diff --git a/src/coreComponents/fileIO/CMakeLists.txt b/src/coreComponents/fileIO/CMakeLists.txt index 2dfde481250..acc739f6df5 100644 --- a/src/coreComponents/fileIO/CMakeLists.txt +++ b/src/coreComponents/fileIO/CMakeLists.txt @@ -37,7 +37,9 @@ set( fileIO_headers timeHistory/BufferedHistoryIO.hpp timeHistory/PackCollection.hpp timeHistory/HDFHistoryIO.hpp - timeHistory/HistoryCollection.hpp ) + timeHistory/HistoryCollection.hpp + timeHistory/ScatterDataProvider.hpp + timeHistory/ScatterDataHistoryCollection.hpp ) # # Specify all sources @@ -54,7 +56,8 @@ set( fileIO_sources timeHistory/HDFFile.cpp timeHistory/HistoryCollectionBase.cpp timeHistory/PackCollection.cpp - timeHistory/HDFHistoryIO.cpp ) + timeHistory/HDFHistoryIO.cpp + timeHistory/ScatterDataHistoryCollection.cpp ) set( dependencyList ${parallelDeps} events mesh HDF5::HDF5 ) if( ENABLE_PYGEOSX ) diff --git a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp new file mode 100644 index 00000000000..660705268be --- /dev/null +++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp @@ -0,0 +1,235 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ScatterDataHistoryCollection.cpp + */ + +#include "ScatterDataHistoryCollection.hpp" + +#include "mesh/DomainPartition.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "common/MpiWrapper.hpp" +#include + +namespace geos +{ + +using namespace dataRepository; + +ScatterDataHistoryCollection::ScatterDataHistoryCollection( string const & name, Group * const parent ) + : HistoryCollectionBase( name, parent ), + m_scatterDataProvider( nullptr ) +{ + registerWrapper( viewKeyStruct::solverNameString(), &m_solverName ) + .setInputFlag( InputFlags::REQUIRED ) + .setDescription( "Name of the physics solver that provides scatter data" ); + + registerWrapper( viewKeyStruct::includeCoordinatesString(), &m_includeCoordinates ) + .setInputFlag( InputFlags::OPTIONAL ) + .setApplyDefaultValue( 1 ) + .setDescription( "Flag to include coordinates in output (1=yes, 0=no)" ); + + registerWrapper( viewKeyStruct::includeMetadataString(), &m_includeMetadata ) + .setInputFlag( InputFlags::OPTIONAL ) + .setApplyDefaultValue( 0 ) + .setDescription( "Flag to include metadata in output (1=yes, 0=no)" ); +} + +void ScatterDataHistoryCollection::initializePostInitialConditionsPreSubGroups() +{ + HistoryCollectionBase::initializePostInitialConditionsPreSubGroups(); + + // Find the domain partition to locate the physics solver + findScatterDataProvider(); + + // Set up the collection count + m_collectionCount = 1; // Main scatter data + + if( m_includeCoordinates && m_scatterDataProvider != nullptr ) + { + array2d< real64 > const & coordinates = m_scatterDataProvider->getScatterCoordinates(); + if( coordinates.size( 0 ) > 0 && coordinates.size( 1 ) > 0 ) + { + m_collectionCount += coordinates.size( 1 ); // Add X, Y, Z coordinates + } + } +} + +void ScatterDataHistoryCollection::findScatterDataProvider() +{ + // Get the problem manager from the hierarchy + Group const & problemManager = this->getGroupByPath( "/Problem" ); + + // Get the physics solver manager + Group const & physicsSolverManager = problemManager.getGroup( "Solvers" ); + + // Try to find the solver by name + Group const * solver = physicsSolverManager.getGroupPointer( m_solverName ); + GEOS_THROW_IF( solver == nullptr, + GEOS_FMT( "Could not find physics solver named '{}'", m_solverName ), + InputError ); + + // Cast to ScatterDataProvider + m_scatterDataProvider = dynamic_cast< ScatterDataProvider const * >( solver ); + GEOS_THROW_IF( m_scatterDataProvider == nullptr, + GEOS_FMT( "Physics solver '{}' does not implement ScatterDataProvider interface", m_solverName ), + InputError ); + + GEOS_LOG_RANK_0( GEOS_FMT( "Found scatter data provider '{}' with {} points", m_solverName, m_scatterDataProvider->getNumScatterPoints()) ); +} + +void ScatterDataHistoryCollection::collect( DomainPartition const & domain ) +{ + GEOS_UNUSED_VAR( domain ); + + if( m_scatterDataProvider == nullptr ) + { + GEOS_THROW( "Scatter data provider not initialized", std::runtime_error ); + } + + // Get the current scatter data values + array1d< real64 > const & scatterData = m_scatterDataProvider->getScatterData(); + localIndex const numPoints = m_scatterDataProvider->getNumScatterPoints(); + + GEOS_THROW_IF( scatterData.size() != numPoints, + GEOS_FMT( "Scatter data size ({}) does not match number of points ({})", scatterData.size(), numPoints ), + std::runtime_error ); +} + +void ScatterDataHistoryCollection::collect( DomainPartition const & domain, + localIndex const collectionIdx, + buffer_unit_type * & buffer ) +{ + GEOS_UNUSED_VAR( domain ); + GEOS_MARK_FUNCTION; + + if( m_scatterDataProvider == nullptr ) + { + return; + } + + localIndex const numPoints = m_scatterDataProvider->getNumScatterPoints(); + + // Collection index 0 is the main scatter data + if( collectionIdx == 0 ) + { + array1d< real64 > const & scatterData = m_scatterDataProvider->getScatterData(); + + // Copy data to buffer + std::memcpy( buffer, scatterData.data(), numPoints * sizeof( real64 ) ); + buffer += numPoints * sizeof( real64 ); + } + // Subsequent indices are coordinates if requested + else if( m_includeCoordinates ) + { + array2d< real64 > const & coordinates = m_scatterDataProvider->getScatterCoordinates(); + localIndex const numDim = coordinates.size( 1 ); + localIndex const coordIdx = collectionIdx - 1; + + if( coordIdx < numDim && coordinates.size( 0 ) == numPoints ) + { + // Extract the specific coordinate component + real64 * coordData = reinterpret_cast< real64 * >( buffer ); + for( localIndex i = 0; i < numPoints; ++i ) + { + coordData[i] = coordinates[i][coordIdx]; + } + buffer += numPoints * sizeof( real64 ); + } + } +} + +localIndex ScatterDataHistoryCollection::numCollectors() const +{ + localIndex count = 1; // Main scatter data + + if( m_includeCoordinates && m_scatterDataProvider != nullptr ) + { + array2d< real64 > const & coordinates = m_scatterDataProvider->getScatterCoordinates(); + if( coordinates.size( 0 ) > 0 && coordinates.size( 1 ) > 0 ) + { + count += coordinates.size( 1 ); // Add X, Y, Z coordinates + } + } + + return count; +} + +const string & ScatterDataHistoryCollection::getTargetName() const +{ + return m_solverName; +} + +HistoryMetadata ScatterDataHistoryCollection::getMetaData( DomainPartition const & domain, localIndex collectionIdx ) const +{ + GEOS_UNUSED_VAR( domain ); + + if( m_scatterDataProvider == nullptr ) + { + return HistoryMetadata(); + } + + localIndex const numPoints = m_scatterDataProvider->getNumScatterPoints(); + + // Collection index 0 is always the main scatter data + if( collectionIdx == 0 ) + { + return HistoryMetadata( "scatterData", numPoints, std::type_index( typeid( real64 ) ) ); + } + + // Subsequent indices are coordinates if requested + if( m_includeCoordinates ) + { + array2d< real64 > const & coordinates = m_scatterDataProvider->getScatterCoordinates(); + localIndex const numDim = coordinates.size( 1 ); + + if( collectionIdx <= numDim ) + { + string coordName; + switch( collectionIdx - 1 ) + { + case 0: coordName = "coordinateX"; break; + case 1: coordName = "coordinateY"; break; + case 2: coordName = "coordinateZ"; break; + default: coordName = "coordinate" + std::to_string( collectionIdx - 1 ); break; + } + return HistoryMetadata( coordName, numPoints, std::type_index( typeid( real64 ) ) ); + } + } + + return HistoryMetadata(); +} + +void ScatterDataHistoryCollection::updateSetsIndices( DomainPartition const & domain ) +{ + GEOS_UNUSED_VAR( domain ); +} + +localIndex ScatterDataHistoryCollection::numMetaDataCollectors() const +{ + return 0; +} + +HistoryCollection & ScatterDataHistoryCollection::getMetaDataCollector( localIndex metaIdx ) +{ + GEOS_UNUSED_VAR( metaIdx ); + GEOS_THROW( "Scatter data collections do not have metadata collectors", + std::runtime_error ); +} + +REGISTER_CATALOG_ENTRY( TaskBase, ScatterDataHistoryCollection, string const &, Group * const ) + +} // namespace geos diff --git a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp new file mode 100644 index 00000000000..66dc12ec1f9 --- /dev/null +++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp @@ -0,0 +1,154 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ScatterDataHistoryCollection.hpp + */ + +#ifndef GEOS_FILEIO_TIMEHISTORY_SCATTERDATAHISTORYCOLLECTION_HPP_ +#define GEOS_FILEIO_TIMEHISTORY_SCATTERDATAHISTORYCOLLECTION_HPP_ + +#include "HistoryCollectionBase.hpp" +#include "ScatterDataProvider.hpp" +#include "common/DataTypes.hpp" + +namespace geos +{ +class DomainPartition; + + +/** + * @brief Collects time history data from scattered points + * + * This class collects time history data from solvers that implement the + * ScatterDataProvider interface. + */ +class ScatterDataHistoryCollection : public HistoryCollectionBase +{ +public: + + /// Constructor + ScatterDataHistoryCollection( string const & name, Group * const parent ); + + /// Destructor + virtual ~ScatterDataHistoryCollection() = default; + + /// Non-copyable + ScatterDataHistoryCollection( ScatterDataHistoryCollection const & ) = delete; + ScatterDataHistoryCollection & operator=( ScatterDataHistoryCollection const & ) = delete; + + /// Movable: default move constructor, deleted move assignment operator + ScatterDataHistoryCollection( ScatterDataHistoryCollection && ) = default; + ScatterDataHistoryCollection & operator=( ScatterDataHistoryCollection && ) = delete; + + + /** + * @brief Get the catalog name for factory registration + * @return The catalog name + */ + static string catalogName() { return "ScatterDataHistoryCollection"; } + + /** + * @brief Initialize the collection + */ + virtual void initializePostInitialConditionsPreSubGroups() override; + + /** + * @brief Get the metadata for a specific collection index + * @param domain The domain partition + * @param collectionIdx The collection index + * @return HistoryMetadata describing the data to be collected + */ + virtual HistoryMetadata getMetaData( DomainPartition const & domain, localIndex collectionIdx ) const override; + + /** + * @brief Get the number of discrete collection operations + * @return Number of collection operations + */ + virtual localIndex numCollectors() const override; + + /** + * @brief Get the name of the target being collected + * @return The target name (solver name) + */ + virtual const string & getTargetName() const override; + + /** + * @brief Get the number of metadata collectors + * @return Number of metadata collectors (always 0 for scatter data) + */ + virtual localIndex numMetaDataCollectors() const override; + + /** + * @brief Get a metadata collector (not applicable for scatter data) + * @param metaIdx The metadata collector index + * @return Reference to metadata collector + */ + virtual HistoryCollection & getMetaDataCollector( localIndex metaIdx ) override; + +protected: + + /** + * @brief Update set indices (required by base class, no-op for scatter data) + * @param domain The domain partition + */ + virtual void updateSetsIndices( DomainPartition const & domain ) override; + + /** + * @brief Collect data from the scatter data provider + * @param domain The domain partition containing the data + */ + virtual void collect( DomainPartition const & domain ); + + /** + * @brief Collect data from a specific collection index into a buffer + * @param domain The domain partition containing the data + * @param collectionIdx The index of the collection operation + * @param buffer The buffer to write data into + */ + virtual void collect( DomainPartition const & domain, + localIndex const collectionIdx, + buffer_unit_type * & buffer ) override; + +private: + + /// Pointer to the scatter data provider (physics solver) + ScatterDataProvider const * m_scatterDataProvider; + + /// Name of the physics solver that provides scatter data + string m_solverName; + + /// Flag to include coordinates in output (default: true) + integer m_includeCoordinates; + + /// Flag to include metadata in output (default: false) + integer m_includeMetadata; + + /** + * @brief Find and validate the scatter data provider + */ + void findScatterDataProvider(); + + struct viewKeyStruct + { + static constexpr char const * solverNameString() { return "solverName"; } + static constexpr char const * includeCoordinatesString() { return "includeCoordinates"; } + static constexpr char const * includeMetadataString() { return "includeMetadata"; } + }; +}; + +} // namespace geos + +#endif /* GEOS_FILEIO_TIMEHISTORY_SCATTERDATAHISTORYCOLLECTION_HPP_ */ diff --git a/src/coreComponents/fileIO/timeHistory/ScatterDataProvider.hpp b/src/coreComponents/fileIO/timeHistory/ScatterDataProvider.hpp new file mode 100644 index 00000000000..025cfc49b47 --- /dev/null +++ b/src/coreComponents/fileIO/timeHistory/ScatterDataProvider.hpp @@ -0,0 +1,86 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ScatterDataProvider.hpp + */ + +#ifndef GEOS_FILEIO_TIMEHISTORY_SCATTERDATAPROVIDER_HPP_ +#define GEOS_FILEIO_TIMEHISTORY_SCATTERDATAPROVIDER_HPP_ + +#include "common/DataTypes.hpp" + +namespace geos +{ + +/** + * @brief Interface for classes that provide scattered data (coordinates and values) + */ +class ScatterDataProvider +{ +public: + + /// Virtual destructor + virtual ~ScatterDataProvider() = default; + + /** + * @brief Get the number of scatter points + * @return Number of points where scattered data is computed + */ + virtual localIndex getNumScatterPoints() const = 0; + + /** + * @brief Get the scattered data values + * @return Array of data values at scatter points + */ + virtual array1d< real64 > const & getScatterData() const = 0; + + /** + * @brief Get the coordinates of scatter points + * @return 2D array of coordinates [nPoints x nDim] + */ + virtual array2d< real64 > const & getScatterCoordinates() const + { + static array2d< real64 > empty; + return empty; + } + + /** + * @brief Get optional metadata for scatter points + * @return String array with metadata for each point (e.g., station names) + */ + virtual string_array const & getScatterMetadata() const + { + static string_array empty; + return empty; + } + +protected: + + /// Protected constructor + ScatterDataProvider() = default; + + /// Non-copyable + ScatterDataProvider( ScatterDataProvider const & ) = delete; + ScatterDataProvider & operator=( ScatterDataProvider const & ) = delete; + + /// Non-movable + ScatterDataProvider( ScatterDataProvider && ) = delete; + ScatterDataProvider & operator=( ScatterDataProvider && ) = delete; +}; + +} // namespace geos + +#endif /* GEOS_FILEIO_TIMEHISTORY_SCATTERDATAPROVIDER_HPP_ */ From 349a044b1bf4a5e9f4108e803d758a7d31e9cb82 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Fri, 22 Aug 2025 22:22:04 -0500 Subject: [PATCH 09/30] Added test. --- .../ScatterDataHistoryCollection.cpp | 19 +- .../ScatterDataHistoryCollection.hpp | 10 +- .../unitTests/fileIOTests/CMakeLists.txt | 3 +- .../testScatterDataHistoryCollection.cpp | 220 ++++++++++++++++++ 4 files changed, 236 insertions(+), 16 deletions(-) create mode 100644 src/coreComponents/unitTests/fileIOTests/testScatterDataHistoryCollection.cpp diff --git a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp index 660705268be..85ea2760059 100644 --- a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp +++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp @@ -100,13 +100,15 @@ void ScatterDataHistoryCollection::collect( DomainPartition const & domain ) GEOS_THROW( "Scatter data provider not initialized", std::runtime_error ); } - // Get the current scatter data values - array1d< real64 > const & scatterData = m_scatterDataProvider->getScatterData(); - localIndex const numPoints = m_scatterDataProvider->getNumScatterPoints(); - - GEOS_THROW_IF( scatterData.size() != numPoints, - GEOS_FMT( "Scatter data size ({}) does not match number of points ({})", scatterData.size(), numPoints ), - std::runtime_error ); + // For each collection, populate the output buffer using the provider + for( localIndex collectionIdx = 0; collectionIdx < m_collectionCount; ++collectionIdx ) + { + // Get the metadata to determine the buffer size needed + HistoryMetadata hmd = this->getMetaData( domain, collectionIdx ); + // Get the output buffer for this collection using the buffer provider callback + buffer_unit_type * buffer = m_bufferProviders[collectionIdx]( hmd.size( 0 )); + collect( domain, collectionIdx, buffer ); + } } void ScatterDataHistoryCollection::collect( DomainPartition const & domain, @@ -127,10 +129,8 @@ void ScatterDataHistoryCollection::collect( DomainPartition const & domain, if( collectionIdx == 0 ) { array1d< real64 > const & scatterData = m_scatterDataProvider->getScatterData(); - // Copy data to buffer std::memcpy( buffer, scatterData.data(), numPoints * sizeof( real64 ) ); - buffer += numPoints * sizeof( real64 ); } // Subsequent indices are coordinates if requested else if( m_includeCoordinates ) @@ -147,7 +147,6 @@ void ScatterDataHistoryCollection::collect( DomainPartition const & domain, { coordData[i] = coordinates[i][coordIdx]; } - buffer += numPoints * sizeof( real64 ); } } } diff --git a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp index 66dc12ec1f9..de5bf23ac08 100644 --- a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp +++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp @@ -53,6 +53,11 @@ class ScatterDataHistoryCollection : public HistoryCollectionBase ScatterDataHistoryCollection( ScatterDataHistoryCollection && ) = default; ScatterDataHistoryCollection & operator=( ScatterDataHistoryCollection && ) = delete; + /** + * @brief Collect data from the scatter data provider + * @param domain The domain partition containing the data + */ + virtual void collect( DomainPartition const & domain ); /** * @brief Get the catalog name for factory registration @@ -106,11 +111,6 @@ class ScatterDataHistoryCollection : public HistoryCollectionBase */ virtual void updateSetsIndices( DomainPartition const & domain ) override; - /** - * @brief Collect data from the scatter data provider - * @param domain The domain partition containing the data - */ - virtual void collect( DomainPartition const & domain ); /** * @brief Collect data from a specific collection index into a buffer diff --git a/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt b/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt index 6dde82ba541..27c91c32e16 100644 --- a/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt @@ -1,7 +1,8 @@ # Specify list of tests set( geosx_fileio_tests testHDFFile.cpp - testMemoryStats.cpp ) + testMemoryStats.cpp + testScatterDataHistoryCollection.cpp ) set( tplDependencyList ${parallelDeps} gtest ) diff --git a/src/coreComponents/unitTests/fileIOTests/testScatterDataHistoryCollection.cpp b/src/coreComponents/unitTests/fileIOTests/testScatterDataHistoryCollection.cpp new file mode 100644 index 00000000000..e6d63d1d316 --- /dev/null +++ b/src/coreComponents/unitTests/fileIOTests/testScatterDataHistoryCollection.cpp @@ -0,0 +1,220 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "fileIO/timeHistory/ScatterDataHistoryCollection.hpp" +#include "fileIO/Outputs/TimeHistoryOutput.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "mesh/DomainPartition.hpp" +#include "mainInterface/initialization.hpp" +#include +#include +#include +#include + +using namespace geos; + +// Minimal mock ScatterDataProvider +class MockScatterDataProvider : public ScatterDataProvider, public dataRepository::Group +{ +public: + MockScatterDataProvider( string const & name, dataRepository::Group *parent ): dataRepository::Group( name, parent ) + { + m_scatterData.resize( 3 ); + m_scatterData[0] = 10; m_scatterData[1] = 20; m_scatterData[2] = 30; + m_coordinates.resize( 3, 3 ); + m_coordinates( 0, 0 ) = 1000.0; m_coordinates( 0, 1 ) = 2000.0; m_coordinates( 0, 2 ) = 100.0; + m_coordinates( 1, 0 ) = 1500.0; m_coordinates( 1, 1 ) = 2500.0; m_coordinates( 1, 2 ) = 200.0; + m_coordinates( 2, 0 ) = 2000.0; m_coordinates( 2, 1 ) = 3000.0; m_coordinates( 2, 2 ) = 300.0; + m_metadata.resize( 3 ); + m_metadata[0] = "Station_A"; m_metadata[1] = "Station_B"; m_metadata[2] = "Station_C"; + } + localIndex getNumScatterPoints() const override { return 3; } + array1d< real64 > const & getScatterData() const override { return m_scatterData; } + array2d< real64 > const & getScatterCoordinates() const override { return m_coordinates; } + string_array const & getScatterMetadata() const override { return m_metadata; } + static string catalogName() { return "MockScatterDataProvider"; } +private: + array1d< real64 > m_scatterData; + array2d< real64 > m_coordinates; + string_array m_metadata; +}; + +REGISTER_CATALOG_ENTRY( dataRepository::Group, MockScatterDataProvider, string const &, dataRepository::Group * const ) + + +class ScatterDataHistoryIntegrationTest : public ::testing::Test +{ +protected: + void SetUp() override + { + m_conduitNode = std::make_unique< conduit::Node >(); + m_rootGroup = std::make_unique< dataRepository::Group >( "Problem", *m_conduitNode ); + auto & solvers = m_rootGroup->registerGroup< PhysicsSolverManager >( "Solvers" ); + m_domain = &m_rootGroup->registerGroup< DomainPartition >( "domain" ); + solvers.registerGroup< MockScatterDataProvider >( "gravitySolver" ); + m_collection = &m_rootGroup->registerGroup< ScatterDataHistoryCollection >( "gravityHistory" ); + m_collection->getReference< string >( "solverName" ) = "gravitySolver"; + m_collection->getReference< integer >( "includeCoordinates" ) = 1; + m_collection->getReference< integer >( "includeMetadata" ) = 0; + } + std::unique_ptr< conduit::Node > m_conduitNode; + std::unique_ptr< dataRepository::Group > m_rootGroup; + DomainPartition *m_domain; + ScatterDataHistoryCollection *m_collection; +}; + +TEST_F( ScatterDataHistoryIntegrationTest, CompleteWorkflow ) +{ + // Test initialization + m_collection->initializePostInitialConditionsPreSubGroups(); + EXPECT_EQ( m_collection->numCollectors(), 4 ); // data + 3 coordinates + EXPECT_EQ( m_collection->getTargetName(), "gravitySolver" ); + EXPECT_EQ( m_collection->catalogName(), "ScatterDataHistoryCollection" ); + EXPECT_EQ( m_collection->numMetaDataCollectors(), 0 ); + EXPECT_THROW( m_collection->getMetaDataCollector( 0 ), std::runtime_error ); + + // Test metadata for all collectors + for( localIndex idx = 0; idx < 4; ++idx ) + { + auto metadata = m_collection->getMetaData( *m_domain, idx ); + EXPECT_EQ( metadata.size(), 3 ); + EXPECT_EQ( metadata.getType(), std::type_index( typeid(real64))); + if( idx == 0 ) + EXPECT_EQ( metadata.getName(), "scatterData" ); + else + EXPECT_TRUE( metadata.getName().find( "coordinate" ) != string::npos ); + } +} + +TEST_F( ScatterDataHistoryIntegrationTest, ErrorHandling ) +{ + m_collection->getReference< string >( "solverName" ) = "invalidSolver"; + EXPECT_THROW( m_collection->initializePostInitialConditionsPreSubGroups(), InputError ); +} + +// HDF5 test +class ScatterDataHistoryHDF5Test : public ::testing::Test +{ +protected: + void SetUp() override + { + m_conduitNode = std::make_unique< conduit::Node >(); + m_rootGroup = std::make_unique< dataRepository::Group >( "Problem", *m_conduitNode ); + auto & solvers = m_rootGroup->registerGroup< PhysicsSolverManager >( "Solvers" ); + m_provider = &solvers.registerGroup< MockScatterDataProvider >( "gravitySolver" ); + m_domain = &m_rootGroup->registerGroup< DomainPartition >( "domain" ); + auto & outputs = m_rootGroup->registerGroup< dataRepository::Group >( "Outputs" ); + m_timeHistoryOutput = &outputs.registerGroup< TimeHistoryOutput >( "timeHistory" ); + m_timeHistoryOutput->getReference< string >( "filename" ) = "testOutput"; // No path, just filename + m_timeHistoryOutput->getReference< string >( "format" ) = "hdf5"; + m_collection = &m_timeHistoryOutput->registerGroup< ScatterDataHistoryCollection >( "gravityHistory" ); + m_collection->getReference< string >( "solverName" ) = "gravitySolver"; + m_collection->getReference< integer >( "includeCoordinates" ) = 1; + // Register the collection as a source for TimeHistoryOutput using string_array + m_timeHistoryOutput->getReference< string_array >( "sources" ).resize( 1 ); + m_timeHistoryOutput->getReference< string_array >( "sources" )[0] = "gravityHistory"; + } + void TearDown() override + { + std::filesystem::remove( "testOutput.hdf5" ); + } + std::unique_ptr< conduit::Node > m_conduitNode; + std::unique_ptr< dataRepository::Group > m_rootGroup; + DomainPartition *m_domain; + TimeHistoryOutput *m_timeHistoryOutput; + ScatterDataHistoryCollection *m_collection; + MockScatterDataProvider *m_provider; +}; + +TEST_F( ScatterDataHistoryHDF5Test, HDF5FileOutput ) +{ + // Set up output directory for TimeHistoryOutput + OutputBase::setOutputDirectory( "." ); + + // Initialize collection first + m_collection->initializePostInitialConditionsPreSubGroups(); + + // Debug: Check the sources array is properly set + auto & sources = m_timeHistoryOutput->getReference< string_array >( "sources" ); + EXPECT_EQ( sources.size(), 1 ); + EXPECT_EQ( sources[0], "gravityHistory" ); + + // Initialize TimeHistoryOutput - this should now work without directory errors + EXPECT_NO_THROW( m_timeHistoryOutput->initializePostInitialConditionsPostSubGroups() ); + + // Verify the collection is configured correctly + EXPECT_EQ( m_collection->numCollectors(), 4 ); // data + 3 coordinates + EXPECT_EQ( m_collection->getTargetName(), "gravitySolver" ); + + // Verify that we can get metadata for all collectors + for( localIndex idx = 0; idx < 4; ++idx ) + { + auto metadata = m_collection->getMetaData( *m_domain, idx ); + EXPECT_EQ( metadata.size(), 3 ); // 3 scatter points + EXPECT_EQ( metadata.getType(), std::type_index( typeid(real64) ) ); + if( idx == 0 ) + EXPECT_EQ( metadata.getName(), "scatterData" ); + else + EXPECT_TRUE( metadata.getName().find( "coordinate" ) != string::npos ); + } + + // Trigger buffer population before writing + m_collection->collect( *m_domain ); + + // Now test the actual I/O execution! + EXPECT_NO_THROW( m_timeHistoryOutput->execute( 0.0, 1.0, 0, 1, 1.0, *m_domain ) ); + + // Verify HDF5 file was created + EXPECT_TRUE( std::filesystem::exists( "testOutput.hdf5" ) ); + + // Read back the HDF5 file and verify the contents + hid_t file_id = H5Fopen( "testOutput.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT ); + EXPECT_GE( file_id, 0 ); + + // Check scatter data values + { + const char * dataset_names[] = { "/scatterData", "/coordinateX", "/coordinateY", "/coordinateZ" }; + double values[3]; + auto & scatterData = m_provider->getScatterData(); + auto & coordinates = m_provider->getScatterCoordinates(); + for( int i = 0; i < 4; ++i ) + { + hid_t dataset_id = H5Dopen2( file_id, dataset_names[i], H5P_DEFAULT ); + EXPECT_GE( dataset_id, 0 ) << "Failed to open dataset " << dataset_names[i]; + herr_t status = H5Dread( dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, values ); + EXPECT_GE( status, 0 ) << "Failed to read dataset " << dataset_names[i]; + for( int j = 0; j < 3; ++j ) + { + double expected = 0.0; + if( i == 0 ) + expected = scatterData[j]; + else + expected = coordinates( j, i-1 ); + EXPECT_DOUBLE_EQ( values[j], expected ) << "Mismatch in " << dataset_names[i] << " at index " << j; + } + H5Dclose( dataset_id ); + } + } + H5Fclose( file_id ); +} + +int main( int argc, char * *argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + basicSetup( argc, argv ); + int result = RUN_ALL_TESTS(); + basicCleanup(); + return result; +} From ad39d3b33325dc86ecc04c650e52c1c11040b813 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Fri, 22 Aug 2025 23:03:55 -0500 Subject: [PATCH 10/30] Added scatter data provider --- .../physicsSolvers/gravity/GravityFE.cpp | 31 +++++++++++++++++++ .../physicsSolvers/gravity/GravityFE.hpp | 12 ++++--- .../gravity/GravitySolverBase.hpp | 3 ++ 3 files changed, 42 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 881bf006849..c293ec438e8 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -23,6 +23,7 @@ #include "finiteElement/FiniteElementDiscretization.hpp" #include "mesh/DomainPartition.hpp" +#include "common/MpiWrapper.hpp" namespace geos @@ -326,6 +327,36 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, return dt; } +localIndex GravityFE::getNumScatterPoints() const +{ + // For gravity applications: Only rank 0 should provide scatter data + // since all ranks compute the same gravity values at the same stations + if( MpiWrapper::commRank( MPI_COMM_GEOS ) != 0 ) + { + return 0; + } + + return m_stationCoordinates.size(0); // Return number of rows (stations) +} + +array1d< real64 > const & GravityFE::getScatterData() const +{ + GEOS_ASSERT( m_gzAtStations.size() == m_stationCoordinates.size(0) ); + return m_gzAtStations; +} + +array2d< real64 > const & GravityFE::getScatterCoordinates() const +{ + return m_stationCoordinates; +} + +string_array const & GravityFE::getScatterMetadata() const +{ + // Just return empty metadata for now + static string_array empty; + return empty; +} + REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE, string const &, dataRepository::Group * const ) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp index 85d9bdf1833..0141feaed1f 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp @@ -23,15 +23,14 @@ #include "GravitySolverBase.hpp" #include "mesh/MeshFields.hpp" - +#include "fileIO/timeHistory/ScatterDataProvider.hpp" namespace geos { constexpr real64 GRAVITATIONAL_CONSTANT = 6.67430e-11; // in m³·kg⁻¹·s⁻² - -class GravityFE : public GravitySolverBase +class GravityFE : public GravitySolverBase, public ScatterDataProvider { public: @@ -81,6 +80,11 @@ class GravityFE : public GravitySolverBase DomainPartition & domain ) override; /**@}*/ + // ScatterDataProvider interface + virtual localIndex getNumScatterPoints() const override; + virtual array1d< real64 > const & getScatterData() const override; + virtual array2d< real64 > const & getScatterCoordinates() const override; + virtual string_array const & getScatterMetadata() const override; protected: @@ -129,4 +133,4 @@ DECLARE_FIELD( Adjoint, } /* namespace geos */ -#endif /* GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ */ +#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp index 478e94e95ab..566add3d828 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp @@ -119,6 +119,9 @@ class GravitySolverBase : public PhysicsSolverBase /// Gz component recorded at stations array1d< real64 > m_gzAtStations; + + /// Station metadata (names/identifiers) + mutable string_array m_stationMetadata; }; } // namespace geos From 1bd4997026016b02c7029269160be7a3233c9ae1 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Sat, 23 Aug 2025 11:11:17 -0500 Subject: [PATCH 11/30] Fixed clang. --- .../fileIO/timeHistory/ScatterDataHistoryCollection.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp index de5bf23ac08..251db4cad10 100644 --- a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp +++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp @@ -49,9 +49,6 @@ class ScatterDataHistoryCollection : public HistoryCollectionBase ScatterDataHistoryCollection( ScatterDataHistoryCollection const & ) = delete; ScatterDataHistoryCollection & operator=( ScatterDataHistoryCollection const & ) = delete; - /// Movable: default move constructor, deleted move assignment operator - ScatterDataHistoryCollection( ScatterDataHistoryCollection && ) = default; - ScatterDataHistoryCollection & operator=( ScatterDataHistoryCollection && ) = delete; /** * @brief Collect data from the scatter data provider From 607bdabc83a074d40f3156a0e8c38a2929a6ce93 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Sat, 23 Aug 2025 11:13:39 -0500 Subject: [PATCH 12/30] Fixed doc. --- .../fileIO/timeHistory/ScatterDataHistoryCollection.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp index 251db4cad10..3a27b48ebec 100644 --- a/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp +++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp @@ -39,7 +39,11 @@ class ScatterDataHistoryCollection : public HistoryCollectionBase { public: - /// Constructor + /** + * @brief Constructor for ScatterDataHistoryCollection + * @param name The name of the collection + * @param parent The parent group in the data repository + */ ScatterDataHistoryCollection( string const & name, Group * const parent ); /// Destructor From 12dc58fe1a1a61f030f8e4dfd959af67cd685b6e Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Sat, 23 Aug 2025 12:51:26 -0500 Subject: [PATCH 13/30] Added test for scatter data provider. --- .../unitTests/fileIOTests/CMakeLists.txt | 1 + .../fileIOTests/testScatterDataProvider.cpp | 90 +++++++++++++++++++ 2 files changed, 91 insertions(+) create mode 100644 src/coreComponents/unitTests/fileIOTests/testScatterDataProvider.cpp diff --git a/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt b/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt index 27c91c32e16..f5f7bd267f7 100644 --- a/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/fileIOTests/CMakeLists.txt @@ -2,6 +2,7 @@ set( geosx_fileio_tests testHDFFile.cpp testMemoryStats.cpp + testScatterDataProvider.cpp testScatterDataHistoryCollection.cpp ) set( tplDependencyList ${parallelDeps} gtest ) diff --git a/src/coreComponents/unitTests/fileIOTests/testScatterDataProvider.cpp b/src/coreComponents/unitTests/fileIOTests/testScatterDataProvider.cpp new file mode 100644 index 00000000000..ce3536d1ae1 --- /dev/null +++ b/src/coreComponents/unitTests/fileIOTests/testScatterDataProvider.cpp @@ -0,0 +1,90 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "fileIO/timeHistory/ScatterDataProvider.hpp" +#include + +using namespace geos; + +class MinimalScatterDataProvider : public ScatterDataProvider +{ +public: + MinimalScatterDataProvider() + { + m_scatterData.resize( 0 ); + m_coordinates.resize( 0, 3 ); + m_metadata.resize( 0 ); + } + localIndex getNumScatterPoints() const override { return 0; } + array1d< real64 > const & getScatterData() const override { return m_scatterData; } + array2d< real64 > const & getScatterCoordinates() const override { return m_coordinates; } + string_array const & getScatterMetadata() const override { return m_metadata; } +private: + array1d< real64 > m_scatterData; + array2d< real64 > m_coordinates; + string_array m_metadata; +}; + +TEST( ScatterDataProviderTest, ZeroScatterPoints ) +{ + MinimalScatterDataProvider provider; + EXPECT_EQ( provider.getNumScatterPoints(), 0 ); + EXPECT_EQ( provider.getScatterData().size(), 0 ); + EXPECT_EQ( provider.getScatterCoordinates().size( 0 ), 0 ); + EXPECT_EQ( provider.getScatterCoordinates().size( 1 ), 3 ); + EXPECT_EQ( provider.getScatterMetadata().size(), 0 ); +} + +class EdgeCaseScatterDataProvider : public ScatterDataProvider +{ +public: + EdgeCaseScatterDataProvider() + { + m_scatterData.resize( 2 ); + m_scatterData[0] = -1e10; m_scatterData[1] = 1e10; + m_coordinates.resize( 2, 3 ); + m_coordinates( 0, 0 ) = -1e5; m_coordinates( 0, 1 ) = 0.0; m_coordinates( 0, 2 ) = 1e5; + m_coordinates( 1, 0 ) = 1e-5; m_coordinates( 1, 1 ) = -1e-5; m_coordinates( 1, 2 ) = 0.0; + m_metadata.resize( 2 ); + m_metadata[0] = ""; + m_metadata[1] = "Special_Station_#1"; + } + localIndex getNumScatterPoints() const override { return 2; } + array1d< real64 > const & getScatterData() const override { return m_scatterData; } + array2d< real64 > const & getScatterCoordinates() const override { return m_coordinates; } + string_array const & getScatterMetadata() const override { return m_metadata; } +private: + array1d< real64 > m_scatterData; + array2d< real64 > m_coordinates; + string_array m_metadata; +}; + +TEST( ScatterDataProviderTest, EdgeCases ) +{ + EdgeCaseScatterDataProvider provider; + EXPECT_EQ( provider.getNumScatterPoints(), 2 ); + EXPECT_DOUBLE_EQ( provider.getScatterData()[0], -1e10 ); + EXPECT_DOUBLE_EQ( provider.getScatterData()[1], 1e10 ); + EXPECT_DOUBLE_EQ( provider.getScatterCoordinates()( 0, 0 ), -1e5 ); + EXPECT_DOUBLE_EQ( provider.getScatterCoordinates()( 0, 2 ), 1e5 ); + EXPECT_EQ( provider.getScatterMetadata()[0], "" ); + EXPECT_EQ( provider.getScatterMetadata()[1], "Special_Station_#1" ); +} + +int main( int argc, char * *argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + return RUN_ALL_TESTS(); +} From 429d43612479abbe9fb591895deda57014b20956 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Sun, 24 Aug 2025 21:03:12 -0500 Subject: [PATCH 14/30] Removed SEP output format. --- .../physicsSolvers/gravity/GravityFE.cpp | 17 ++----- .../gravity/GravitySolverBase.cpp | 47 ------------------- .../gravity/GravitySolverBase.hpp | 9 ---- .../gravityTests/testGravitySolver.cpp | 5 +- 4 files changed, 8 insertions(+), 70 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index c293ec438e8..204d59246ef 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -105,8 +105,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, DomainPartition & domain ) { GEOS_MARK_FUNCTION; - - int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); + GEOS_UNUSED_VAR( time_n, cycleNumber ); array1d< real64 > localGzAtStations( m_stationCoordinates.size( 0 ) ); localGzAtStations.setValues< parallelHostPolicy >( 0. ); @@ -224,12 +223,6 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, } } - // Dump result to disk... - if( (rank==0) && ( this->m_outputGz == 1 )) - { - GravitySolverBase::saveGz( time_n, cycleNumber, this->m_outputGzBasename, gzAtStations ); - } - return dt; } @@ -335,14 +328,14 @@ localIndex GravityFE::getNumScatterPoints() const { return 0; } - - return m_stationCoordinates.size(0); // Return number of rows (stations) + + return m_stationCoordinates.size( 0 ); // Return number of rows (stations) } array1d< real64 > const & GravityFE::getScatterData() const { - GEOS_ASSERT( m_gzAtStations.size() == m_stationCoordinates.size(0) ); - return m_gzAtStations; + GEOS_ASSERT( m_gzAtStations.size() == m_stationCoordinates.size( 0 ) ); + return m_gzAtStations; } array2d< real64 > const & GravityFE::getScatterCoordinates() const diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp index 5989ef81c31..3e5027657a2 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp @@ -39,16 +39,6 @@ GravitySolverBase::GravitySolverBase( const std::string & name, Group * const pa .setSizedFromParent( 0 ) .setDescription( "Coordinates (x,y,z) of the gravimeter stations" ); - registerWrapper( viewKeyStruct::outputGzString(), &m_outputGz ) - .setInputFlag( InputFlags::OPTIONAL ) - .setApplyDefaultValue( 1 ) - .setDescription( "Flag to dump to disk field recorded at gravimeters" ); - - registerWrapper( viewKeyStruct::outputGzBasenameString(), &m_outputGzBasename ) - .setInputFlag( InputFlags::OPTIONAL ) - .setApplyDefaultValue( "gz" ) - .setDescription( "Basename to dump to disk field recorded at gravimeters" ); - registerWrapper( viewKeyStruct::residueString(), &m_residue ) .setInputFlag( InputFlags::OPTIONAL ) .setSizedFromParent( 0 ) @@ -135,8 +125,6 @@ void GravitySolverBase::initializePostInitialConditionsPreSubGroups() GEOS_LOG_RANK_0( "Name: " << getName()); GEOS_LOG_RANK_0( "Mode: " << m_modeString ); GEOS_LOG_RANK_0( "Number of stations: " << m_stationCoordinates.size( 0 )); - GEOS_LOG_RANK_0( "Output Gz to file: " << yesno( m_outputGz )); - GEOS_LOG_RANK_0( "Output Gz basename: " << m_outputGzBasename ); GEOS_LOG_RANK_0( "Log level: " << getLogLevel()); GEOS_LOG_RANK_0( " Output Gz to logs: " << yesno( getLogLevel() > 1 )); GEOS_LOG_RANK_0( " Output Adjoint to logs: " << yesno( getLogLevel() > 2 )); @@ -181,39 +169,4 @@ real64 GravitySolverBase::explicitStep( real64 const & time_n, } } -void GravitySolverBase::saveGz( real64 const & time_n, - integer const cycleNumber, - string const basename, - arrayView1d< real64 > const & gzAtStations ) -{ - // Convert to float32 - std::vector< float > tmp32( gzAtStations.size()); - std::transform( gzAtStations.begin(), gzAtStations.end(), tmp32.begin(), - [] ( real64 val )-> float { return static_cast< float >(val); } ); - - // Binary data - std::filesystem::path dataFilename = basename + GEOS_FMT( "_{:015}.H@", time_n ); - std::ofstream fout( dataFilename, std::ios::binary ); - if( !fout ) - throw std::ios_base::failure( "Failed to open data file" ); - fout.write( reinterpret_cast< const char * >(tmp32.data()), tmp32.size() * sizeof(float)); - - // Header file - std::filesystem::path headerFilename = basename + GEOS_FMT( "_{:015}.H", time_n ); - std::ofstream fheader( headerFilename ); - if( !fheader ) - throw std::ios_base::failure( "Failed to open header file" ); - - fheader << "n1=" << gzAtStations.size() << '\n' - << "d1=1.\n" - << "o1=0.\n" - << "label1=STATION\n" - << "esize=4\n" - << "data_format=native_float\n" - << "data_label=GZ\n" - << "in=" << dataFilename << '\n' - << GEOS_FMT( "time={:015}", time_n ) << '\n' - << "cycle=" << cycleNumber << '\n'; -} - } // namespace geos diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp index 566add3d828..3019be12bc0 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp @@ -102,18 +102,9 @@ class GravitySolverBase : public PhysicsSolverBase integer const cycleNumber, DomainPartition & domain ) = 0; - void saveGz( real64 const & time_n, - integer const cycleNumber, - string const basename, - arrayView1d< real64 > const & gzAtStations ); - /// Coordinates of the gravimeter stations array2d< real64 > m_stationCoordinates; - /// Dump vertical component Gz to disk - localIndex m_outputGz; - string m_outputGzBasename; - /// Residue observed at station (only for adjoint computation) array1d< real64 > m_residue; diff --git a/src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp b/src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp index 01b00f051f0..0cde3413cfe 100644 --- a/src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp +++ b/src/coreComponents/unitTests/gravityTests/testGravitySolver.cpp @@ -52,6 +52,8 @@ double analyticPrismGz( const std::array< double, 3 > & station, double density_contrast ) { const double x = station[0]; + + const double y = station[1]; const double z = station[2]; @@ -129,8 +131,7 @@ class GravitySolverTest : public ::testing::Test << "\n" << " \n" << " \n" + << " mode=\"modeling\" stationCoordinates=\"" << coords.str() <<"\"/>\n" << " \n" << " \n" << " Date: Sun, 24 Aug 2025 21:05:09 -0500 Subject: [PATCH 15/30] Added modeling input files. --- .../analytic_prism_comparison.py | 209 +++ .../basicRectPrism_C3D8.xml | 1369 +++++++++++++++++ .../extract_gravity_data.py | 107 ++ 3 files changed, 1685 insertions(+) create mode 100644 inputFiles/gravityModeling/basicRectPrismC3D8/analytic_prism_comparison.py create mode 100644 inputFiles/gravityModeling/basicRectPrismC3D8/basicRectPrism_C3D8.xml create mode 100644 inputFiles/gravityModeling/basicRectPrismC3D8/extract_gravity_data.py diff --git a/inputFiles/gravityModeling/basicRectPrismC3D8/analytic_prism_comparison.py b/inputFiles/gravityModeling/basicRectPrismC3D8/analytic_prism_comparison.py new file mode 100644 index 00000000000..8cf0bfe3567 --- /dev/null +++ b/inputFiles/gravityModeling/basicRectPrismC3D8/analytic_prism_comparison.py @@ -0,0 +1,209 @@ +#!/usr/bin/env python3 +""" +Compare computed gravity data from HDF5 file to analytic rectangular prism solution. + +Analytic solution reference: +Nagy, D. (1966). "The gravitational attraction of a right rectangular prism." Geophysics, 31(2), 362-371. https://doi.org/10.1190/1.1439779 +""" + +import h5py +import numpy as np +import sys + +# Analytic solution for gravity anomaly (gz) due to a rectangular prism +def analyticPrismGz(station, prism_bounds, density_contrast): + x, y, z = station + x0, x1, y0, y1, z0, z1 = prism_bounds + epsilon = 1e-16 + val = 0.0 + for i in range(2): + for j in range(2): + for k in range(2): + xi = x0 if i == 0 else x1 + yj = y0 if j == 0 else y1 + zk = z0 if k == 0 else z1 + dx = xi - x + dy = yj - y + dz = zk - z + r = np.sqrt(dx*dx + dy*dy + dz*dz) + if r > epsilon: + sign = 1.0 if (i + j + k) % 2 == 0 else -1.0 + arg1 = np.arctan2(dx * dy, dz * r + epsilon) + arg2 = np.log(r + dy + epsilon) + arg3 = np.log(r + dx + epsilon) + val += sign * (dz * arg1 - dx * arg2 - dy * arg3) + G = 6.67430e-11 # m^3 kg^-1 s^-2 + return G * density_contrast * val + + +def write_vtk_surface(filename, x_unique, y_unique, z_val, gravity_2d, analytic_2d=None): + """Write gravity data as a VTK RECTILINEAR_GRID surface.""" + nx, ny = len(x_unique), len(y_unique) + with open(filename, 'w') as f: + f.write('# vtk DataFile Version 3.0\n') + f.write('Gravity surface data\n') + f.write('ASCII\n') + f.write(f'DATASET RECTILINEAR_GRID\n') + f.write(f'DIMENSIONS {nx} {ny} 1\n') + f.write(f'X_COORDINATES {nx} float\n') + for xi in x_unique: + f.write(f'{xi} ') + f.write('\n') + f.write(f'Y_COORDINATES {ny} float\n') + for yi in y_unique: + f.write(f'{yi} ') + f.write('\n') + f.write(f'Z_COORDINATES 1 float\n') + f.write(f'{z_val}\n') + f.write(f'\nPOINT_DATA {nx*ny}\n') + f.write('SCALARS gravity float 1\n') + f.write('LOOKUP_TABLE default\n') + for iy in range(ny): + for ix in range(nx): + f.write(f'{gravity_2d[ix, iy]}\n') + if analytic_2d is not None: + f.write('\nSCALARS analytic float 1\n') + f.write('LOOKUP_TABLE default\n') + for iy in range(ny): + for ix in range(nx): + f.write(f'{analytic_2d[ix, iy]}\n') + + +def compare_gravity_data(hdf5_file, prism_bounds, density_contrast, save_fig=None): + """Extract gravity data and coordinates from HDF5 file and compare to analytic solution""" + try: + with h5py.File(hdf5_file, 'r') as f: + # Read coordinate data + coord_x = f['/coordinateX'][:] + coord_y = f['/coordinateY'][:] + coord_z = f['/coordinateZ'][:] + # Read gravity data + gravity_data = f['/scatterData'][:] + # Read time data + time_data = f['/gravity Time'][:] + # Assume single time step for simplicity + if len(coord_x.shape) > 1: + current_x = coord_x[0, :] + current_y = coord_y[0, :] + current_z = coord_z[0, :] + current_gravity = gravity_data[0, :] + else: + current_x = coord_x + current_y = coord_y + current_z = coord_z + current_gravity = gravity_data + stations = np.stack([current_x, current_y, -current_z], axis=-1) + # Compute analytic solution + gz_analytic = np.array([analyticPrismGz(station, prism_bounds, density_contrast) for station in stations]) + # Compare + diff = current_gravity - gz_analytic + for i, (g_comp, g_ana, d) in enumerate(zip(current_gravity, gz_analytic, diff)): + print(f"Station {i:4d}: Computed={g_comp:15.6e} Analytic={g_ana:15.6e} Diff={d:15.6e}") + print("Comparison to analytic solution:") + print(f" Max abs difference: {np.max(np.abs(diff)):.6e}") + print(f" Mean abs difference: {np.mean(np.abs(diff)):.6e}") + + # Plot results as 2D arrays + try: + import matplotlib.pyplot as plt + # Infer grid shape from coordinate data + unique_x = np.sort(np.unique(current_x)) + unique_y = np.sort(np.unique(current_y)) + nx = unique_x.size + ny = unique_y.size + # Build 2D arrays using coordinate mapping + gz_analytic_2d = np.zeros((nx, ny)) + current_gravity_2d = np.zeros((nx, ny)) + diff_2d = np.zeros((nx, ny)) + x_grid = np.zeros((nx, ny)) + y_grid = np.zeros((nx, ny)) + # Map each station to grid index + for idx in range(len(current_x)): + ix = np.searchsorted(unique_x, current_x[idx]) + iy = np.searchsorted(unique_y, current_y[idx]) + gz_analytic_2d[ix, iy] = gz_analytic[idx] + current_gravity_2d[ix, iy] = current_gravity[idx] + diff_2d[ix, iy] = diff[idx] + x_grid[ix, iy] = current_x[idx] + y_grid[ix, iy] = current_y[idx] + rel_diff_2d = 100.0 * np.abs(diff_2d / (gz_analytic_2d + 1e-16)) # percentage error + # Use min/max for extent + extent = [y_grid.min(), y_grid.max(), x_grid.min(), x_grid.max()] + + # Write VTK surface output + z_val = current_z[0] if len(current_z) > 1 else current_z + write_vtk_surface('gravity_surface.vtk', unique_x, unique_y, z_val, current_gravity_2d, analytic_2d=gz_analytic_2d) + + # Set common colorbar range for analytic and computed + vmin = min(gz_analytic_2d.min(), current_gravity_2d.min()) + vmax = max(gz_analytic_2d.max(), current_gravity_2d.max()) + fig, axs = plt.subplots(1, 3, figsize=(18, 6)) + + # Prism outline coordinates (define here before plotting) + x0, x1, y0, y1, z0, z1 = prism_bounds + prism_x = [y0, y1, y1, y0, y0] + prism_y = [x0, x0, x1, x1, x0] + + im0 = axs[0].imshow(gz_analytic_2d, aspect='equal', cmap='viridis', extent=extent, origin='lower', vmin=vmin, vmax=vmax) + axs[0].set_title('Analytic Solution') + plt.colorbar(im0, ax=axs[0], label='gz (m/s²)') + axs[0].plot(prism_x, prism_y, color='red', linewidth=2, label='Prism Outline') + + im1 = axs[1].imshow(current_gravity_2d, aspect='equal', cmap='viridis', extent=extent, origin='lower', vmin=vmin, vmax=vmax) + axs[1].set_title('Computed (HDF5)') + plt.colorbar(im1, ax=axs[1], label='gz (m/s²)') + axs[1].plot(prism_x, prism_y, color='red', linewidth=2, label='Prism Outline') + + im2 = axs[2].imshow(rel_diff_2d, aspect='equal', cmap='magma', extent=extent, origin='lower') + axs[2].set_title('Relative Difference (%)') + cbar = plt.colorbar(im2, ax=axs[2]) + cbar.set_label('Percent Error (%)') + axs[2].plot(prism_x, prism_y, color='red', linewidth=2, label='Prism Outline') + + for ax in axs: + ax.set_xlabel('Y Coordinate') + ax.set_ylabel('X Coordinate') + ax.set_aspect('equal') + ax.legend() + plt.suptitle('Gravity Comparison: Analytic vs Computed') + plt.tight_layout(rect=[0, 0.03, 1, 0.95]) + if save_fig: + plt.savefig(save_fig, dpi=300) + print(f"Figure saved to {save_fig}") + plt.show() + + # Numerical accuracy test (after plotting) + max_error = np.max(rel_diff_2d) + mean_error = np.mean(rel_diff_2d) + threshold_max = 1.0 # percent + print(f"Max percent error: {max_error:.4f}%") + passed = (max_error < threshold_max) + if passed: + print("Numerical accuracy test: PASS") + else: + print("Numerical accuracy test: FAIL") + return passed + except ImportError: + print("matplotlib not installed. Skipping plot.") + return False + except Exception as e: + print(f"Error reading HDF5 file: {e}") + return False + +if __name__ == "__main__": + if len(sys.argv) < 2: + input_file = "gravity_history.hdf5" + else: + input_file = sys.argv[1] + density_contrast = 900.0 + prism_bounds = [10000, 14000, 8000, 11000, 2000, 2200] + + print(f"Comparing gravity data from: {input_file}") + print(f"Density contrast: {density_contrast}") + print(f"Prism bounds: {prism_bounds}") + save_fig = "gravity_comparison.png" + success = compare_gravity_data(input_file, prism_bounds, density_contrast, save_fig) + if success: + print("Analysis completed successfully.") + else: + print("Analysis failed.") diff --git a/inputFiles/gravityModeling/basicRectPrismC3D8/basicRectPrism_C3D8.xml b/inputFiles/gravityModeling/basicRectPrismC3D8/basicRectPrism_C3D8.xml new file mode 100644 index 00000000000..1578e4122fd --- /dev/null +++ b/inputFiles/gravityModeling/basicRectPrismC3D8/basicRectPrism_C3D8.xml @@ -0,0 +1,1369 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/gravityModeling/basicRectPrismC3D8/extract_gravity_data.py b/inputFiles/gravityModeling/basicRectPrismC3D8/extract_gravity_data.py new file mode 100644 index 00000000000..e68962c096a --- /dev/null +++ b/inputFiles/gravityModeling/basicRectPrismC3D8/extract_gravity_data.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python3 +""" +Extract gravity data from HDF5 file and convert to ASCII format +""" + +import h5py +import numpy as np +import sys + +def extract_gravity_data(hdf5_file, output_file): + """Extract gravity data and coordinates from HDF5 file""" + + try: + with h5py.File(hdf5_file, 'r') as f: + # Read time data + time_data = f['/gravity Time'][:] + + # Read coordinate data + coord_x = f['/coordinateX'][:] + coord_y = f['/coordinateY'][:] + coord_z = f['/coordinateZ'][:] + + # Read gravity data + gravity_data = f['/scatterData'][:] + + print(f"Data shape information:") + print(f" Time steps: {time_data.shape}") + print(f" Coordinates X: {coord_x.shape}") + print(f" Coordinates Y: {coord_y.shape}") + print(f" Coordinates Z: {coord_z.shape}") + print(f" Gravity data: {gravity_data.shape}") + + # Write to ASCII file + with open(output_file, 'w') as out: + # Write header + out.write("# GEOS Gravity Data Export\n") + out.write(f"# Number of time steps: {time_data.shape[0] if len(time_data.shape) > 0 else 1}\n") + out.write(f"# Number of stations: {coord_x.shape[1] if len(coord_x.shape) > 1 else coord_x.shape[0]}\n") + out.write("# Format: Station_ID, X, Y, Z, Time, Gravity\n") + out.write("#\n") + + # Get number of stations + if len(coord_x.shape) > 1: + n_stations = coord_x.shape[1] + n_times = time_data.shape[0] + else: + n_stations = coord_x.shape[0] + n_times = 1 + + # Write data for each time step + for t_idx in range(n_times): + if n_times > 1: + current_time = time_data[t_idx, 0] if len(time_data.shape) > 1 else time_data[t_idx] + current_gravity = gravity_data[t_idx, :] + current_x = coord_x[t_idx, :] + current_y = coord_y[t_idx, :] + current_z = coord_z[t_idx, :] + else: + current_time = time_data[0, 0] if len(time_data.shape) > 1 else time_data[0] + current_gravity = gravity_data[0, :] if len(gravity_data.shape) > 1 else gravity_data + current_x = coord_x[0, :] if len(coord_x.shape) > 1 else coord_x + current_y = coord_y[0, :] if len(coord_y.shape) > 1 else coord_y + current_z = coord_z[0, :] if len(coord_z.shape) > 1 else coord_z + + out.write(f"# Time step {t_idx + 1}, time = {current_time}\n") + + # Write station data + for station_id in range(n_stations): + out.write(f"{station_id:6d} {current_x[station_id]:12.3f} {current_y[station_id]:12.3f} {current_z[station_id]:12.3f} {current_time:12.3f} {current_gravity[station_id]:15.6e}\n") + + out.write("\n") # Blank line between time steps + + print(f"Successfully exported data to: {output_file}") + + # Print summary statistics + print(f"\nSummary:") + print(f" Time range: {np.min(time_data):.3f} to {np.max(time_data):.3f}") + if len(coord_x.shape) > 1: + print(f" X coordinate range: {np.min(coord_x):.3f} to {np.max(coord_x):.3f}") + print(f" Y coordinate range: {np.min(coord_y):.3f} to {np.max(coord_y):.3f}") + print(f" Z coordinate range: {np.min(coord_z):.3f} to {np.max(coord_z):.3f}") + print(f" Gravity range: {np.min(gravity_data):.6e} to {np.max(gravity_data):.6e}") + else: + print(f" X coordinate range: {np.min(current_x):.3f} to {np.max(current_x):.3f}") + print(f" Y coordinate range: {np.min(current_y):.3f} to {np.max(current_y):.3f}") + print(f" Z coordinate range: {np.min(current_z):.3f} to {np.max(current_z):.3f}") + print(f" Gravity range: {np.min(current_gravity):.6e} to {np.max(current_gravity):.6e}") + + except Exception as e: + print(f"Error reading HDF5 file: {e}") + return False + + return True + +if __name__ == "__main__": + if len(sys.argv) < 2: + print("Usage: python extract_gravity_data.py [output.txt]") + sys.exit(1) + + input_file = sys.argv[1] + output_file = sys.argv[2] if len(sys.argv) > 2 else "gravity_data.txt" + + print(f"Extracting gravity data from: {input_file}") + print(f"Output file: {output_file}") + + success = extract_gravity_data(input_file, output_file) + sys.exit(0 if success else 1) From 86353cb8944ab189bdc5f38c79d8bb813ef79199 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Sun, 24 Aug 2025 23:01:27 -0500 Subject: [PATCH 16/30] Fixed clang. --- src/coreComponents/physicsSolvers/gravity/GravityFE.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp index 0141feaed1f..caa0e28757c 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp @@ -47,7 +47,7 @@ class GravityFE : public GravitySolverBase, public ScatterDataProvider GravityFE( GravityFE const & ) = delete; - GravityFE( GravityFE && ) = default; + GravityFE( GravityFE && ) = delete; GravityFE & operator=( GravityFE const & ) = delete; GravityFE & operator=( GravityFE && ) = delete; From 8e47628aae37988a77c03f0f3b1cf0d9bf387c27 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Tue, 26 Aug 2025 18:40:22 -0500 Subject: [PATCH 17/30] Moved scatter to base. --- .../physicsSolvers/gravity/GravityFE.cpp | 42 +++---------------- .../physicsSolvers/gravity/GravityFE.hpp | 13 ++---- .../gravity/GravityFEKernel.hpp | 2 +- .../gravity/GravitySolverBase.cpp | 30 +++++++++++++ .../gravity/GravitySolverBase.hpp | 13 ++++-- 5 files changed, 50 insertions(+), 50 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 204d59246ef..f5e08cd7468 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -59,6 +59,9 @@ void GravityFE::initializePreSubGroups() GEOS_THROW_IF( feDiscretization == nullptr, getName() << ": FE discretization not found: " << m_discretizationName, InputError ); + GEOS_THROW_IF( feDiscretization->getOrder() >1, + getName() << ": FE discretization order should be 1, but we get: " << feDiscretization->getOrder(), + InputError ); } @@ -69,8 +72,6 @@ void GravityFE::registerDataOnMesh( Group & meshBodies ) string_array const & ) { NodeManager & nodeManager = mesh.getNodeManager(); - - nodeManager.registerField< fields::VolumeIntegral >( this->getName() ); ElementRegionManager & elemManager = mesh.getElemManager(); @@ -109,6 +110,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, array1d< real64 > localGzAtStations( m_stationCoordinates.size( 0 ) ); localGzAtStations.setValues< parallelHostPolicy >( 0. ); + auto localGzAtStationsView = localGzAtStations.toView(); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -190,14 +192,12 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, gz += GRAVITATIONAL_CONSTANT * volumeIntegral[a] * dz / r3; } ); - localGzAtStations[iStation]=gz.get(); + localGzAtStationsView[iStation]=gz.get(); } // Loop station } ); // Loop mesh // Brute force reduce... - // Actually using "allReduce" here, since the wrapper for "Reduce" only works on scalars, and not arrays. - // In place: Source==Destination. MpiWrapper::allReduce( localGzAtStations, localGzAtStations, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); arrayView1d< real64 > gzAtStations = m_gzAtStations.toView(); @@ -320,37 +320,7 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, return dt; } -localIndex GravityFE::getNumScatterPoints() const -{ - // For gravity applications: Only rank 0 should provide scatter data - // since all ranks compute the same gravity values at the same stations - if( MpiWrapper::commRank( MPI_COMM_GEOS ) != 0 ) - { - return 0; - } - - return m_stationCoordinates.size( 0 ); // Return number of rows (stations) -} - -array1d< real64 > const & GravityFE::getScatterData() const -{ - GEOS_ASSERT( m_gzAtStations.size() == m_stationCoordinates.size( 0 ) ); - return m_gzAtStations; -} - -array2d< real64 > const & GravityFE::getScatterCoordinates() const -{ - return m_stationCoordinates; -} - -string_array const & GravityFE::getScatterMetadata() const -{ - // Just return empty metadata for now - static string_array empty; - return empty; -} - REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE, string const &, dataRepository::Group * const ) -} /* namespace geos */ +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp index caa0e28757c..76942191586 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp @@ -23,14 +23,12 @@ #include "GravitySolverBase.hpp" #include "mesh/MeshFields.hpp" -#include "fileIO/timeHistory/ScatterDataProvider.hpp" namespace geos { -constexpr real64 GRAVITATIONAL_CONSTANT = 6.67430e-11; // in m³·kg⁻¹·s⁻² -class GravityFE : public GravitySolverBase, public ScatterDataProvider +class GravityFE : public GravitySolverBase { public: @@ -80,11 +78,7 @@ class GravityFE : public GravitySolverBase, public ScatterDataProvider DomainPartition & domain ) override; /**@}*/ - // ScatterDataProvider interface - virtual localIndex getNumScatterPoints() const override; - virtual array1d< real64 > const & getScatterData() const override; - virtual array2d< real64 > const & getScatterCoordinates() const override; - virtual string_array const & getScatterMetadata() const override; + protected: @@ -130,7 +124,6 @@ DECLARE_FIELD( Adjoint, "Adjoint field." ); } - -} /* namespace geos */ +} // namespace geos #endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp index bf2db96bd07..3238f4e154e 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp @@ -39,7 +39,7 @@ struct ForwardVolumeIntegralKernel {} /** - * @brief Launches the precomputation of the DensityVolumeIntegral matrix. + * @brief Launches the precomputation of the VolumeIntegral matrix. * @tparam EXEC_POLICY the execution policy * @tparam ATOMIC_POLICY the atomic policy * @param[in] size Number of elements in the subregion diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp index 3e5027657a2..71ef4808517 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp @@ -169,4 +169,34 @@ real64 GravitySolverBase::explicitStep( real64 const & time_n, } } +localIndex GravitySolverBase::getNumScatterPoints() const +{ + // For gravity applications: Only rank 0 should provide scatter data + // since all ranks compute the same gravity values at the same stations + if( MpiWrapper::commRank( MPI_COMM_GEOS ) != 0 ) + { + return 0; + } + + return m_stationCoordinates.size( 0 ); // Return number of rows (stations) +} + +array1d< real64 > const & GravitySolverBase::getScatterData() const +{ + GEOS_ASSERT( m_gzAtStations.size() == m_stationCoordinates.size( 0 ) ); + return m_gzAtStations; +} + +array2d< real64 > const & GravitySolverBase::getScatterCoordinates() const +{ + return m_stationCoordinates; +} + +string_array const & GravitySolverBase::getScatterMetadata() const +{ + // Just return empty metadata for now + static string_array empty; + return empty; +} + } // namespace geos diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp index 3019be12bc0..fe111f84b29 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp @@ -21,12 +21,15 @@ #define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYSOLVERBASE_HPP_ #include "physicsSolvers/PhysicsSolverBase.hpp" +#include "fileIO/timeHistory/ScatterDataProvider.hpp" namespace geos { +constexpr real64 GRAVITATIONAL_CONSTANT = 6.67430e-11; // in m³·kg⁻¹·s⁻² -class GravitySolverBase : public PhysicsSolverBase + +class GravitySolverBase : public PhysicsSolverBase, public ScatterDataProvider { public: @@ -68,8 +71,6 @@ class GravitySolverBase : public PhysicsSolverBase { static constexpr char const * modeString() { return "mode"; } static constexpr char const * stationCoordinatesString() { return "stationCoordinates"; } - static constexpr char const * outputGzString() { return "outputGz"; } - static constexpr char const * outputGzBasenameString() { return "outputGzBasename"; } static constexpr char const * residueString() { return "residue"; } static constexpr char const * gzAtStationsString() { return "gzAtStations"; } }; @@ -77,6 +78,12 @@ class GravitySolverBase : public PhysicsSolverBase void reinit() override final; + // ScatterDataProvider interface + virtual localIndex getNumScatterPoints() const override; + virtual array1d< real64 > const & getScatterData() const override; + virtual array2d< real64 > const & getScatterCoordinates() const override; + virtual string_array const & getScatterMetadata() const override; + protected: enum class GravityMode { Modeling, Adjoint }; From e0cc3e6780307e2638f1651f607a29a9d2563d5e Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Tue, 26 Aug 2025 19:38:12 -0500 Subject: [PATCH 18/30] Fixed clang. --- src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp index fe111f84b29..dfd206888b5 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp @@ -43,7 +43,7 @@ class GravitySolverBase : public PhysicsSolverBase, public ScatterDataProvider virtual ~GravitySolverBase() override; GravitySolverBase( GravitySolverBase const & ) = delete; - GravitySolverBase( GravitySolverBase && ) = default; + GravitySolverBase( GravitySolverBase && ) = delete; GravitySolverBase & operator=( GravitySolverBase const & ) = delete; GravitySolverBase & operator=( GravitySolverBase && ) = delete; From 0262245c3bbbaa4a5e304cf6a2fb1b5d84e62fac Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 27 Aug 2025 10:33:06 -0500 Subject: [PATCH 19/30] Added view. --- src/coreComponents/physicsSolvers/gravity/GravityFE.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index f5e08cd7468..fb74343b870 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -174,7 +174,9 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, // Hardcode lowest order... i.e. take advantage that mesh vertices and quadrature nodes are collocated. // Warning: this code will fail for higher order. - auto const & nodePosition = nodeManager.referencePosition(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = + nodeManager.referencePosition().toViewConst(); + for( localIndex iStation=0; iStation const nodePosition = + nodeManager.referencePosition().toViewConst(); elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & elementSubRegion ) From 98a856fee51b53e56ec13d08f11c596324892ab2 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 27 Aug 2025 11:30:25 -0500 Subject: [PATCH 20/30] Removed brackets --- .../physicsSolvers/gravity/GravityFE.cpp | 40 ++++++++++++------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index fb74343b870..1638d79bf2b 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -176,21 +176,25 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, // Warning: this code will fail for higher order. arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition().toViewConst(); + auto const station = m_stationCoordinates.toViewConst(); for( localIndex iStation=0; iStation gz( 0.0 ); forAll< parallelDevicePolicy<> >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) { - real64 dx = nodePosition[a][0] - coords[0]; - real64 dy = nodePosition[a][1] - coords[1]; - real64 dz = nodePosition[a][2] - coords[2]; - real64 r2 = dx*dx + dy*dy + dz*dz; - real64 r3 = std::sqrt( r2 ) * r2; + const real64 dx = nodePosition( a, 0 ) - cx; + const real64 dy = nodePosition( a, 1 ) - cy; + const real64 dz = nodePosition( a, 2 ) - cz; + const real64 r2 = dx*dx + dy*dy + dz*dz; + const real64 r3 = std::sqrt( r2 ) * r2; gz += GRAVITATIONAL_CONSTANT * volumeIntegral[a] * dz / r3; } ); @@ -256,10 +260,14 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, { arrayView1d< real64 const > const residue = m_residue.toViewConst(); + auto const station = m_stationCoordinates.toViewConst(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); arrayView2d< real64 > volumeIntegral2d = elementSubRegion.getReference< array2d< real64 > >( fields::VolumeIntegral2d::key()); arrayView1d< real64 > adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()); + + volumeIntegral2d.setValues< parallelHostPolicy >( 0.0 ); adjoint.zero(); arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition().toViewConst(); @@ -288,21 +296,23 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, for( localIndex iStation=0; iStation( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) { for( localIndex iLoc=0; iLoc < numSupportPoints; ++iLoc ) { - localIndex a = elemsToNodes[k][iLoc]; - - real64 dx = nodePosition[a][0] - coords[0]; - real64 dy = nodePosition[a][1] - coords[1]; - real64 dz = nodePosition[a][2] - coords[2]; - real64 r2 = dx*dx + dy*dy + dz*dz; - real64 r3 = sqrt( r2 ) * r2; - adjoint[k] += GRAVITATIONAL_CONSTANT * volumeIntegral2d[k][iLoc] * res * dz / r3; + const localIndex a = elemsToNodes( k, iLoc ); + const real64 dx = nodePosition( a, 0 ) - cx; + const real64 dy = nodePosition( a, 1 ) - cy; + const real64 dz = nodePosition( a, 2 ) - cz; + const real64 r2 = dx*dx + dy*dy + dz*dz; + const real64 r3 = std::sqrt( r2 ) * r2; + adjoint[k] += GRAVITATIONAL_CONSTANT * volumeIntegral2d( k, iLoc ) * res * dz / r3; } } ); // Loop elem } //Loop station From 3000036ddfb7168f58ab146963f32f2f1b5c2da4 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 27 Aug 2025 12:35:26 -0500 Subject: [PATCH 21/30] Addressed CI --- .../physicsSolvers/gravity/GravityFE.cpp | 25 +++++++++++++------ 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 1638d79bf2b..699a9ddd702 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -119,10 +119,10 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, { // Step 1: Precompute volumeIntegral once and for all. NodeManager & nodeManager = mesh.getNodeManager(); - arrayView1d< integer const > const & nodeGhostRank = nodeManager.ghostRank(); + auto const nodeGhostRank = nodeManager.ghostRank().toViewConst(); // VolumeIntegral matrix to be computed in this function. - arrayView1d< real64 > const volumeIntegral = nodeManager.getField< fields::VolumeIntegral >(); + auto volumeIntegral = nodeManager.getField< fields::VolumeIntegral >().toView(); volumeIntegral.setValues< parallelHostPolicy >( 0. ); // Loop over all sub-regions in regions of type "CellElements". @@ -194,8 +194,12 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, const real64 dy = nodePosition( a, 1 ) - cy; const real64 dz = nodePosition( a, 2 ) - cz; const real64 r2 = dx*dx + dy*dy + dz*dz; - const real64 r3 = std::sqrt( r2 ) * r2; - gz += GRAVITATIONAL_CONSTANT * volumeIntegral[a] * dz / r3; + real64 inv_r3 = 0.0; + if( r2 > 0.0 ) + { + inv_r3 = 1.0 / ( std::sqrt( r2 ) * r2 ); + } + gz += GRAVITATIONAL_CONSTANT * volumeIntegral[a] * dz * inv_r3; } ); localGzAtStationsView[iStation]=gz.get(); @@ -264,8 +268,8 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); - arrayView2d< real64 > volumeIntegral2d = elementSubRegion.getReference< array2d< real64 > >( fields::VolumeIntegral2d::key()); - arrayView1d< real64 > adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()); + arrayView2d< real64 > volumeIntegral2d = elementSubRegion.getReference< array2d< real64 > >( fields::VolumeIntegral2d::key()).toView(); + arrayView1d< real64 > adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()).toView(); volumeIntegral2d.setValues< parallelHostPolicy >( 0.0 ); adjoint.zero(); @@ -311,8 +315,13 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, const real64 dy = nodePosition( a, 1 ) - cy; const real64 dz = nodePosition( a, 2 ) - cz; const real64 r2 = dx*dx + dy*dy + dz*dz; - const real64 r3 = std::sqrt( r2 ) * r2; - adjoint[k] += GRAVITATIONAL_CONSTANT * volumeIntegral2d( k, iLoc ) * res * dz / r3; + real64 inv_r3 = 0.0; + if( r2 > 0.0 ) + { + inv_r3 = 1.0 / ( std::sqrt( r2 ) * r2 ); + } + + adjoint[k] += GRAVITATIONAL_CONSTANT * volumeIntegral2d( k, iLoc ) * res * dz * inv_r3; } } ); // Loop elem } //Loop station From 730b4d40607554e82b7e542ae78612ff9034dffc Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 27 Aug 2025 13:44:44 -0500 Subject: [PATCH 22/30] Added coupling --- .../physicsSolvers/gravity/CMakeLists.txt | 2 + .../physicsSolvers/gravity/GravityFE.cpp | 70 ++-- .../GravityFE_CompositionalMultiphaseFVM.cpp | 386 ++++++++++++++++++ .../GravityFE_CompositionalMultiphaseFVM.hpp | 150 +++++++ 4 files changed, 579 insertions(+), 29 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp create mode 100644 src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp diff --git a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt index a05504c0211..2b7c1d7eae8 100644 --- a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt @@ -20,11 +20,13 @@ Contains: set( gravitySolvers_headers GravitySolverBase.hpp GravityFE.hpp + GravityFE_CompositionalMultiphaseFVM.hpp GravityFEKernel.hpp ) # Specify solver sources set( gravitySolvers_sources GravitySolverBase.cpp + GravityFE_CompositionalMultiphaseFVM.cpp GravityFE.cpp ) set( dependencyList ${parallelDeps} physicsSolversBase ) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 699a9ddd702..c8237389b5a 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -123,6 +123,8 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, // VolumeIntegral matrix to be computed in this function. auto volumeIntegral = nodeManager.getField< fields::VolumeIntegral >().toView(); + + volumeIntegral.setValues< parallelHostPolicy >( 0. ); // Loop over all sub-regions in regions of type "CellElements". @@ -171,13 +173,19 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, // Step #2: Compute contribution to all stations. - + // Make volumeIntegral const + arrayView1d< real64 const > const volumeIntegralConst = volumeIntegral.toViewConst(); + // Hardcode lowest order... i.e. take advantage that mesh vertices and quadrature nodes are collocated. // Warning: this code will fail for higher order. arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition().toViewConst(); auto const station = m_stationCoordinates.toViewConst(); + // Constants + constexpr real64 GRAVITATIONAL_CONSTANT_LOCAL = GRAVITATIONAL_CONSTANT; + constexpr real64 EPSILON = 1e-15; + for( localIndex iStation=0; iStation 0.0 ) + if( r2 > EPSILON ) { - inv_r3 = 1.0 / ( std::sqrt( r2 ) * r2 ); + inv_r3 = 1.0 / ( r2 * std::sqrt( r2 ) ); } - gz += GRAVITATIONAL_CONSTANT * volumeIntegral[a] * dz * inv_r3; + gz += GRAVITATIONAL_CONSTANT_LOCAL * volumeIntegralConst[a] * dz * inv_r3; } ); localGzAtStationsView[iStation]=gz.get(); @@ -248,7 +256,6 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, array1d< real64 > localGzAtStations( m_stationCoordinates.size( 0 ) ); localGzAtStations.setValues< parallelHostPolicy >( 0. ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, string_array const & regionNames ) @@ -262,27 +269,26 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & elementSubRegion ) { - + // Const views for input data arrayView1d< real64 const > const residue = m_residue.toViewConst(); - auto const station = m_stationCoordinates.toViewConst(); - + auto const station = m_stationCoordinates.toViewConst(); + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition().toViewConst(); - arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes = elementSubRegion.nodeList(); + // Output arrays arrayView2d< real64 > volumeIntegral2d = elementSubRegion.getReference< array2d< real64 > >( fields::VolumeIntegral2d::key()).toView(); arrayView1d< real64 > adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()).toView(); volumeIntegral2d.setValues< parallelHostPolicy >( 0.0 ); adjoint.zero(); - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition().toViewConst(); - finiteElement::FiniteElementBase const & fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); localIndex const numSupportPoints = fe.getNumSupportPoints(); + localIndex const numStations = m_stationCoordinates.size( 0 ); GEOS_ERROR_IF_GT_MSG( numSupportPoints, MAX_SUPPORT_POINTS, GEOS_FMT( "Maximum number of SupportPoints is {}", MAX_SUPPORT_POINTS ) ); - finiteElement::dispatchlowOrder3D( fe, [&] ( auto const finiteElement ) { using FE_TYPE = TYPEOFREF( finiteElement ); @@ -296,19 +302,26 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, volumeIntegral2d ); } ); + // Make volumeIntegral2d const + arrayView2d< real64 const > const volumeIntegral2dConst = volumeIntegral2d.toViewConst(); - for( localIndex iStation=0; iStation( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + // Optimized memory access: compute all stations for each element in one kernel + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + adjoint[k] = 0.0; // Initialize once + + for( localIndex iStation = 0; iStation < numStations; ++iStation ) { - for( localIndex iLoc=0; iLoc < numSupportPoints; ++iLoc ) + const real64 cx = station( iStation, 0 ); + const real64 cy = station( iStation, 1 ); + const real64 cz = station( iStation, 2 ); + const real64 res = residue[iStation]; + + for( localIndex iLoc = 0; iLoc < numSupportPoints; ++iLoc ) { const localIndex a = elemsToNodes( k, iLoc ); const real64 dx = nodePosition( a, 0 ) - cx; @@ -316,15 +329,15 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, const real64 dz = nodePosition( a, 2 ) - cz; const real64 r2 = dx*dx + dy*dy + dz*dz; real64 inv_r3 = 0.0; - if( r2 > 0.0 ) + if( r2 > EPSILON ) { - inv_r3 = 1.0 / ( std::sqrt( r2 ) * r2 ); + inv_r3 = 1.0 / ( r2 * std::sqrt( r2 ) ); } - - adjoint[k] += GRAVITATIONAL_CONSTANT * volumeIntegral2d( k, iLoc ) * res * dz * inv_r3; + + adjoint[k] += GRAVITATIONAL_CONSTANT_LOCAL * volumeIntegral2dConst( k, iLoc ) * res * dz * inv_r3; } - } ); // Loop elem - } //Loop station + } + } ); // Loop elem adjoint.move( LvArray::MemorySpace::host, true ); @@ -342,7 +355,6 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, return dt; } - REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE, string const &, dataRepository::Group * const ) } // namespace geos diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp new file mode 100644 index 00000000000..b6c228f4c19 --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp @@ -0,0 +1,386 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 GravityFE_CompositionalMultiphaseFVM.cpp + */ +#include + +#include "GravityFE_CompositionalMultiphaseFVM.hpp" +#include "GravityFEKernel.hpp" + + +#include "finiteElement/FiniteElementDiscretization.hpp" +#include "mesh/DomainPartition.hpp" +#include "discretizationMethods/NumericalMethodsManager.hpp" +#include "finiteElement/FiniteElementDiscretizationManager.hpp" +#include "common/MpiWrapper.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/ElementType.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" + +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/solid/CoupledSolid.hpp" +#include "constitutive/solid/SolidBase.hpp" +#include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" + +namespace geos +{ + +using namespace constitutive; + + +using namespace dataRepository; + +GravityFE_CompositionalMultiphaseFVM::GravityFE_CompositionalMultiphaseFVM( const std::string & name, + Group * const parent ): + GravitySolverBase( name, parent ) +{ + + registerWrapper( viewKeyStruct::flowSolverNameString(), &m_flowSolverName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the flow solver to dynamically compute the density values" ); + + registerWrapper( viewKeyStruct::useRockDensityString(), &m_useRockDensity ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1 ). + setDescription( "Flag to use the rock density along with the fluid density" ); + + registerWrapper( viewKeyStruct::usePorosityString(), &m_usePorosity ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "Flag to use the true rock porosity (changing with pressure, etc.)" ); + + registerWrapper( viewKeyStruct::useReferencePorosityString(), &m_useReferencePorosity ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1 ). + setDescription( "Flag to use the reference porosity" ); +} + +GravityFE_CompositionalMultiphaseFVM::~GravityFE_CompositionalMultiphaseFVM() +{ + // TODO Auto-generated destructor stub +} + + +void GravityFE_CompositionalMultiphaseFVM::initializePreSubGroups() +{ + GravitySolverBase::initializePreSubGroups(); + DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteElementDiscretizationManager const & + feDiscretizationManager = numericalMethodManager.getFiniteElementDiscretizationManager(); + + FiniteElementDiscretization const * const + feDiscretization = feDiscretizationManager.getGroupPointer< FiniteElementDiscretization >( m_discretizationName ); + GEOS_THROW_IF( feDiscretization == nullptr, + getName() << ": FE discretization not found: " << m_discretizationName, + InputError ); + GEOS_THROW_IF( feDiscretization->getOrder() >1, + getName() << ": FE discretization order should be 1, but we get: " << feDiscretization->getOrder(), + InputError ); +} + + +void GravityFE_CompositionalMultiphaseFVM::registerDataOnMesh( Group & meshBodies ) +{ + + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + const string_array & ) + { + NodeManager & nodeManager = mesh.getNodeManager(); + + nodeManager.registerField< fields::VolumeIntegral >( this->getName() ); + + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion ) + { + subRegion.registerField< fields::MediumDensity >( this->getName() ); + subRegion.registerField< fields::FluidDensity >( this->getName() ); + subRegion.registerField< fields::RockDensity >( this->getName() ); + subRegion.registerField< fields::Porosity >( this->getName() ); + } ); + } ); +} + + +void GravityFE_CompositionalMultiphaseFVM::postInputInitialization() +{ + GravitySolverBase::postInputInitialization(); + + // Create a pointer to the flow solver. + m_flowSolver = &this->getParent().getGroup< CompositionalMultiphaseFVM >( m_flowSolverName ); +} + + + +void GravityFE_CompositionalMultiphaseFVM::initializePostInitialConditionsPreSubGroups() +{ + GravitySolverBase::initializePostInitialConditionsPreSubGroups(); +} + + + +real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( time_n, cycleNumber ); + + array1d< real64 > localGzAtStations( m_stationCoordinates.size( 0 ) ); + localGzAtStations.setValues< parallelHostPolicy >( 0. ); + auto localGzAtStationsView = localGzAtStations.toView(); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), + [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + // Step 1: Precompute volumeIntegral once and for all. + NodeManager & nodeManager = mesh.getNodeManager(); + auto const nodeGhostRank = nodeManager.ghostRank().toViewConst(); + + // VolumeIntegral matrix to be computed in this function. + auto volumeIntegral = nodeManager.getField< fields::VolumeIntegral >().toView(); + volumeIntegral.setValues< parallelHostPolicy >( 0. ); + + // Loop over all sub-regions in regions of type "CellElements". + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + CellElementSubRegion & elementSubRegion ) + { + Group const & constitutiveModels = elementSubRegion.getGroup( ConstitutiveManager::groupKeyStruct::constitutiveModelsString() ); + + // Retrieve fluid density. + string const & fluidName = elementSubRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::fluidNamesString() ); + MultiFluidBase const & fluid = constitutiveModels.getGroup< MultiFluidBase >( fluidName ); + arrayView2d< real64 const, multifluid::USD_FLUID > const reservoirFluidDensity = fluid.totalDensity().toViewConst(); + + arrayView1d< real64 > const fluidDensity = elementSubRegion.getReference< array1d< real64 > >( fields::FluidDensity::key()).toView(); + + // Fluid density assignment + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + fluidDensity[i] = reservoirFluidDensity[i][0]; + } ); + + // Retrieve porosity. + string const & solidName = elementSubRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); + arrayView1d< real64 > const porosity = elementSubRegion.getReference< array1d< real64 > >( fields::Porosity::key()).toView(); + + if( this->m_useReferencePorosity == 1 ) + { + GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Use referencePorosity" ); + arrayView1d< real64 const > const reservoirPorosity = solid.getReferencePorosity().toViewConst(); + + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + porosity[i] = reservoirPorosity[i]; + } ); + } + else if( this->m_usePorosity == 1 ) + { + GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Use Porosity" ); + arrayView2d< real64 const > const reservoirPorosity = solid.getPorosity().toViewConst(); + + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + porosity[i] = reservoirPorosity[i][0]; + } ); + } + else + { + GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Set Porosity to 1" ); + + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + porosity[i] = 1.0; + } ); + } + + // Retrieve rock density. + arrayView1d< real64 > const rockDensity = elementSubRegion.getReference< array1d< real64 > >( fields::RockDensity::key()).toView(); + rockDensity.setValues< parallelHostPolicy >( 0. ); + + if( this->m_useRockDensity == 1 ) + { + string const & solidModelName = elementSubRegion.getReference< string >( SolidMechanicsLagrangianFEM::viewKeyStruct::solidMaterialNamesString()); + arrayView2d< real64 const > const reservoirRockDensity = elementSubRegion.getConstitutiveModel( solidModelName ).getReference< array2d< real64 > >( SolidBase::viewKeyStruct::densityString() ).toViewConst(); + + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + rockDensity[i] = reservoirRockDensity[i][0]; + } ); + } + + // Build medium density. + arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()).toView(); + + // Make input arrays const for density calculation + arrayView1d< real64 const > const porosityConst = porosity.toViewConst(); + arrayView1d< real64 const > const fluidDensityConst = fluidDensity.toViewConst(); + arrayView1d< real64 const > const rockDensityConst = rockDensity.toViewConst(); + + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + density[i] = porosityConst[i] * fluidDensityConst[i] + (1.0 - porosityConst[i]) * rockDensityConst[i]; + } ); + + // Debug output + if( this->getLogLevel() > 3 ) + { + // Move arrays to host for logging + density.move( LvArray::MemorySpace::host, true ); + fluidDensityConst.move( LvArray::MemorySpace::host, true ); + rockDensityConst.move( LvArray::MemorySpace::host, true ); + porosityConst.move( LvArray::MemorySpace::host, true ); + + for( localIndex i = 0; i < elementSubRegion.size(); ++i ) + { + GEOS_LOG( "GravityFE_CompositionalMultiphaseFVM: Cell[" << i << "], density= " << density[i] + << ", fluidDensity= " << fluidDensityConst[i] << ", rockDensity= " << rockDensityConst[i] + << ", porosity= " << porosityConst[i] ); + } + } + + // Make density const after computation + arrayView1d< real64 const > const densityConst = density.toViewConst(); + + // Compute volume integral at each mesh node. + arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition().toViewConst(); + + finiteElement::FiniteElementBase const & + fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() ); + + finiteElement::dispatchlowOrder3D( fe, [&] ( auto const finiteElement ) + { + using FE_TYPE = TYPEOFREF( finiteElement ); + + gravityFEKernel:: + ForwardVolumeIntegralKernel< FE_TYPE > kernel( finiteElement ); + kernel.template launch< EXEC_POLICY, ATOMIC_POLICY > + ( elementSubRegion.size(), + X, + elemsToNodes, + densityConst, + volumeIntegral ); + } ); + + // Zero out ghost node. + forAll< parallelDevicePolicy<> >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + if( nodeGhostRank[a] >= 0 ) + { + volumeIntegral[a] = 0.; + } + } ); + } ); // loop on cellElements + + + // Step #2: Compute contribution to all stations. + + // Make volumeIntegral const + arrayView1d< real64 const > const volumeIntegralConst = volumeIntegral.toViewConst(); + // Hardcode lowest order... i.e. take advantage that mesh vertices and quadrature nodes are collocated. + // Warning: this code will fail for higher order. + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = + nodeManager.referencePosition().toViewConst(); + auto const station = m_stationCoordinates.toViewConst(); + + // Constants + constexpr real64 GRAVITATIONAL_CONSTANT_LOCAL = GRAVITATIONAL_CONSTANT; + constexpr real64 EPSILON = 1e-15; + + for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) + { + // Deal with one station. + const real64 cx = station( iStation, 0 ); + const real64 cy = station( iStation, 1 ); + const real64 cz = station( iStation, 2 ); + + RAJA::ReduceSum< parallelDeviceReduce, real64 > gz( 0.0 ); + + forAll< parallelDevicePolicy<> >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + const real64 dx = nodePosition( a, 0 ) - cx; + const real64 dy = nodePosition( a, 1 ) - cy; + const real64 dz = nodePosition( a, 2 ) - cz; + const real64 r2 = dx*dx + dy*dy + dz*dz; + real64 inv_r3 = 0.0; + if( r2 > EPSILON ) + { + inv_r3 = 1.0 / ( r2 * std::sqrt( r2 ) ); + } + gz += GRAVITATIONAL_CONSTANT_LOCAL * volumeIntegralConst[a] * dz * inv_r3; + } ); + + localGzAtStationsView[iStation] = gz.get(); + } // Loop station + } ); // Loop mesh + + // Brute force reduce... + MpiWrapper::allReduce( localGzAtStations, localGzAtStations, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); + + arrayView1d< real64 > gzAtStations = m_gzAtStations.toView(); + for( localIndex i = 0; i < gzAtStations.size(); ++i ) + { + gzAtStations[i] = localGzAtStations[i]; + } + + if( this->getLogLevel() > 1 ) + { + for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) + { + auto const & coords = m_stationCoordinates[iStation]; + std::ostringstream logStream; + logStream << std::fixed << std::setprecision( 2 ); + logStream << "GravityFE_CompositionalMultiphaseFVM: station[" << std::setw( 5 ) << iStation << "] " + << std::setw( 15 ) << coords[0] << " " + << std::setw( 15 ) << coords[1] << " " + << std::setw( 10 ) << coords[2] << " " + << std::scientific << std::setprecision( 6 ) + << std::setw( 14 ) << gzAtStations[iStation]; + GEOS_LOG_RANK_0( logStream.str() ); + } + } + + return dt; +} + + +real64 GravityFE_CompositionalMultiphaseFVM::explicitStepAdjoint( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) +{ + GEOS_UNUSED_VAR( time_n, cycleNumber, domain ); + GEOS_ERROR( "Adjoint computation is not available when Gravity is coupled with Flow simulator." ); + return dt; +} + + + +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE_CompositionalMultiphaseFVM, string const &, dataRepository::Group * const ) + +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp new file mode 100644 index 00000000000..e72de7577b6 --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp @@ -0,0 +1,150 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 GravityFE_CompositionalMultiphaseFVM.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ +#define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ + +#include "GravitySolverBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" + +namespace geos +{ + +class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase +{ + +public: + using EXEC_POLICY = parallelDevicePolicy< 32 >; + using ATOMIC_POLICY = parallelDeviceAtomic; + + GravityFE_CompositionalMultiphaseFVM() = delete; + GravityFE_CompositionalMultiphaseFVM( const std::string & name, + Group * const parent ); + virtual ~GravityFE_CompositionalMultiphaseFVM() override; + + GravityFE_CompositionalMultiphaseFVM( GravityFE_CompositionalMultiphaseFVM const & ) = delete; + GravityFE_CompositionalMultiphaseFVM( GravityFE_CompositionalMultiphaseFVM && ) = delete; + GravityFE_CompositionalMultiphaseFVM & operator=( GravityFE_CompositionalMultiphaseFVM const & ) = delete; + GravityFE_CompositionalMultiphaseFVM & operator=( GravityFE_CompositionalMultiphaseFVM && ) = delete; + + static string catalogName() { return "GravityFE_CompositionalMultiphaseFVM"; } + string getCatalogName() const override { return catalogName(); } + + virtual void initializePreSubGroups() override; + virtual void registerDataOnMesh( Group & meshBodies ) override final; + virtual void postInputInitialization() override final; + + + /** + * @defgroup Solver Interface Functions + * + * These functions provide the primary interface that is required for derived classes + */ + /**@{*/ + virtual + real64 explicitStepModeling( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override; + + virtual + real64 explicitStepAdjoint( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override; + + + + /**@}*/ + + struct viewKeyStruct : GravitySolverBase::viewKeyStruct + { + static constexpr char const * useRockDensityString() { return "useRockDensity"; } + static constexpr char const * usePorosityString() { return "usePorosity"; } + static constexpr char const * useReferencePorosityString() { return "useReferencePorosity"; } + constexpr static char const * flowSolverNameString() { return "flowSolverName"; } + } gravityViewKeys; + + +protected: + + + + virtual void initializePostInitialConditionsPreSubGroups() override final; + + /// pointer to the flow solver + CompositionalMultiphaseFVM * m_flowSolver; + string m_flowSolverName; + + /// Use rock density in addition to fluid density + localIndex m_useRockDensity; + localIndex m_useReferencePorosity; + localIndex m_usePorosity; + +private: + + + +}; + +namespace fields +{ + +DECLARE_FIELD( MediumDensity, + "mediumDensity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Medium density of the cell" ); +DECLARE_FIELD( FluidDensity, + "fluidDensity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Fluid density of the cell" ); +DECLARE_FIELD( RockDensity, + "rockDensity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Rock density of the cell" ); +DECLARE_FIELD( Porosity, + "porosity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Porosity of the cell" ); +DECLARE_FIELD( VolumeIntegral, + "volumeIntegral", + array1d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "VolumeIntegral of the cell." ); + +} + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ From e61558842d2cac7be4c129c990156480a1f18c73 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 27 Aug 2025 15:06:55 -0500 Subject: [PATCH 23/30] Fixed uncrustify and dependency --- .../physicsSolvers/gravity/CMakeLists.txt | 3 +-- .../physicsSolvers/gravity/GravityFE.cpp | 8 +++--- .../GravityFE_CompositionalMultiphaseFVM.cpp | 27 ++++++++++--------- .../GravityFE_CompositionalMultiphaseFVM.hpp | 4 +-- 4 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt index 2b7c1d7eae8..cf81fb71290 100644 --- a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt @@ -29,8 +29,7 @@ set( gravitySolvers_sources GravityFE_CompositionalMultiphaseFVM.cpp GravityFE.cpp ) -set( dependencyList ${parallelDeps} physicsSolversBase ) - +set( dependencyList ${parallelDeps} physicsSolversBase fluidFlowSolvers) geos_decorate_link_dependencies( LIST decoratedDependencies DEPENDENCIES ${dependencyList} ) blt_add_library( NAME gravitySolvers diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index c8237389b5a..9eb4fd1b200 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -175,7 +175,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, // Step #2: Compute contribution to all stations. // Make volumeIntegral const arrayView1d< real64 const > const volumeIntegralConst = volumeIntegral.toViewConst(); - + // Hardcode lowest order... i.e. take advantage that mesh vertices and quadrature nodes are collocated. // Warning: this code will fail for higher order. arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = @@ -313,14 +313,14 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) { adjoint[k] = 0.0; // Initialize once - + for( localIndex iStation = 0; iStation < numStations; ++iStation ) { const real64 cx = station( iStation, 0 ); const real64 cy = station( iStation, 1 ); const real64 cz = station( iStation, 2 ); const real64 res = residue[iStation]; - + for( localIndex iLoc = 0; iLoc < numSupportPoints; ++iLoc ) { const localIndex a = elemsToNodes( k, iLoc ); @@ -333,7 +333,7 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, { inv_r3 = 1.0 / ( r2 * std::sqrt( r2 ) ); } - + adjoint[k] += GRAVITATIONAL_CONSTANT_LOCAL * volumeIntegral2dConst( k, iLoc ) * res * dz * inv_r3; } } diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp index b6c228f4c19..2e266f06b7b 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp @@ -94,7 +94,7 @@ void GravityFE_CompositionalMultiphaseFVM::initializePreSubGroups() InputError ); GEOS_THROW_IF( feDiscretization->getOrder() >1, getName() << ": FE discretization order should be 1, but we get: " << feDiscretization->getOrder(), - InputError ); + InputError ); } @@ -192,7 +192,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const { GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Use referencePorosity" ); arrayView1d< real64 const > const reservoirPorosity = solid.getReferencePorosity().toViewConst(); - + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { porosity[i] = reservoirPorosity[i]; @@ -202,7 +202,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const { GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Use Porosity" ); arrayView2d< real64 const > const reservoirPorosity = solid.getPorosity().toViewConst(); - + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { porosity[i] = reservoirPorosity[i][0]; @@ -211,7 +211,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const else { GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Set Porosity to 1" ); - + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { porosity[i] = 1.0; @@ -221,12 +221,13 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const // Retrieve rock density. arrayView1d< real64 > const rockDensity = elementSubRegion.getReference< array1d< real64 > >( fields::RockDensity::key()).toView(); rockDensity.setValues< parallelHostPolicy >( 0. ); - + if( this->m_useRockDensity == 1 ) { string const & solidModelName = elementSubRegion.getReference< string >( SolidMechanicsLagrangianFEM::viewKeyStruct::solidMaterialNamesString()); - arrayView2d< real64 const > const reservoirRockDensity = elementSubRegion.getConstitutiveModel( solidModelName ).getReference< array2d< real64 > >( SolidBase::viewKeyStruct::densityString() ).toViewConst(); - + arrayView2d< real64 const > const reservoirRockDensity = + elementSubRegion.getConstitutiveModel( solidModelName ).getReference< array2d< real64 > >( SolidBase::viewKeyStruct::densityString() ).toViewConst(); + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { rockDensity[i] = reservoirRockDensity[i][0]; @@ -235,12 +236,12 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const // Build medium density. arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()).toView(); - + // Make input arrays const for density calculation arrayView1d< real64 const > const porosityConst = porosity.toViewConst(); arrayView1d< real64 const > const fluidDensityConst = fluidDensity.toViewConst(); arrayView1d< real64 const > const rockDensityConst = rockDensity.toViewConst(); - + forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { density[i] = porosityConst[i] * fluidDensityConst[i] + (1.0 - porosityConst[i]) * rockDensityConst[i]; @@ -254,12 +255,12 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const fluidDensityConst.move( LvArray::MemorySpace::host, true ); rockDensityConst.move( LvArray::MemorySpace::host, true ); porosityConst.move( LvArray::MemorySpace::host, true ); - + for( localIndex i = 0; i < elementSubRegion.size(); ++i ) { - GEOS_LOG( "GravityFE_CompositionalMultiphaseFVM: Cell[" << i << "], density= " << density[i] - << ", fluidDensity= " << fluidDensityConst[i] << ", rockDensity= " << rockDensityConst[i] - << ", porosity= " << porosityConst[i] ); + GEOS_LOG( "GravityFE_CompositionalMultiphaseFVM: Cell[" << i << "], density= " << density[i] + << ", fluidDensity= " << fluidDensityConst[i] << ", rockDensity= " << rockDensityConst[i] + << ", porosity= " << porosityConst[i] ); } } diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp index e72de7577b6..eaebfed5d51 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp @@ -36,7 +36,7 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase GravityFE_CompositionalMultiphaseFVM() = delete; GravityFE_CompositionalMultiphaseFVM( const std::string & name, - Group * const parent ); + Group * const parent ); virtual ~GravityFE_CompositionalMultiphaseFVM() override; GravityFE_CompositionalMultiphaseFVM( GravityFE_CompositionalMultiphaseFVM const & ) = delete; @@ -147,4 +147,4 @@ DECLARE_FIELD( VolumeIntegral, } // namespace geos -#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ +#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ From cfaa28634e704e0670ead7aeca20fcb92bc40363 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 27 Aug 2025 15:35:19 -0500 Subject: [PATCH 24/30] Made flow conditional --- .../physicsSolvers/gravity/CMakeLists.txt | 16 +++++++++++++--- .../GravityFE_CompositionalMultiphaseFVM.cpp | 3 ++- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt index cf81fb71290..c4663d1824d 100644 --- a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt @@ -20,16 +20,26 @@ Contains: set( gravitySolvers_headers GravitySolverBase.hpp GravityFE.hpp - GravityFE_CompositionalMultiphaseFVM.hpp GravityFEKernel.hpp ) # Specify solver sources set( gravitySolvers_sources GravitySolverBase.cpp - GravityFE_CompositionalMultiphaseFVM.cpp GravityFE.cpp ) -set( dependencyList ${parallelDeps} physicsSolversBase fluidFlowSolvers) +set( dependencyList ${parallelDeps} physicsSolversBase ) + +# Conditionally add compositional multiphase support if fluid flow is enabled +if( GEOS_ENABLE_FLUIDFLOW AND TARGET fluidFlowSolvers ) + list( APPEND gravitySolvers_headers GravityFE_CompositionalMultiphaseFVM.hpp ) + list( APPEND gravitySolvers_sources GravityFE_CompositionalMultiphaseFVM.cpp ) + list( APPEND dependencyList fluidFlowSolvers ) + + # Define a preprocessor macro so the code knows fluid flow is available + add_definitions( -DGEOS_GRAVITY_WITH_FLUIDFLOW ) +endif() + + geos_decorate_link_dependencies( LIST decoratedDependencies DEPENDENCIES ${dependencyList} ) blt_add_library( NAME gravitySolvers diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp index 2e266f06b7b..e0afe68472a 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp @@ -381,7 +381,8 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepAdjoint( real64 const & } - +#ifdef GEOS_GRAVITY_WITH_FLUIDFLOW REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE_CompositionalMultiphaseFVM, string const &, dataRepository::Group * const ) +#endif } // namespace geos From ac613ad151d948914f24fe6f75420e77861bb71c Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 18 Sep 2025 14:34:08 -0500 Subject: [PATCH 25/30] clean pt --- .../gravity/GravityFE_CompositionalMultiphaseFVM.cpp | 5 ----- .../gravity/GravityFE_CompositionalMultiphaseFVM.hpp | 5 ----- 2 files changed, 10 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp index e0afe68472a..341051555fb 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp @@ -52,9 +52,6 @@ GravityFE_CompositionalMultiphaseFVM::GravityFE_CompositionalMultiphaseFVM( cons GravitySolverBase( name, parent ) { - registerWrapper( viewKeyStruct::flowSolverNameString(), &m_flowSolverName ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Name of the flow solver to dynamically compute the density values" ); registerWrapper( viewKeyStruct::useRockDensityString(), &m_useRockDensity ). setInputFlag( InputFlags::OPTIONAL ). @@ -126,8 +123,6 @@ void GravityFE_CompositionalMultiphaseFVM::postInputInitialization() { GravitySolverBase::postInputInitialization(); - // Create a pointer to the flow solver. - m_flowSolver = &this->getParent().getGroup< CompositionalMultiphaseFVM >( m_flowSolverName ); } diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp index eaebfed5d51..4ce0d908803 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp @@ -79,7 +79,6 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase static constexpr char const * useRockDensityString() { return "useRockDensity"; } static constexpr char const * usePorosityString() { return "usePorosity"; } static constexpr char const * useReferencePorosityString() { return "useReferencePorosity"; } - constexpr static char const * flowSolverNameString() { return "flowSolverName"; } } gravityViewKeys; @@ -89,10 +88,6 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase virtual void initializePostInitialConditionsPreSubGroups() override final; - /// pointer to the flow solver - CompositionalMultiphaseFVM * m_flowSolver; - string m_flowSolverName; - /// Use rock density in addition to fluid density localIndex m_useRockDensity; localIndex m_useReferencePorosity; From 55833d4cecb760b95e83234911b979db0ed7aa11 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 18 Sep 2025 17:21:39 -0500 Subject: [PATCH 26/30] Used isLogLevelActive --- .../physicsSolvers/gravity/CMakeLists.txt | 1 + .../physicsSolvers/gravity/GravityFE.cpp | 7 ++++--- .../gravity/GravityFE_CompositionalMultiphaseFVM.cpp | 5 +++-- .../physicsSolvers/gravity/GravitySolverBase.cpp | 10 +++++----- 4 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt index c4663d1824d..8509a068963 100644 --- a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt @@ -18,6 +18,7 @@ Contains: # Specify solver headers set( gravitySolvers_headers + GravityLogLevelsInfo.hpp GravitySolverBase.hpp GravityFE.hpp GravityFEKernel.hpp ) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 9eb4fd1b200..5195bfa2491 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -20,6 +20,7 @@ #include "GravityFE.hpp" #include "GravityFEKernel.hpp" +#include "GravityLogLevelsInfo.hpp" #include "finiteElement/FiniteElementDiscretization.hpp" #include "mesh/DomainPartition.hpp" @@ -132,7 +133,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, CellElementSubRegion & elementSubRegion ) { arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()); - if( this->getLogLevel()>3 ) + if( geos::isLogLevelActive< geos::gravity::GravityPropertiesDebug >( getLogLevel() ) ) { for( localIndex i=0; igetLogLevel()>1 ) + if( geos::isLogLevelActive< geos::gravity::GravityComponentDebug >( getLogLevel() ) ) { for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) { @@ -341,7 +342,7 @@ real64 GravityFE::explicitStepAdjoint( real64 const & time_n, adjoint.move( LvArray::MemorySpace::host, true ); - if( this->getLogLevel() > 2 ) + if( geos::isLogLevelActive< geos::gravity::GravityAdjointDebug >( getLogLevel() ) ) { for( localIndex i=0; igetLogLevel() > 3 ) + if( geos::isLogLevelActive< geos::gravity::GravityPropertiesDebug >( getLogLevel() ) ) { // Move arrays to host for logging density.move( LvArray::MemorySpace::host, true ); @@ -344,7 +345,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const gzAtStations[i] = localGzAtStations[i]; } - if( this->getLogLevel() > 1 ) + if( geos::isLogLevelActive< geos::gravity::GravityComponentDebug >( getLogLevel() ) ) { for( localIndex iStation = 0; iStation < m_stationCoordinates.size( 0 ); ++iStation ) { diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp index 71ef4808517..69bd27f560f 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp @@ -18,6 +18,7 @@ */ #include "GravitySolverBase.hpp" +#include "GravityLogLevelsInfo.hpp" #include @@ -115,8 +116,7 @@ void GravitySolverBase::initializePostInitialConditionsPreSubGroups() ); } - - if( this->getLogLevel()>0 ) + if( geos::isLogLevelActive< geos::gravity::GravitySolverStatus >( getLogLevel() ) ) { constexpr auto yesno = [] (int flag) noexcept->const char * { return flag ? "yes" : "no"; }; @@ -126,9 +126,9 @@ void GravitySolverBase::initializePostInitialConditionsPreSubGroups() GEOS_LOG_RANK_0( "Mode: " << m_modeString ); GEOS_LOG_RANK_0( "Number of stations: " << m_stationCoordinates.size( 0 )); GEOS_LOG_RANK_0( "Log level: " << getLogLevel()); - GEOS_LOG_RANK_0( " Output Gz to logs: " << yesno( getLogLevel() > 1 )); - GEOS_LOG_RANK_0( " Output Adjoint to logs: " << yesno( getLogLevel() > 2 )); - GEOS_LOG_RANK_0( " Output Properties to logs: " << yesno( getLogLevel() > 3 )); + GEOS_LOG_RANK_0( " Output Gz to logs: " << yesno( geos::isLogLevelActive< geos::gravity::GravityComponentDebug >( getLogLevel() ) ) ); + GEOS_LOG_RANK_0( " Output Adjoint to logs: " << yesno( geos::isLogLevelActive< geos::gravity::GravityAdjointDebug >( getLogLevel() ) ) ); + GEOS_LOG_RANK_0( " Output Properties to logs: " << yesno( geos::isLogLevelActive< geos::gravity::GravityPropertiesDebug >( getLogLevel() ) ) ); gravitySolverLog.end(); } } From a5104d1d35c457a4fa4820d265a32e1a4b1c3ea4 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 18 Sep 2025 18:15:59 -0500 Subject: [PATCH 27/30] Seperated fields --- .../gravityTests/testGravitySolverAdjoint.cpp | 5 +- .../physicsSolvers/gravity/CMakeLists.txt | 1 + .../physicsSolvers/gravity/GravityFE.cpp | 31 +++---- .../physicsSolvers/gravity/GravityFE.hpp | 44 ---------- .../GravityFE_CompositionalMultiphaseFVM.cpp | 43 +++------ .../GravityFE_CompositionalMultiphaseFVM.hpp | 49 +---------- .../physicsSolvers/gravity/GravityFields.hpp | 88 +++++++++++++++++++ .../gravity/GravityLogLevelsInfo.hpp | 75 ++++++++++++++++ 8 files changed, 189 insertions(+), 147 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/gravity/GravityFields.hpp create mode 100644 src/coreComponents/physicsSolvers/gravity/GravityLogLevelsInfo.hpp diff --git a/src/coreComponents/integrationTests/gravityTests/testGravitySolverAdjoint.cpp b/src/coreComponents/integrationTests/gravityTests/testGravitySolverAdjoint.cpp index f5fbd8e47ba..cca63de83d9 100644 --- a/src/coreComponents/integrationTests/gravityTests/testGravitySolverAdjoint.cpp +++ b/src/coreComponents/integrationTests/gravityTests/testGravitySolverAdjoint.cpp @@ -5,6 +5,7 @@ #include "mainInterface/ProblemManager.hpp" #include "physicsSolvers/PhysicsSolverManager.hpp" #include "physicsSolvers/gravity/GravityFE.hpp" +#include "physicsSolvers/gravity/GravityFields.hpp" #include "mesh/DomainPartition.hpp" #include "integrationTests/fluidFlowTests/testCompFlowUtils.hpp" // For setupProblemFromXML #include "physicsSolvers/PhysicsSolverBase.hpp" @@ -129,7 +130,7 @@ TEST( GravityFEKernelTest, AdjointConsistencyWithGravityFE ) mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & elementSubRegion ) { - arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()); + arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::mediumDensity::key()); for( localIndex i = 0; i < density.size(); ++i ) { const_cast< real64 & >(density[i]) = xVec[i]; @@ -167,7 +168,7 @@ TEST( GravityFEKernelTest, AdjointConsistencyWithGravityFE ) mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & elementSubRegion ) { - arrayView1d< real64 > const adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()); + arrayView1d< real64 > const adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::adjoint::key()); for( localIndex i = 0; i < adjoint.size(); ++i ) { adjointValues.push_back( adjoint[i] ); diff --git a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt index 8509a068963..6b869b11cad 100644 --- a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt @@ -19,6 +19,7 @@ Contains: # Specify solver headers set( gravitySolvers_headers GravityLogLevelsInfo.hpp + GravityFields.hpp GravitySolverBase.hpp GravityFE.hpp GravityFEKernel.hpp ) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index 5195bfa2491..f4040f87a7c 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -19,6 +19,7 @@ */ #include "GravityFE.hpp" +#include "GravityFields.hpp" #include "GravityFEKernel.hpp" #include "GravityLogLevelsInfo.hpp" @@ -73,34 +74,22 @@ void GravityFE::registerDataOnMesh( Group & meshBodies ) string_array const & ) { NodeManager & nodeManager = mesh.getNodeManager(); - nodeManager.registerField< fields::VolumeIntegral >( this->getName() ); + nodeManager.registerField< fields::volumeIntegral >( this->getName() ); ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion ) { - subRegion.registerField< fields::MediumDensity >( this->getName() ); - subRegion.registerField< fields::Adjoint >( this->getName() ); - subRegion.registerField< fields::VolumeIntegral2d >( this->getName() ); + subRegion.registerField< fields::mediumDensity >( this->getName() ); + subRegion.registerField< fields::adjoint >( this->getName() ); + subRegion.registerField< fields::volumeIntegral2d >( this->getName() ); // Assume the maximum number of points per cell is MAX_SUPPORT_POINTS... - subRegion.getField< fields::VolumeIntegral2d >().resizeDimension< 1 >( MAX_SUPPORT_POINTS ); + subRegion.getField< fields::volumeIntegral2d >().resizeDimension< 1 >( MAX_SUPPORT_POINTS ); } ); } ); } -void GravityFE::postInputInitialization() -{ - GravitySolverBase::postInputInitialization(); -} - - -void GravityFE::initializePostInitialConditionsPreSubGroups() -{ - GravitySolverBase::initializePostInitialConditionsPreSubGroups(); -} - - real64 GravityFE::explicitStepModeling( real64 const & time_n, real64 const & dt, integer const cycleNumber, @@ -123,7 +112,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, auto const nodeGhostRank = nodeManager.ghostRank().toViewConst(); // VolumeIntegral matrix to be computed in this function. - auto volumeIntegral = nodeManager.getField< fields::VolumeIntegral >().toView(); + auto volumeIntegral = nodeManager.getField< fields::volumeIntegral >().toView(); volumeIntegral.setValues< parallelHostPolicy >( 0. ); @@ -132,7 +121,7 @@ real64 GravityFE::explicitStepModeling( real64 const & time_n, mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, CellElementSubRegion & elementSubRegion ) { - arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()); + arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::mediumDensity::key()); if( geos::isLogLevelActive< geos::gravity::GravityPropertiesDebug >( getLogLevel() ) ) { for( localIndex i=0; i const X = nodeManager.referencePosition().toViewConst(); // Output arrays - arrayView2d< real64 > volumeIntegral2d = elementSubRegion.getReference< array2d< real64 > >( fields::VolumeIntegral2d::key()).toView(); - arrayView1d< real64 > adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::Adjoint::key()).toView(); + arrayView2d< real64 > volumeIntegral2d = elementSubRegion.getReference< array2d< real64 > >( fields::volumeIntegral2d::key()).toView(); + arrayView1d< real64 > adjoint = elementSubRegion.getReference< array1d< real64 > >( fields::adjoint::key()).toView(); volumeIntegral2d.setValues< parallelHostPolicy >( 0.0 ); adjoint.zero(); diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp index 76942191586..66160a8e78c 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.hpp @@ -78,52 +78,8 @@ class GravityFE : public GravitySolverBase DomainPartition & domain ) override; /**@}*/ - - -protected: - - virtual void postInputInitialization() override final; - - virtual void initializePostInitialConditionsPreSubGroups() override final; - }; - -namespace fields -{ -DECLARE_FIELD( MediumDensity, - "mediumDensity", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "Medium density of the cell" ); - -DECLARE_FIELD( VolumeIntegral, - "volumeIntegral", - array1d< real64 >, - 0, - NOPLOT, - WRITE_AND_READ, - "VolumeIntegral of the cell." ); - -DECLARE_FIELD( VolumeIntegral2d, - "volumeIntegral2d", - array2d< real64 >, - 0, - NOPLOT, - WRITE_AND_READ, - "VolumeIntegral for adjoint computation." ); - -DECLARE_FIELD( Adjoint, - "adjoint", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "Adjoint field." ); -} - } // namespace geos #endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_HPP_ diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp index 7f97e5ccc88..ca847df0e4e 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp @@ -20,6 +20,7 @@ #include #include "GravityFE_CompositionalMultiphaseFVM.hpp" +#include "GravityFields.hpp" #include "GravityFEKernel.hpp" #include "GravityLogLevelsInfo.hpp" @@ -52,8 +53,6 @@ GravityFE_CompositionalMultiphaseFVM::GravityFE_CompositionalMultiphaseFVM( cons Group * const parent ): GravitySolverBase( name, parent ) { - - registerWrapper( viewKeyStruct::useRockDensityString(), &m_useRockDensity ). setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 1 ). @@ -70,11 +69,6 @@ GravityFE_CompositionalMultiphaseFVM::GravityFE_CompositionalMultiphaseFVM( cons setDescription( "Flag to use the reference porosity" ); } -GravityFE_CompositionalMultiphaseFVM::~GravityFE_CompositionalMultiphaseFVM() -{ - // TODO Auto-generated destructor stub -} - void GravityFE_CompositionalMultiphaseFVM::initializePreSubGroups() { @@ -105,36 +99,21 @@ void GravityFE_CompositionalMultiphaseFVM::registerDataOnMesh( Group & meshBodie { NodeManager & nodeManager = mesh.getNodeManager(); - nodeManager.registerField< fields::VolumeIntegral >( this->getName() ); + nodeManager.registerField< fields::volumeIntegral >( this->getName() ); ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion ) { - subRegion.registerField< fields::MediumDensity >( this->getName() ); - subRegion.registerField< fields::FluidDensity >( this->getName() ); - subRegion.registerField< fields::RockDensity >( this->getName() ); - subRegion.registerField< fields::Porosity >( this->getName() ); + subRegion.registerField< fields::mediumDensity >( this->getName() ); + subRegion.registerField< fields::fluidDensity >( this->getName() ); + subRegion.registerField< fields::rockDensity >( this->getName() ); + subRegion.registerField< fields::porosity >( this->getName() ); } ); } ); } -void GravityFE_CompositionalMultiphaseFVM::postInputInitialization() -{ - GravitySolverBase::postInputInitialization(); - -} - - - -void GravityFE_CompositionalMultiphaseFVM::initializePostInitialConditionsPreSubGroups() -{ - GravitySolverBase::initializePostInitialConditionsPreSubGroups(); -} - - - real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const & time_n, real64 const & dt, integer const cycleNumber, @@ -157,7 +136,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const auto const nodeGhostRank = nodeManager.ghostRank().toViewConst(); // VolumeIntegral matrix to be computed in this function. - auto volumeIntegral = nodeManager.getField< fields::VolumeIntegral >().toView(); + auto volumeIntegral = nodeManager.getField< fields::volumeIntegral >().toView(); volumeIntegral.setValues< parallelHostPolicy >( 0. ); // Loop over all sub-regions in regions of type "CellElements". @@ -171,7 +150,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const MultiFluidBase const & fluid = constitutiveModels.getGroup< MultiFluidBase >( fluidName ); arrayView2d< real64 const, multifluid::USD_FLUID > const reservoirFluidDensity = fluid.totalDensity().toViewConst(); - arrayView1d< real64 > const fluidDensity = elementSubRegion.getReference< array1d< real64 > >( fields::FluidDensity::key()).toView(); + arrayView1d< real64 > const fluidDensity = elementSubRegion.getReference< array1d< real64 > >( fields::fluidDensity::key()).toView(); // Fluid density assignment forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) @@ -182,7 +161,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const // Retrieve porosity. string const & solidName = elementSubRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::solidNamesString() ); CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); - arrayView1d< real64 > const porosity = elementSubRegion.getReference< array1d< real64 > >( fields::Porosity::key()).toView(); + arrayView1d< real64 > const porosity = elementSubRegion.getReference< array1d< real64 > >( fields::porosity::key()).toView(); if( this->m_useReferencePorosity == 1 ) { @@ -215,7 +194,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const } // Retrieve rock density. - arrayView1d< real64 > const rockDensity = elementSubRegion.getReference< array1d< real64 > >( fields::RockDensity::key()).toView(); + arrayView1d< real64 > const rockDensity = elementSubRegion.getReference< array1d< real64 > >( fields::rockDensity::key()).toView(); rockDensity.setValues< parallelHostPolicy >( 0. ); if( this->m_useRockDensity == 1 ) @@ -231,7 +210,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const } // Build medium density. - arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::MediumDensity::key()).toView(); + arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::mediumDensity::key()).toView(); // Make input arrays const for density calculation arrayView1d< real64 const > const porosityConst = porosity.toViewConst(); diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp index 4ce0d908803..8d74c09e06f 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp @@ -37,7 +37,6 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase GravityFE_CompositionalMultiphaseFVM() = delete; GravityFE_CompositionalMultiphaseFVM( const std::string & name, Group * const parent ); - virtual ~GravityFE_CompositionalMultiphaseFVM() override; GravityFE_CompositionalMultiphaseFVM( GravityFE_CompositionalMultiphaseFVM const & ) = delete; GravityFE_CompositionalMultiphaseFVM( GravityFE_CompositionalMultiphaseFVM && ) = delete; @@ -47,9 +46,8 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase static string catalogName() { return "GravityFE_CompositionalMultiphaseFVM"; } string getCatalogName() const override { return catalogName(); } - virtual void initializePreSubGroups() override; virtual void registerDataOnMesh( Group & meshBodies ) override final; - virtual void postInputInitialization() override final; + virtual void initializePreSubGroups() override; /** @@ -84,10 +82,6 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase protected: - - - virtual void initializePostInitialConditionsPreSubGroups() override final; - /// Use rock density in addition to fluid density localIndex m_useRockDensity; localIndex m_useReferencePorosity; @@ -99,47 +93,6 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase }; -namespace fields -{ - -DECLARE_FIELD( MediumDensity, - "mediumDensity", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "Medium density of the cell" ); -DECLARE_FIELD( FluidDensity, - "fluidDensity", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "Fluid density of the cell" ); -DECLARE_FIELD( RockDensity, - "rockDensity", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "Rock density of the cell" ); -DECLARE_FIELD( Porosity, - "porosity", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "Porosity of the cell" ); -DECLARE_FIELD( VolumeIntegral, - "volumeIntegral", - array1d< real64 >, - 0, - NOPLOT, - WRITE_AND_READ, - "VolumeIntegral of the cell." ); - -} - } // namespace geos #endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFields.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFields.hpp new file mode 100644 index 00000000000..e2804a8d87f --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravityFields.hpp @@ -0,0 +1,88 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * (Copy the appropriate license block from one of the other files) + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFIELDS_HPP_ +#define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFIELDS_HPP_ + +#include "mesh/MeshFields.hpp" // Required for the DECLARE_FIELD macro + +namespace geos +{ +namespace fields +{ + +//================================================================================ +// Fields shared by multiple gravity solvers +//================================================================================ + +DECLARE_FIELD( mediumDensity, + "mediumDensity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Medium density of the cell" ); + +DECLARE_FIELD( volumeIntegral, + "volumeIntegral", + array1d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "VolumeIntegral of the cell." ); + +//================================================================================ +// Fields specific to the GravityFE solver +//================================================================================ + +DECLARE_FIELD( volumeIntegral2d, + "volumeIntegral2d", + array2d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "VolumeIntegral for adjoint computation." ); + +DECLARE_FIELD( adjoint, + "adjoint", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Adjoint field." ); + +//================================================================================ +// Fields specific to the GravityFE_CompositionalMultiphaseFVM solver +//================================================================================ + +DECLARE_FIELD( fluidDensity, + "fluidDensity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Fluid density of the cell" ); + +DECLARE_FIELD( rockDensity, + "rockDensity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Rock density of the cell" ); + +DECLARE_FIELD( porosity, + "porosity", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Porosity of the cell" ); + +} // namespace fields +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFIELDS_HPP_ diff --git a/src/coreComponents/physicsSolvers/gravity/GravityLogLevelsInfo.hpp b/src/coreComponents/physicsSolvers/gravity/GravityLogLevelsInfo.hpp new file mode 100644 index 00000000000..6961a1586bf --- /dev/null +++ b/src/coreComponents/physicsSolvers/gravity/GravityLogLevelsInfo.hpp @@ -0,0 +1,75 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 GravityLogLevels.hpp + */ + + + +#ifndef GEOS_SOLVERS_GRAVITY_GRAVITYLOGLEVELSINFO_HPP +#define GEOS_SOLVERS_GRAVITY_GRAVITYLOGLEVELSINFO_HPP + +#include "common/DataTypes.hpp" + +namespace geos +{ +namespace gravity +{ +/// Log basic solver status and configuration. +struct GravitySolverStatus +{ + static constexpr integer getMinLogLevel() { return 1; } + static constexpr std::string_view getDescription() + { + return "Log basic Gravity solver status and configuration on initialization."; + } +}; + +/// Log the computed vertical gravity component (Gz). +struct GravityComponentDebug +{ + static constexpr integer getMinLogLevel() { return 2; } + static constexpr std::string_view getDescription() + { + return "Log the calculated vertical gravity component (Gz) for debugging."; + } +}; + +/// Log adjoint-related computations. +struct GravityAdjointDebug +{ + static constexpr integer getMinLogLevel() { return 3; } + static constexpr std::string_view getDescription() + { + return "Log the adjoint calculations for debugging."; + } +}; + +/// Log detailed, per-element properties like density. +struct GravityPropertiesDebug +{ + static constexpr integer getMinLogLevel() { return 4; } + static constexpr std::string_view getDescription() + { + return "Log detailed per-element properties (e.g., density)."; + } +}; + +} // namespace gravity +} // namespace geos + +#endif // GEOS_SOLVERS_GRAVITY_GRAVITYLOGLEVELSINFO_HPP From 77feca234bb951b1a0fd46c0dc4b4d362ede6294 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 18 Sep 2025 19:11:27 -0500 Subject: [PATCH 28/30] Renamed class --- .../physicsSolvers/gravity/CMakeLists.txt | 4 +-- ...GravityFE_CompositionalMultiphaseBase.cpp} | 28 +++++++++---------- ...GravityFE_CompositionalMultiphaseBase.hpp} | 25 ++++++++--------- 3 files changed, 28 insertions(+), 29 deletions(-) rename src/coreComponents/physicsSolvers/gravity/{GravityFE_CompositionalMultiphaseFVM.cpp => GravityFE_CompositionalMultiphaseBase.cpp} (92%) rename src/coreComponents/physicsSolvers/gravity/{GravityFE_CompositionalMultiphaseFVM.hpp => GravityFE_CompositionalMultiphaseBase.hpp} (76%) diff --git a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt index 6b869b11cad..76bf059865e 100644 --- a/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt @@ -33,8 +33,8 @@ set( dependencyList ${parallelDeps} physicsSolversBase ) # Conditionally add compositional multiphase support if fluid flow is enabled if( GEOS_ENABLE_FLUIDFLOW AND TARGET fluidFlowSolvers ) - list( APPEND gravitySolvers_headers GravityFE_CompositionalMultiphaseFVM.hpp ) - list( APPEND gravitySolvers_sources GravityFE_CompositionalMultiphaseFVM.cpp ) + list( APPEND gravitySolvers_headers GravityFE_CompositionalMultiphaseBase.hpp ) + list( APPEND gravitySolvers_sources GravityFE_CompositionalMultiphaseBase.cpp ) list( APPEND dependencyList fluidFlowSolvers ) # Define a preprocessor macro so the code knows fluid flow is available diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp similarity index 92% rename from src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp rename to src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp index ca847df0e4e..e817afe14b2 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp @@ -15,11 +15,11 @@ /** - * @file GravityFE_CompositionalMultiphaseFVM.cpp + * @file GravityFE_CompositionalMultiphaseBase.cpp */ #include -#include "GravityFE_CompositionalMultiphaseFVM.hpp" +#include "GravityFE_CompositionalMultiphaseBase.hpp" #include "GravityFields.hpp" #include "GravityFEKernel.hpp" #include "GravityLogLevelsInfo.hpp" @@ -35,7 +35,7 @@ #include "mesh/ElementType.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/solid/CoupledSolid.hpp" #include "constitutive/solid/SolidBase.hpp" @@ -49,7 +49,7 @@ using namespace constitutive; using namespace dataRepository; -GravityFE_CompositionalMultiphaseFVM::GravityFE_CompositionalMultiphaseFVM( const std::string & name, +GravityFE_CompositionalMultiphaseBase::GravityFE_CompositionalMultiphaseBase( const std::string & name, Group * const parent ): GravitySolverBase( name, parent ) { @@ -70,7 +70,7 @@ GravityFE_CompositionalMultiphaseFVM::GravityFE_CompositionalMultiphaseFVM( cons } -void GravityFE_CompositionalMultiphaseFVM::initializePreSubGroups() +void GravityFE_CompositionalMultiphaseBase::initializePreSubGroups() { GravitySolverBase::initializePreSubGroups(); DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); @@ -90,7 +90,7 @@ void GravityFE_CompositionalMultiphaseFVM::initializePreSubGroups() } -void GravityFE_CompositionalMultiphaseFVM::registerDataOnMesh( Group & meshBodies ) +void GravityFE_CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) { forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, @@ -114,7 +114,7 @@ void GravityFE_CompositionalMultiphaseFVM::registerDataOnMesh( Group & meshBodie } -real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const & time_n, +real64 GravityFE_CompositionalMultiphaseBase::explicitStepModeling( real64 const & time_n, real64 const & dt, integer const cycleNumber, DomainPartition & domain ) @@ -165,7 +165,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const if( this->m_useReferencePorosity == 1 ) { - GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Use referencePorosity" ); + GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseBase: Use referencePorosity" ); arrayView1d< real64 const > const reservoirPorosity = solid.getReferencePorosity().toViewConst(); forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) @@ -175,7 +175,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const } else if( this->m_usePorosity == 1 ) { - GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Use Porosity" ); + GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseBase: Use Porosity" ); arrayView2d< real64 const > const reservoirPorosity = solid.getPorosity().toViewConst(); forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) @@ -185,7 +185,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const } else { - GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseFVM: Set Porosity to 1" ); + GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseBase: Set Porosity to 1" ); forAll< EXEC_POLICY >( elementSubRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const i ) { @@ -233,7 +233,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const for( localIndex i = 0; i < elementSubRegion.size(); ++i ) { - GEOS_LOG( "GravityFE_CompositionalMultiphaseFVM: Cell[" << i << "], density= " << density[i] + GEOS_LOG( "GravityFE_CompositionalMultiphaseBase: Cell[" << i << "], density= " << density[i] << ", fluidDensity= " << fluidDensityConst[i] << ", rockDensity= " << rockDensityConst[i] << ", porosity= " << porosityConst[i] ); } @@ -331,7 +331,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const auto const & coords = m_stationCoordinates[iStation]; std::ostringstream logStream; logStream << std::fixed << std::setprecision( 2 ); - logStream << "GravityFE_CompositionalMultiphaseFVM: station[" << std::setw( 5 ) << iStation << "] " + logStream << "GravityFE_CompositionalMultiphaseBase: station[" << std::setw( 5 ) << iStation << "] " << std::setw( 15 ) << coords[0] << " " << std::setw( 15 ) << coords[1] << " " << std::setw( 10 ) << coords[2] << " " @@ -345,7 +345,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepModeling( real64 const } -real64 GravityFE_CompositionalMultiphaseFVM::explicitStepAdjoint( real64 const & time_n, +real64 GravityFE_CompositionalMultiphaseBase::explicitStepAdjoint( real64 const & time_n, real64 const & dt, integer const cycleNumber, DomainPartition & domain ) @@ -357,7 +357,7 @@ real64 GravityFE_CompositionalMultiphaseFVM::explicitStepAdjoint( real64 const & #ifdef GEOS_GRAVITY_WITH_FLUIDFLOW -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE_CompositionalMultiphaseFVM, string const &, dataRepository::Group * const ) +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE_CompositionalMultiphaseBase, string const &, dataRepository::Group * const ) #endif } // namespace geos diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp similarity index 76% rename from src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp rename to src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp index 8d74c09e06f..b27d687248b 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp @@ -15,35 +15,34 @@ /** - * @file GravityFE_CompositionalMultiphaseFVM.hpp + * @file GravityFE_CompositionalMultiphaseBase.hpp */ -#ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ -#define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ +#ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEBASE_HPP_ +#define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEBASE_HPP_ #include "GravitySolverBase.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" namespace geos { -class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase +class GravityFE_CompositionalMultiphaseBase : public GravitySolverBase { public: using EXEC_POLICY = parallelDevicePolicy< 32 >; using ATOMIC_POLICY = parallelDeviceAtomic; - GravityFE_CompositionalMultiphaseFVM() = delete; - GravityFE_CompositionalMultiphaseFVM( const std::string & name, + GravityFE_CompositionalMultiphaseBase() = delete; + GravityFE_CompositionalMultiphaseBase( const std::string & name, Group * const parent ); - GravityFE_CompositionalMultiphaseFVM( GravityFE_CompositionalMultiphaseFVM const & ) = delete; - GravityFE_CompositionalMultiphaseFVM( GravityFE_CompositionalMultiphaseFVM && ) = delete; - GravityFE_CompositionalMultiphaseFVM & operator=( GravityFE_CompositionalMultiphaseFVM const & ) = delete; - GravityFE_CompositionalMultiphaseFVM & operator=( GravityFE_CompositionalMultiphaseFVM && ) = delete; + GravityFE_CompositionalMultiphaseBase( GravityFE_CompositionalMultiphaseBase const & ) = delete; + GravityFE_CompositionalMultiphaseBase( GravityFE_CompositionalMultiphaseBase && ) = delete; + GravityFE_CompositionalMultiphaseBase & operator=( GravityFE_CompositionalMultiphaseBase const & ) = delete; + GravityFE_CompositionalMultiphaseBase & operator=( GravityFE_CompositionalMultiphaseBase && ) = delete; - static string catalogName() { return "GravityFE_CompositionalMultiphaseFVM"; } + static string catalogName() { return "GravityFE_CompositionalMultiphaseBase"; } string getCatalogName() const override { return catalogName(); } virtual void registerDataOnMesh( Group & meshBodies ) override final; @@ -95,4 +94,4 @@ class GravityFE_CompositionalMultiphaseFVM : public GravitySolverBase } // namespace geos -#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEFVM_HPP_ +#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEBASE_HPP_ From 03e049fd4d25a28c3143b41ae3a8ad3cfe9b101b Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 18 Sep 2025 19:19:00 -0500 Subject: [PATCH 29/30] Removed this. --- .../physicsSolvers/gravity/GravityFE.cpp | 8 ++++---- .../GravityFE_CompositionalMultiphaseBase.cpp | 16 ++++++++-------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp index f4040f87a7c..fe38bfffa69 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp @@ -74,14 +74,14 @@ void GravityFE::registerDataOnMesh( Group & meshBodies ) string_array const & ) { NodeManager & nodeManager = mesh.getNodeManager(); - nodeManager.registerField< fields::volumeIntegral >( this->getName() ); + nodeManager.registerField< fields::volumeIntegral >( getName() ); ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion ) { - subRegion.registerField< fields::mediumDensity >( this->getName() ); - subRegion.registerField< fields::adjoint >( this->getName() ); - subRegion.registerField< fields::volumeIntegral2d >( this->getName() ); + subRegion.registerField< fields::mediumDensity >( getName() ); + subRegion.registerField< fields::adjoint >( getName() ); + subRegion.registerField< fields::volumeIntegral2d >( getName() ); // Assume the maximum number of points per cell is MAX_SUPPORT_POINTS... subRegion.getField< fields::volumeIntegral2d >().resizeDimension< 1 >( MAX_SUPPORT_POINTS ); diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp index e817afe14b2..8cc55046da7 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp @@ -99,16 +99,16 @@ void GravityFE_CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodi { NodeManager & nodeManager = mesh.getNodeManager(); - nodeManager.registerField< fields::volumeIntegral >( this->getName() ); + nodeManager.registerField< fields::volumeIntegral >( getName() ); ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion ) { - subRegion.registerField< fields::mediumDensity >( this->getName() ); - subRegion.registerField< fields::fluidDensity >( this->getName() ); - subRegion.registerField< fields::rockDensity >( this->getName() ); - subRegion.registerField< fields::porosity >( this->getName() ); + subRegion.registerField< fields::mediumDensity >( getName() ); + subRegion.registerField< fields::fluidDensity >( getName() ); + subRegion.registerField< fields::rockDensity >( getName() ); + subRegion.registerField< fields::porosity >( getName() ); } ); } ); } @@ -163,7 +163,7 @@ real64 GravityFE_CompositionalMultiphaseBase::explicitStepModeling( real64 const CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); arrayView1d< real64 > const porosity = elementSubRegion.getReference< array1d< real64 > >( fields::porosity::key()).toView(); - if( this->m_useReferencePorosity == 1 ) + if( m_useReferencePorosity == 1 ) { GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseBase: Use referencePorosity" ); arrayView1d< real64 const > const reservoirPorosity = solid.getReferencePorosity().toViewConst(); @@ -173,7 +173,7 @@ real64 GravityFE_CompositionalMultiphaseBase::explicitStepModeling( real64 const porosity[i] = reservoirPorosity[i]; } ); } - else if( this->m_usePorosity == 1 ) + else if( m_usePorosity == 1 ) { GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseBase: Use Porosity" ); arrayView2d< real64 const > const reservoirPorosity = solid.getPorosity().toViewConst(); @@ -197,7 +197,7 @@ real64 GravityFE_CompositionalMultiphaseBase::explicitStepModeling( real64 const arrayView1d< real64 > const rockDensity = elementSubRegion.getReference< array1d< real64 > >( fields::rockDensity::key()).toView(); rockDensity.setValues< parallelHostPolicy >( 0. ); - if( this->m_useRockDensity == 1 ) + if( m_useRockDensity == 1 ) { string const & solidModelName = elementSubRegion.getReference< string >( SolidMechanicsLagrangianFEM::viewKeyStruct::solidMaterialNamesString()); arrayView2d< real64 const > const reservoirRockDensity = From b5c0692e1d8712addc5cbdd8bea64437590dceb8 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 18 Sep 2025 19:19:38 -0500 Subject: [PATCH 30/30] Uncrustify --- .../GravityFE_CompositionalMultiphaseBase.cpp | 18 +++++++++--------- .../GravityFE_CompositionalMultiphaseBase.hpp | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp index 8cc55046da7..d20824e76fe 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp @@ -50,7 +50,7 @@ using namespace constitutive; using namespace dataRepository; GravityFE_CompositionalMultiphaseBase::GravityFE_CompositionalMultiphaseBase( const std::string & name, - Group * const parent ): + Group * const parent ): GravitySolverBase( name, parent ) { registerWrapper( viewKeyStruct::useRockDensityString(), &m_useRockDensity ). @@ -115,9 +115,9 @@ void GravityFE_CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodi real64 GravityFE_CompositionalMultiphaseBase::explicitStepModeling( real64 const & time_n, - real64 const & dt, - integer const cycleNumber, - DomainPartition & domain ) + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) { GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( time_n, cycleNumber ); @@ -234,8 +234,8 @@ real64 GravityFE_CompositionalMultiphaseBase::explicitStepModeling( real64 const for( localIndex i = 0; i < elementSubRegion.size(); ++i ) { GEOS_LOG( "GravityFE_CompositionalMultiphaseBase: Cell[" << i << "], density= " << density[i] - << ", fluidDensity= " << fluidDensityConst[i] << ", rockDensity= " << rockDensityConst[i] - << ", porosity= " << porosityConst[i] ); + << ", fluidDensity= " << fluidDensityConst[i] << ", rockDensity= " << rockDensityConst[i] + << ", porosity= " << porosityConst[i] ); } } @@ -346,9 +346,9 @@ real64 GravityFE_CompositionalMultiphaseBase::explicitStepModeling( real64 const real64 GravityFE_CompositionalMultiphaseBase::explicitStepAdjoint( real64 const & time_n, - real64 const & dt, - integer const cycleNumber, - DomainPartition & domain ) + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) { GEOS_UNUSED_VAR( time_n, cycleNumber, domain ); GEOS_ERROR( "Adjoint computation is not available when Gravity is coupled with Flow simulator." ); diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp index b27d687248b..f147cbefe82 100644 --- a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp @@ -35,7 +35,7 @@ class GravityFE_CompositionalMultiphaseBase : public GravitySolverBase GravityFE_CompositionalMultiphaseBase() = delete; GravityFE_CompositionalMultiphaseBase( const std::string & name, - Group * const parent ); + Group * const parent ); GravityFE_CompositionalMultiphaseBase( GravityFE_CompositionalMultiphaseBase const & ) = delete; GravityFE_CompositionalMultiphaseBase( GravityFE_CompositionalMultiphaseBase && ) = delete;