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)
diff --git a/src/cmake/GeosxOptions.cmake b/src/cmake/GeosxOptions.cmake
index fe3adabcc4b..e32c1a647e3 100644
--- a/src/cmake/GeosxOptions.cmake
+++ b/src/cmake/GeosxOptions.cmake
@@ -128,6 +128,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/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..d9fcf73fdcc 100644
--- a/src/coreComponents/dataRepository/python/PyWrapper.hpp
+++ b/src/coreComponents/dataRepository/python/PyWrapper.hpp
@@ -24,6 +24,26 @@ namespace geos
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;
+};
+
+
/**
*
*/
@@ -34,6 +54,26 @@ 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 );
+
+
} // namespace python
} // namespace geos
diff --git a/src/coreComponents/dataRepository/unitTests/CMakeLists.txt b/src/coreComponents/dataRepository/unitTests/CMakeLists.txt
index 578071e570d..8d3b9e8e981 100644
--- a/src/coreComponents/dataRepository/unitTests/CMakeLists.txt
+++ b/src/coreComponents/dataRepository/unitTests/CMakeLists.txt
@@ -9,6 +9,10 @@ set( dataRepository_tests
testObjectCatalog.cpp
testWrapperHelpers.cpp )
+if( ENABLE_PYGEOSX )
+ list( APPEND dataRepository_tests testPyWrapper_setValue.cpp )
+endif()
+
set( dependencyList mainInterface gtest ${parallelDeps} )
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..aef65a795ea
--- /dev/null
+++ b/src/coreComponents/dataRepository/unitTests/testPyWrapper_setValue.cpp
@@ -0,0 +1,131 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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"
+#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 );
+}
diff --git a/src/coreComponents/fileIO/CMakeLists.txt b/src/coreComponents/fileIO/CMakeLists.txt
index eedca99b65c..84c1f0f3218 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..85ea2760059
--- /dev/null
+++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.cpp
@@ -0,0 +1,234 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 );
+ }
+
+ // 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,
+ 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 ) );
+ }
+ // 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];
+ }
+ }
+ }
+}
+
+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..3a27b48ebec
--- /dev/null
+++ b/src/coreComponents/fileIO/timeHistory/ScatterDataHistoryCollection.hpp
@@ -0,0 +1,155 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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:
+
+ /**
+ * @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
+ virtual ~ScatterDataHistoryCollection() = default;
+
+ /// Non-copyable
+ ScatterDataHistoryCollection( ScatterDataHistoryCollection const & ) = delete;
+ ScatterDataHistoryCollection & operator=( ScatterDataHistoryCollection const & ) = 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
+ * @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 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_ */
diff --git a/src/coreComponents/fileIO/timeHistory/unitTests/CMakeLists.txt b/src/coreComponents/fileIO/timeHistory/unitTests/CMakeLists.txt
index 6f5aa42c175..d1da6926e6f 100644
--- a/src/coreComponents/fileIO/timeHistory/unitTests/CMakeLists.txt
+++ b/src/coreComponents/fileIO/timeHistory/unitTests/CMakeLists.txt
@@ -1,7 +1,9 @@
# Specify list of tests
set( gtest_geosx_tests
testHDFFile.cpp
- testHDFParallelFile.cpp )
+ testHDFParallelFile.cpp
+ testScatterDataProvider.cpp
+ testScatterDataHistoryCollection.cpp )
set( dependencyList mainInterface gtest ${parallelDeps} HDF5::HDF5 )
diff --git a/src/coreComponents/fileIO/timeHistory/unitTests/testScatterDataHistoryCollection.cpp b/src/coreComponents/fileIO/timeHistory/unitTests/testScatterDataHistoryCollection.cpp
new file mode 100644
index 00000000000..e6d63d1d316
--- /dev/null
+++ b/src/coreComponents/fileIO/timeHistory/unitTests/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;
+}
diff --git a/src/coreComponents/fileIO/timeHistory/unitTests/testScatterDataProvider.cpp b/src/coreComponents/fileIO/timeHistory/unitTests/testScatterDataProvider.cpp
new file mode 100644
index 00000000000..ce3536d1ae1
--- /dev/null
+++ b/src/coreComponents/fileIO/timeHistory/unitTests/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();
+}
diff --git a/src/coreComponents/integrationTests/CMakeLists.txt b/src/coreComponents/integrationTests/CMakeLists.txt
index 38686825966..07a0bca4ced 100644
--- a/src/coreComponents/integrationTests/CMakeLists.txt
+++ b/src/coreComponents/integrationTests/CMakeLists.txt
@@ -32,3 +32,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/integrationTests/gravityTests/CMakeLists.txt b/src/coreComponents/integrationTests/gravityTests/CMakeLists.txt
new file mode 100644
index 00000000000..629d5307b08
--- /dev/null
+++ b/src/coreComponents/integrationTests/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/integrationTests/gravityTests/testGravitySolver.cpp b/src/coreComponents/integrationTests/gravityTests/testGravitySolver.cpp
new file mode 100644
index 00000000000..8ce7c155e13
--- /dev/null
+++ b/src/coreComponents/integrationTests/gravityTests/testGravitySolver.cpp
@@ -0,0 +1,198 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 "integrationTests/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/integrationTests/gravityTests/testGravitySolverAdjoint.cpp b/src/coreComponents/integrationTests/gravityTests/testGravitySolverAdjoint.cpp
new file mode 100644
index 00000000000..cca63de83d9
--- /dev/null
+++ b/src/coreComponents/integrationTests/gravityTests/testGravitySolverAdjoint.cpp
@@ -0,0 +1,203 @@
+
+#include
+#include "mainInterface/GeosxState.hpp"
+#include "mainInterface/initialization.hpp"
+#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"
+#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;
+}
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..76bf059865e
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/gravity/CMakeLists.txt
@@ -0,0 +1,61 @@
+# 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
+ GravityLogLevelsInfo.hpp
+ GravityFields.hpp
+ GravitySolverBase.hpp
+ GravityFE.hpp
+ GravityFEKernel.hpp )
+
+# Specify solver sources
+set( gravitySolvers_sources
+ GravitySolverBase.cpp
+ GravityFE.cpp )
+
+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_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
+ add_definitions( -DGEOS_GRAVITY_WITH_FLUIDFLOW )
+endif()
+
+
+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..fe38bfffa69
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/gravity/GravityFE.cpp
@@ -0,0 +1,350 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 "GravityFE.hpp"
+#include "GravityFields.hpp"
+#include "GravityFEKernel.hpp"
+#include "GravityLogLevelsInfo.hpp"
+
+#include "finiteElement/FiniteElementDiscretization.hpp"
+#include "mesh/DomainPartition.hpp"
+#include "common/MpiWrapper.hpp"
+
+
+namespace geos
+{
+
+using namespace constitutive;
+using namespace dataRepository;
+
+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 );
+ GEOS_THROW_IF( feDiscretization->getOrder() >1,
+ getName() << ": FE discretization order should be 1, but we get: " << feDiscretization->getOrder(),
+ InputError );
+}
+
+
+void GravityFE::registerDataOnMesh( Group & meshBodies )
+{
+ forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
+ MeshLevel & mesh,
+ string_array const & )
+ {
+ NodeManager & nodeManager = mesh.getNodeManager();
+ nodeManager.registerField< fields::volumeIntegral >( getName() );
+
+ ElementRegionManager & elemManager = mesh.getElemManager();
+ elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion )
+ {
+ 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 );
+ } );
+ } );
+}
+
+
+real64 GravityFE::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 )
+ {
+ arrayView1d< real64 > const density = elementSubRegion.getReference< array1d< real64 > >( fields::mediumDensity::key());
+ if( geos::isLogLevelActive< geos::gravity::GravityPropertiesDebug >( getLogLevel() ) )
+ {
+ for( localIndex i=0; i 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();
+
+ arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition =
+ nodeManager.referencePosition().toViewConst();
+
+ 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();
+ arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
+ arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > 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();
+
+ volumeIntegral2d.setValues< parallelHostPolicy >( 0.0 );
+ adjoint.zero();
+
+ 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 );
+
+ gravityFEKernel::
+ AdjointVolumeIntegralKernel< FE_TYPE > kernel( finiteElement );
+ kernel.template launch< EXEC_POLICY, ATOMIC_POLICY >
+ ( elementSubRegion.size(),
+ X,
+ elemsToNodes,
+ volumeIntegral2d );
+ } );
+
+ // Make volumeIntegral2d const
+ arrayView2d< real64 const > const volumeIntegral2dConst = volumeIntegral2d.toViewConst();
+
+ // Constants
+ constexpr real64 GRAVITATIONAL_CONSTANT_LOCAL = GRAVITATIONAL_CONSTANT;
+ constexpr real64 EPSILON = 1e-15;
+
+ // 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 )
+ {
+ 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;
+ 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 ) );
+ }
+
+ adjoint[k] += GRAVITATIONAL_CONSTANT_LOCAL * volumeIntegral2dConst( k, iLoc ) * res * dz * inv_r3;
+ }
+ }
+ } ); // Loop elem
+
+ adjoint.move( LvArray::MemorySpace::host, true );
+
+ if( geos::isLogLevelActive< geos::gravity::GravityAdjointDebug >( getLogLevel() ) )
+ {
+ 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 && ) = delete;
+
+ 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;
+ /**@}*/
+
+};
+
+} // 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..3238f4e154e
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/gravity/GravityFEKernel.hpp
@@ -0,0 +1,163 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+
+namespace geos
+{
+
+/// Namespace to contain the gravity kernels.
+namespace gravityFEKernel
+{
+
+template< typename FE_TYPE >
+struct ForwardVolumeIntegralKernel
+{
+ explicit ForwardVolumeIntegralKernel( FE_TYPE const & finiteElement )
+ : m_finiteElement( finiteElement )
+ {}
+
+ /**
+ * @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
+ * @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 AdjointVolumeIntegralKernel
+{
+ explicit AdjointVolumeIntegralKernel( 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/GravityFE_CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp
new file mode 100644
index 00000000000..d20824e76fe
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.cpp
@@ -0,0 +1,363 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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_CompositionalMultiphaseBase.cpp
+ */
+#include
+
+#include "GravityFE_CompositionalMultiphaseBase.hpp"
+#include "GravityFields.hpp"
+#include "GravityFEKernel.hpp"
+#include "GravityLogLevelsInfo.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/CompositionalMultiphaseBase.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_CompositionalMultiphaseBase::GravityFE_CompositionalMultiphaseBase( const std::string & name,
+ Group * const parent ):
+ GravitySolverBase( name, parent )
+{
+ 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" );
+}
+
+
+void GravityFE_CompositionalMultiphaseBase::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_CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies )
+{
+
+ forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
+ MeshLevel & mesh,
+ const string_array & )
+ {
+ NodeManager & nodeManager = mesh.getNodeManager();
+
+ nodeManager.registerField< fields::volumeIntegral >( getName() );
+
+ ElementRegionManager & elemManager = mesh.getElemManager();
+
+ elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion )
+ {
+ subRegion.registerField< fields::mediumDensity >( getName() );
+ subRegion.registerField< fields::fluidDensity >( getName() );
+ subRegion.registerField< fields::rockDensity >( getName() );
+ subRegion.registerField< fields::porosity >( getName() );
+ } );
+ } );
+}
+
+
+real64 GravityFE_CompositionalMultiphaseBase::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( m_useReferencePorosity == 1 )
+ {
+ 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 )
+ {
+ porosity[i] = reservoirPorosity[i];
+ } );
+ }
+ else if( m_usePorosity == 1 )
+ {
+ 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 )
+ {
+ porosity[i] = reservoirPorosity[i][0];
+ } );
+ }
+ else
+ {
+ GEOS_LOG_RANK_0( "GravityFE_CompositionalMultiphaseBase: 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( 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( geos::isLogLevelActive< geos::gravity::GravityPropertiesDebug >( getLogLevel() ) )
+ {
+ // 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_CompositionalMultiphaseBase: 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( geos::isLogLevelActive< geos::gravity::GravityComponentDebug >( getLogLevel() ) )
+ {
+ 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_CompositionalMultiphaseBase: 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_CompositionalMultiphaseBase::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;
+}
+
+
+#ifdef GEOS_GRAVITY_WITH_FLUIDFLOW
+REGISTER_CATALOG_ENTRY( PhysicsSolverBase, GravityFE_CompositionalMultiphaseBase, string const &, dataRepository::Group * const )
+#endif
+
+} // namespace geos
diff --git a/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp
new file mode 100644
index 00000000000..f147cbefe82
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/gravity/GravityFE_CompositionalMultiphaseBase.hpp
@@ -0,0 +1,97 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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_CompositionalMultiphaseBase.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEBASE_HPP_
+#define GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEBASE_HPP_
+
+#include "GravitySolverBase.hpp"
+
+namespace geos
+{
+
+class GravityFE_CompositionalMultiphaseBase : public GravitySolverBase
+{
+
+public:
+ using EXEC_POLICY = parallelDevicePolicy< 32 >;
+ using ATOMIC_POLICY = parallelDeviceAtomic;
+
+ GravityFE_CompositionalMultiphaseBase() = delete;
+ GravityFE_CompositionalMultiphaseBase( const std::string & name,
+ Group * const parent );
+
+ 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_CompositionalMultiphaseBase"; }
+ string getCatalogName() const override { return catalogName(); }
+
+ virtual void registerDataOnMesh( Group & meshBodies ) override final;
+ virtual void initializePreSubGroups() override;
+
+
+ /**
+ * @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"; }
+ } gravityViewKeys;
+
+
+protected:
+
+ /// Use rock density in addition to fluid density
+ localIndex m_useRockDensity;
+ localIndex m_useReferencePorosity;
+ localIndex m_usePorosity;
+
+private:
+
+
+
+};
+
+} // namespace geos
+
+#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYFE_COMPOSITIONALMULTUPHASEBASE_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
diff --git a/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp
new file mode 100644
index 00000000000..69bd27f560f
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.cpp
@@ -0,0 +1,202 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 "GravityLogLevelsInfo.hpp"
+#include
+
+
+namespace geos
+{
+
+using namespace dataRepository;
+
+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::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;
+
+
+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();
+}
+
+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 )
+ {
+
+
+ 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
+ );
+ }
+
+ if( geos::isLogLevelActive< geos::gravity::GravitySolverStatus >( getLogLevel() ) )
+ {
+ 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( "Log level: " << getLogLevel());
+ 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();
+ }
+}
+
+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 );
+ }
+}
+
+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
new file mode 100644
index 00000000000..dfd206888b5
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/gravity/GravitySolverBase.hpp
@@ -0,0 +1,127 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 "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, public ScatterDataProvider
+{
+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 && ) = delete;
+
+ 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 * residueString() { return "residue"; }
+ static constexpr char const * gzAtStationsString() { return "gzAtStations"; }
+ };
+
+ 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 };
+
+ 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,
+ integer const cycleNumber,
+ DomainPartition & domain ) = 0;
+
+ virtual real64 explicitStepAdjoint( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain ) = 0;
+
+ /// Coordinates of the gravimeter stations
+ array2d< real64 > m_stationCoordinates;
+
+ /// Residue observed at station (only for adjoint computation)
+ array1d< real64 > m_residue;
+
+ /// Gz component recorded at stations
+ array1d< real64 > m_gzAtStations;
+
+ /// Station metadata (names/identifiers)
+ mutable string_array m_stationMetadata;
+};
+
+} // namespace geos
+
+#endif // GEOS_PHYSICSSOLVERS_GRAVITY_GRAVITYSOLVERBASE_HPP_