diff --git a/.github/workflows/cmake-test.yml b/.github/workflows/cmake-test.yml index d736a952..e18de25e 100644 --- a/.github/workflows/cmake-test.yml +++ b/.github/workflows/cmake-test.yml @@ -67,7 +67,7 @@ jobs: if: matrix.python_api == 'ON' run: | sudo apt-get install -yq python3 python3-pip python3-dev python3-pybind11 pybind11-dev - pip3 install numpy + pip3 install numpy pytest - uses: actions/checkout@v4 @@ -288,12 +288,8 @@ jobs: run: | export PYTHONPATH=${{ runner.temp }}/build-pcms/install/lib/python3.12/site-packages:$PYTHONPATH export PYTHONPATH=${{ runner.temp }}/build-omega_h/install/lib/python/dist-packages:$PYTHONPATH - cd ${{ github.workspace }}/src/pcms/pythonapi - python3 test_file_io.py - python3 test_omega_h_field.py - python3 test_uniform_grid_field.py - python3 test_field_copy.py - python3 test_mls_interpolation.py + cd ${{ github.workspace }} + pytest -v --tb=short -p no:cacheprovider - name: Test PCMS Installation if: matrix.python_api == 'OFF' diff --git a/.github/workflows/self-hosted.yml b/.github/workflows/self-hosted.yml index c4dc85b8..ea5f2757 100644 --- a/.github/workflows/self-hosted.yml +++ b/.github/workflows/self-hosted.yml @@ -48,6 +48,7 @@ jobs: #module load openblas/0.3.23-wqm7iud module load netlib-lapack/3.11.0-b22mgwg #netlib-lapack includes blas module load netlib-scalapack/2.2.0-fzd4jvl + module load python py-pybind11 py-pytest py-numpy set -e EOF @@ -98,13 +99,13 @@ jobs: cmake --build $bdir -j 4 ctest --test-dir $bdir --output-on-failure - # Build pcms with PETSc + # Build pcms with PETSc and python enabled cmake -S ${{github.workspace}}/pcms_${{ github.event.id }} -B $bdir \ -DCMAKE_CXX_COMPILER=`which mpicxx` \ -DCMAKE_C_COMPILER=`which mpicc` \ -DCMAKE_Fortran_COMPILER=`which mpifort`\ -DPCMS_TIMEOUT=20 \ - -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_BUILD_TYPE=DEBUG \ -DCatch2_DIR=$DEPENDENCY_DIR/Catch2/install/lib64/cmake/Catch2/ \ -Dmeshfields_DIR=$DEPENDENCY_DIR/${DEVICE_ARCH}/meshFields/install/lib64/cmake/meshfields \ -DOmega_h_DIR=$DEPENDENCY_DIR/${DEVICE_ARCH}/omega_h/install/lib64/cmake/Omega_h/ \ @@ -116,6 +117,7 @@ jobs: -DKokkos_DIR=$DEPENDENCY_DIR/${DEVICE_ARCH}/kokkos/install/lib64/cmake/Kokkos/ \ -DKokkosKernels_DIR=$DEPENDENCY_DIR/${DEVICE_ARCH}/kokkos-kernels/install/lib64/cmake/KokkosKernels/ \ -DBUILD_TESTING=ON \ + -DPCMS_ENABLE_Python=ON \ -DPCMS_ENABLE_PETSC=ON \ -DPETSC_ARCH="" \ -DPETSC_DIR=$DEPENDENCY_DIR/${DEVICE_ARCH}/petsc/install \ @@ -124,6 +126,10 @@ jobs: cmake --build $bdir -j 4 ctest --test-dir $bdir --output-on-failure + export PYTHONPATH=$DEPENDENCY_DIR/${DEVICE_ARCH}/omega_h/install/lib64/python3.11/site-packages:$bdir/src/pcms/pythonapi:$PYTHONPATH + export LD_LIBRARY_PATH=$DEPENDENCY_DIR/${DEVICE_ARCH}/omega_h/install/lib64:$LD_LIBRARY_PATH + pytest -v -p no:cacheprovider ${{github.workspace}}/pcms_${{ github.event.id }}/test/python + - name: Save Result Link if: ${{ !cancelled() }} #prepare report unless the job was cancelled diff --git a/pyproject.toml b/pyproject.toml index 949bcbcc..25dc90ec 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,3 +8,10 @@ version = "0.4.0" [tool.scikit-build.cmake.define] PCMS_ENABLE_Python = true + +[tool.pytest.ini_options] +testpaths = ["test/python"] +python_files = "test_*.py" +python_classes = "Test*" +python_functions = "test_*" +addopts = "--strict-markers" diff --git a/src/pcms/pythonapi/test_field_copy.py b/src/pcms/pythonapi/test_field_copy.py deleted file mode 100644 index 8609ea4f..00000000 --- a/src/pcms/pythonapi/test_field_copy.py +++ /dev/null @@ -1,83 +0,0 @@ -import pcms -import PyOmega_h as omega_h -import numpy as np - -def test_copy(world, dim, order, num_components): - """Test copying omega_h field data.""" - nx = 100 - ny = 100 if dim > 1 else 0 - nz = 100 if dim > 2 else 0 - print(f"\nStarting test: dim={dim}, order={order}, num_components={num_components}") - - # Build mesh - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, - nx, ny, nz, - False - ) - print(f" Built {dim}D mesh with {nx}x{ny}x{nz} elements") - print(f" Mesh type: {type(mesh)}") - print(f" Mesh object: {mesh}") - - # Create factory and layout - print(" About to create factory...") - factory = pcms.LagrangeFunctionSpace.from_mesh( - mesh, order, num_components, pcms.CoordinateSystem.Cartesian - ) - print(f"Testing dim={dim}, order={order}, num_components={num_components}...") - - # Create original field and set data - original = factory.create_field() - print(" Created original field") - ndata = original.get_num_dof_holders() * original.get_num_components() - print(f" Number of data points: {ndata}") - - # Create sequential array of IDs - ids = np.arange(ndata, dtype=np.float64) - print(f" Created array of IDs from 0 to {ndata-1}") - original.set_dof_holder_data(ids) - print(" Set data in original field") - - # Create copied field and copy data - copied = factory.create_field() - copier = pcms.Copy(factory, factory) - copier.apply(original, copied) - print(" Copied data to new field") - - # Get copied data - copied_array = copied.get_dof_holder_data() - - # Verify the copy - assert len(copied_array) == ndata, f"Expected {ndata} elements, got {len(copied_array)}" - - # Check that all values match - differences = np.abs(ids - copied_array) - num_matches = np.sum(differences < 1e-12) - - assert num_matches == ndata, f"Only {num_matches}/{ndata} elements matched" - - print(f"✓ Test passed: dim={dim}, order={order}, num_components={num_components}") - -def main(): - """Run all test cases""" - print("Testing copy omega_h field data...") - - # Initialize Omega_h library - lib = omega_h.OmegaHLibrary() - world = lib.world() - print("Initialized Omega_h library and world") - - # Run test cases - test_copy(world, 2, 1, 1) - print("first passed") - test_copy(world, 2, 2, 1) - - print("\nAll tests passed!") - - del world - del lib - -if __name__ == "__main__": - main() \ No newline at end of file diff --git a/src/pcms/pythonapi/test_file_io.py b/src/pcms/pythonapi/test_file_io.py deleted file mode 100644 index 7366af04..00000000 --- a/src/pcms/pythonapi/test_file_io.py +++ /dev/null @@ -1,442 +0,0 @@ -#!/usr/bin/env python3 -""" -Test Omega_h file I/O Python bindings -""" - -import pcms -import PyOmega_h as omega_h -import numpy as np -import os -import shutil -import gc -import tempfile - -def test_binary_io(lib, world): - """Test binary file format I/O""" - print("\n=== Testing Binary Format I/O ===") - - # Build a simple 3D box mesh - print("Creating test mesh...") - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, # x, y, z dimensions - 4, 4, 4, # nx, ny, nz divisions - False # symmetric - ) - - # Get initial mesh properties - nverts_orig = mesh.nverts() - nelems_orig = mesh.nelems() - dim_orig = mesh.dim() - print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") - - # Create unique temporary directory for test files - test_dir = tempfile.mkdtemp(prefix="pcms_test_binary_") - print(f"Using temporary directory: {test_dir}") - try: - # Test binary write/read - binary_file = os.path.join(test_dir, "test_mesh.osh") - print(f"Writing to binary file: {binary_file}") - omega_h.write_mesh_binary(binary_file, mesh) - - print(f"Reading from binary file: {binary_file}") - mesh_read = omega_h.read_mesh_binary(binary_file, lib) - - # Verify mesh properties - assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" - assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" - assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" - - print("✓ Binary I/O test passed") - - # Explicitly delete mesh objects before cleanup - del mesh_read - del mesh - gc.collect() # Force garbage collection - - finally: - # Clean up temporary files - shutil.rmtree(test_dir, ignore_errors=True) - print(f"Cleaned up test directory: {test_dir}") - - -def test_gmsh_io(lib, world): - """Test Gmsh file format I/O""" - print("\n=== Testing Gmsh Format I/O ===") - - # Build a simple 2D box mesh - print("Creating test mesh...") - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 2.0, 1.0, 0.0, # x, y, z dimensions (z=0 for 2D) - 5, 3, 0, # nx, ny, nz divisions (nz=0 for 2D) - False # symmetric - ) - - # Get initial mesh properties - nverts_orig = mesh.nverts() - nelems_orig = mesh.nelems() - dim_orig = mesh.dim() - print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") - - # Create unique temporary directory for test files - test_dir = tempfile.mkdtemp(prefix="pcms_test_gmsh_") - print(f"Using temporary directory: {test_dir}") - try: - # Test Gmsh write/read - gmsh_file = os.path.join(test_dir, "test_mesh.msh") - print(f"Writing to Gmsh file: {gmsh_file}") - omega_h.write_mesh_gmsh(gmsh_file, mesh) - - print(f"Reading from Gmsh file: {gmsh_file}") - mesh_read = omega_h.read_mesh_gmsh(gmsh_file, world) - - # Verify mesh properties - assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" - assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" - assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" - - print("✓ Gmsh I/O test passed") - - # Explicitly delete mesh objects before cleanup - del mesh_read - del mesh - gc.collect() # Force garbage collection - - finally: - # Clean up temporary files - shutil.rmtree(test_dir, ignore_errors=True) - print(f"Cleaned up test directory: {test_dir}") - - -def test_vtk_io(lib, world): - """Test VTK file format I/O""" - print("\n=== Testing VTK Format I/O ===") - - # Build a simple 3D box mesh - print("Creating test mesh...") - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.5, 1.0, 0.5, # x, y, z dimensions - 3, 3, 2, # nx, ny, nz divisions - False # symmetric - ) - - # Get initial mesh properties - nverts_orig = mesh.nverts() - nelems_orig = mesh.nelems() - dim_orig = mesh.dim() - coords_orig = mesh.coords() - print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") - print(f"Coordinates shape: {coords_orig.shape}") - - # Create unique temporary directory for test files - test_dir = tempfile.mkdtemp(prefix="pcms_test_vtk_") - print(f"Using temporary directory: {test_dir}") - try: - # Test VTU write (VTU is write-only in the API, typically for visualization) - vtu_file = os.path.join(test_dir, "test_mesh.vtu") - print(f"Writing to VTU file: {vtu_file}") - omega_h.write_mesh_vtu(vtu_file, mesh, compress=False) - - # Verify file was created - assert os.path.exists(vtu_file), f"VTU file was not created: {vtu_file}" - file_size = os.path.getsize(vtu_file) - print(f"VTU file created successfully (size: {file_size} bytes)") - - # Test compressed VTU write - vtu_compressed = os.path.join(test_dir, "test_mesh_compressed.vtu") - print(f"Writing compressed VTU file: {vtu_compressed}") - omega_h.write_mesh_vtu(vtu_compressed, mesh, compress=True) - - assert os.path.exists(vtu_compressed), f"Compressed VTU file was not created" - compressed_size = os.path.getsize(vtu_compressed) - print(f"Compressed VTU file created (size: {compressed_size} bytes)") - - print("✓ VTK I/O test passed") - - # Explicitly delete mesh objects before cleanup - del mesh - gc.collect() # Force garbage collection - - finally: - # Clean up temporary files - shutil.rmtree(test_dir, ignore_errors=True) - print(f"Cleaned up test directory: {test_dir}") - - -def test_meshb_io(lib, world): - """Test MESHB file format I/O""" - print("\n=== Testing MESHB Format I/O ===") - - # Check if MESHB support is available - if not hasattr(omega_h, 'write_mesh_meshb'): - print("⊘ MESHB support not available (OMEGA_H_USE_LIBMESHB not enabled)") - return - - # Build a simple 3D box mesh - print("Creating test mesh...") - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, # x, y, z dimensions - 3, 3, 3, # nx, ny, nz divisions - False # symmetric - ) - - # Get initial mesh properties - nverts_orig = mesh.nverts() - nelems_orig = mesh.nelems() - dim_orig = mesh.dim() - print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") - - # Create unique temporary directory for test files - test_dir = tempfile.mkdtemp(prefix="pcms_test_meshb_") - print(f"Using temporary directory: {test_dir}") - try: - # Test MESHB write/read - meshb_file = os.path.join(test_dir, "test_mesh.mesh") - print(f"Writing to MESHB file: {meshb_file}") - omega_h.write_mesh_meshb(mesh, meshb_file, version=2) - - print(f"Reading from MESHB file: {meshb_file}") - # MESHB read requires pre-created mesh object - mesh_read = omega_h.Mesh(lib) - mesh_read.set_comm(world) - omega_h.read_mesh_meshb(mesh_read, meshb_file) - - # Verify mesh properties - assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" - assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" - assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" - - print("✓ MESHB I/O test passed") - - # Explicitly delete mesh objects before cleanup - del mesh_read - del mesh - gc.collect() # Force garbage collection - - finally: - # Clean up temporary files - shutil.rmtree(test_dir, ignore_errors=True) - print(f"Cleaned up test directory: {test_dir}") - - -def test_exodus_io(lib, world): - """Test Exodus file format I/O""" - print("\n=== Testing Exodus Format I/O ===") - - # Check if Exodus support is available - if not hasattr(omega_h, 'write_mesh_exodus'): - print("⊘ Exodus support not available (OMEGA_H_USE_SEACASEXODUS not enabled)") - return - - # Build a simple 3D box mesh - print("Creating test mesh...") - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, # x, y, z dimensions - 3, 3, 3, # nx, ny, nz divisions - False # symmetric - ) - - # Get initial mesh properties - nverts_orig = mesh.nverts() - nelems_orig = mesh.nelems() - dim_orig = mesh.dim() - print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") - - # Create unique temporary directory for test files - test_dir = tempfile.mkdtemp(prefix="pcms_test_exodus_") - print(f"Using temporary directory: {test_dir}") - try: - # Test Exodus write - exodus_file = os.path.join(test_dir, "test_mesh.exo") - print(f"Writing to Exodus file: {exodus_file}") - omega_h.write_mesh_exodus(exodus_file, mesh, verbose=False) - - # Test Exodus file handle API - if hasattr(omega_h, 'exodus_open'): - print(f"Testing Exodus file handle API with: {exodus_file}") - - # Open the file - exo_handle = omega_h.exodus_open(exodus_file, verbose=False) - print(f"Opened Exodus file with handle: {exo_handle}") - - # Get number of time steps - num_steps = omega_h.exodus_get_num_time_steps(exo_handle) - print(f"Number of time steps: {num_steps}") - - # Read mesh using file handle - mesh_from_handle = omega_h.OmegaHMesh(lib) - mesh_from_handle.set_comm(world) - omega_h.read_mesh_exodus(exo_handle, mesh_from_handle, verbose=False) - - # Verify mesh properties - assert mesh_from_handle.dim() == dim_orig, f"Dimension mismatch: {mesh_from_handle.dim()} != {dim_orig}" - assert mesh_from_handle.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_from_handle.nverts()} != {nverts_orig}" - assert mesh_from_handle.nelems() == nelems_orig, f"Element count mismatch: {mesh_from_handle.nelems()} != {nelems_orig}" - - # Close the file - omega_h.exodus_close(exo_handle) - # Explicitly delete mesh object from handle - del mesh_from_handle - - print("✓ Exodus I/O test passed") - - # Explicitly delete mesh objects before cleanup - del mesh - gc.collect() # Force garbage collection - - finally: - # Clean up temporary files - shutil.rmtree(test_dir, ignore_errors=True) - print(f"Cleaned up test directory: {test_dir}") - - -def test_adios2_io(lib, world): - """Test ADIOS2 file format I/O""" - print("\n=== Testing ADIOS2 Format I/O ===") - - # Check if ADIOS2 support is available - if not hasattr(omega_h, 'write_mesh_adios2'): - print("⊘ ADIOS2 support not available (OMEGA_H_USE_ADIOS2 not enabled)") - return - - # Build a simple 3D box mesh - print("Creating test mesh...") - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, # x, y, z dimensions - 3, 3, 3, # nx, ny, nz divisions - False # symmetric - ) - - # Get initial mesh properties - nverts_orig = mesh.nverts() - nelems_orig = mesh.nelems() - dim_orig = mesh.dim() - print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") - - # Create unique temporary directory for test files - test_dir = tempfile.mkdtemp(prefix="pcms_test_adios2_") - print(f"Using temporary directory: {test_dir}") - try: - # Test ADIOS2 write/read - adios2_file = os.path.join(test_dir, "test_mesh.bp") - print(f"Writing to ADIOS2 file: {adios2_file}") - omega_h.write_mesh_adios2(adios2_file, mesh, prefix="") - - print(f"Reading from ADIOS2 file: {adios2_file}") - mesh_read = omega_h.read_mesh_adios2(adios2_file, lib, prefix="") - - # Verify mesh properties - assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" - assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" - assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" - - print("✓ ADIOS2 I/O test passed") - - # Explicitly delete mesh objects before cleanup - del mesh_read - del mesh - gc.collect() # Force garbage collection - - finally: - # Clean up temporary files - shutil.rmtree(test_dir, ignore_errors=True) - print(f"Cleaned up test directory: {test_dir}") - - -def test_read_mesh_file_auto_detect(lib, world): - """Test automatic format detection with read_mesh_file""" - print("\n=== Testing Automatic Format Detection ===") - - # Build a test mesh - print("Creating test mesh...") - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, - 3, 3, 3, - False - ) - - nverts_orig = mesh.nverts() - nelems_orig = mesh.nelems() - - # Create unique temporary directory for test files - test_dir = tempfile.mkdtemp(prefix="pcms_test_autodetect_") - print(f"Using temporary directory: {test_dir}") - try: - # Write in different formats - binary_file = os.path.join(test_dir, "mesh.osh") - gmsh_file = os.path.join(test_dir, "mesh.msh") - - omega_h.write_mesh_binary(binary_file, mesh) - omega_h.write_mesh_gmsh(gmsh_file, mesh) - - # Test auto-detection for binary format - print(f"Auto-detecting binary file: {binary_file}") - mesh_binary = omega_h.read_mesh_file(binary_file, world) - assert mesh_binary.nverts() == nverts_orig - assert mesh_binary.nelems() == nelems_orig - print("✓ Binary auto-detection passed") - - # Test auto-detection for Gmsh format - print(f"Auto-detecting Gmsh file: {gmsh_file}") - mesh_gmsh = omega_h.read_mesh_file(gmsh_file, world) - assert mesh_gmsh.nverts() == nverts_orig - assert mesh_gmsh.nelems() == nelems_orig - print("✓ Gmsh auto-detection passed") - - # Explicitly delete mesh objects before cleanup - del mesh_binary - del mesh_gmsh - del mesh - gc.collect() # Force garbage collection - - finally: - # Clean up temporary files - shutil.rmtree(test_dir, ignore_errors=True) - print(f"Cleaned up test directory: {test_dir}") - -if __name__ == "__main__": - print("=" * 60) - print("Omega_h File I/O Python Binding Tests") - print("=" * 60) - - # Create library and world communicator once for all tests - lib = omega_h.OmegaHLibrary() - world = lib.world() - - try: - test_binary_io(lib, world) - test_gmsh_io(lib, world) - test_vtk_io(lib, world) - test_meshb_io(lib, world) - test_exodus_io(lib, world) - test_adios2_io(lib, world) - test_read_mesh_file_auto_detect(lib, world) - - print("\n" + "=" * 60) - print("✓ All tests passed!") - print("=" * 60) - - except Exception as e: - print("\n" + "=" * 60) - print(f"✗ Test failed with error: {e}") - print("=" * 60) - import traceback - traceback.print_exc() - exit(1) - finally: - # Explicitly delete objects to avoid an MPI finalizing issue - del world - del lib \ No newline at end of file diff --git a/src/pcms/pythonapi/test_mls_interpolation.py b/src/pcms/pythonapi/test_mls_interpolation.py deleted file mode 100644 index df3f116b..00000000 --- a/src/pcms/pythonapi/test_mls_interpolation.py +++ /dev/null @@ -1,108 +0,0 @@ -#!/usr/bin/env python3 -""" -Tests for MLS-backed PolynomialReconstructionFunctionSpace using the -Field-based API. -""" - -import numpy as np -import pcms -import PyOmega_h as omega_h - - -def poly_value(x, y, degree): - if degree == 0: - return 3.0 - if degree == 1: - return x + y - if degree == 2: - return x**2 + y**2 - if degree == 3: - return x**3 + y**3 - raise ValueError(f"Unsupported polynomial degree: {degree}") - - -def test_mls_interpolation_polynomial_reproduction(): - """ - Verify that MLS reproduces polynomials up to the configured degree. - - Source points are a 4x4 grid on [0,1]^2. Target points are six - arbitrary interior locations. For each interpolation degree d, MLS - must exactly reproduce any polynomial of degree <= d (up to the - configured tolerance). - """ - tolerance = 5e-4 - - grid_vals = np.linspace(0.0, 1.0, 4) - source_xy = np.array([(x, y) for x in grid_vals for y in grid_vals], - dtype=np.float64) - target_xy = np.array( - [ - [0.11, 0.19], - [0.23, 0.77], - [0.51, 0.49], - [0.87, 0.26], - [0.69, 0.91], - [0.35, 0.62], - ], - dtype=np.float64, - ) - num_targets = target_xy.shape[0] - - print("Testing MLS polynomial reproduction via PolynomialReconstructionFunctionSpace...") - - for interp_degree in range(1, 4): - # Use a radius that covers the entire unit square from any target so - # all 16 source points are always in support (equivalent to the - # full-support test in the old low-level API). - opts = pcms.MLSOptions() - opts.degree = interp_degree - opts.radius = 1.5 - opts.adapt_radius = False - opts.basis = pcms.RadialBasisFunction.NO_OP - - factory = pcms.PolynomialReconstructionFunctionSpace.from_coords( - source_xy, pcms.CoordinateSystem.Cartesian, opts - ) - field = factory.create_field() - request = pcms.EvaluationRequest.from_coordinates(target_xy) - evaluator = factory.create_point_evaluator(request) - results = np.zeros((num_targets, 1), dtype=np.float64) - - for func_degree in range(interp_degree, -1, -1): - source_values = np.array( - [poly_value(x, y, func_degree) for (x, y) in source_xy], - dtype=np.float64, - ) - exact_target_values = np.array( - [poly_value(x, y, func_degree) for (x, y) in target_xy], - dtype=np.float64, - ) - - field.set_dof_holder_data(source_values) - evaluator.evaluate(field, results) - approx_target_values = results[:, 0] - - max_abs_err = np.max(np.abs(exact_target_values - approx_target_values)) - assert max_abs_err < tolerance, ( - f"MLS failed for interp_degree={interp_degree}, " - f"func_degree={func_degree}: max_abs_err={max_abs_err}" - ) - - print("MLS polynomial reproduction test passed.") - - -if __name__ == "__main__": - lib = omega_h.OmegaHLibrary() - world = lib.world() - - try: - test_mls_interpolation_polynomial_reproduction() - print("MLS interpolation test passed") - except Exception as e: - print(f"MLS interpolation test failed: {e}") - import traceback - traceback.print_exc() - exit(1) - finally: - del world - del lib \ No newline at end of file diff --git a/src/pcms/pythonapi/test_omega_h_field.py b/src/pcms/pythonapi/test_omega_h_field.py deleted file mode 100644 index 4ef1db70..00000000 --- a/src/pcms/pythonapi/test_omega_h_field.py +++ /dev/null @@ -1,339 +0,0 @@ -#!/usr/bin/env python3 -""" -Field-centric Python examples for Omega_h-backed function spaces. -""" - -import numpy as np -import pcms -import PyOmega_h as omega_h - - -def test_field_methods(world, dim, order, num_components): - """Create an Omega_h-backed Field and exercise the public Field API.""" - nx = 10 - ny = 10 if dim > 1 else 0 - nz = 10 if dim > 2 else 0 - - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, - nx, ny, nz, - False - ) - - factory = pcms.LagrangeFunctionSpace.from_mesh( - mesh, order, num_components, pcms.CoordinateSystem.Cartesian - ) - field = factory.create_field() - - assert field.get_num_components() == num_components - assert field.get_num_dof_holders() > 0 - - coords = field.get_dof_holder_coordinates() - assert coords.shape[0] == field.get_num_dof_holders() - assert coords.shape[1] == dim - - ndata = field.get_num_dof_holders() * num_components - test_data = np.arange(ndata, dtype=np.float64) - field.set_dof_holder_data(test_data) - np.testing.assert_allclose(field.get_dof_holder_data(), test_data) - - -def test_field_transfer(world, dim, order, num_components): - """Transfer data between identical Omega_h-backed function spaces.""" - if num_components != 1: - return - - nx = 10 - ny = 10 if dim > 1 else 0 - nz = 10 if dim > 2 else 0 - - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, - nx, ny, nz, - False - ) - - source_space = pcms.LagrangeFunctionSpace.from_mesh( - mesh, order, num_components, pcms.CoordinateSystem.Cartesian - ) - target_space = pcms.LagrangeFunctionSpace.from_mesh( - mesh, order, num_components, pcms.CoordinateSystem.Cartesian - ) - - source = source_space.create_field() - target = target_space.create_field() - - coords = source.get_dof_holder_coordinates() - source_data = np.zeros(source.get_num_dof_holders(), dtype=np.float64) - for i in range(source.get_num_dof_holders()): - source_data[i] = np.sum(coords[i, :dim]) - source.set_dof_holder_data(source_data) - - interp = pcms.Interpolator(source_space, target_space) - interp.apply(source, target) - - np.testing.assert_allclose(target.get_dof_holder_data(), source_data, atol=1e-14) - - -def test_field_evaluation(world, dim, order, num_components): - """Evaluate an Omega_h-backed field at explicit query points.""" - nx = 10 - ny = 10 if dim > 1 else 0 - nz = 10 if dim > 2 else 0 - - print( - f"\nTesting field evaluation: dim={dim}, order={order}, " - f"num_components={num_components}" - ) - - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, - nx, ny, nz, - False - ) - - factory = pcms.LagrangeFunctionSpace.from_mesh( - mesh, order, num_components, pcms.CoordinateSystem.Cartesian - ) - field = factory.create_field() - - print(" Setting up field data...") - mesh_coords = mesh.coords() - num_verts = mesh.nverts() - print(f" Mesh has {num_verts} vertices") - - num_owned = field.get_num_dof_holders() - ndata = num_owned * num_components - print(f" Setting up field data for {ndata} DOFs") - - def test_func(x, y, z, component): - return np.sin(5.0 * x * y) + float(component) - - test_data = np.zeros(ndata, dtype=np.float64) - for i in range(min(num_verts, num_owned)): - x = mesh_coords[i * dim + 0] if dim >= 1 else 0.0 - y = mesh_coords[i * dim + 1] if dim >= 2 else 0.0 - z = mesh_coords[i * dim + 2] if dim >= 3 else 0.0 - for c in range(num_components): - idx = i * num_components + c - test_data[idx] = test_func(x, y, z, c) - - if order == 2 and num_owned > num_verts: - for i in range(num_verts, num_owned): - for c in range(num_components): - idx = i * num_components + c - test_data[idx] = float(c) + 0.5 - - field.set_dof_holder_data(test_data) - print(" Set field data with test function") - - if dim == 2: - eval_coords = np.array([ - [0.5, 0.5], - [0.25, 0.25], - [0.75, 0.75], - [0.1, 0.9], - [0.9, 0.1], - ], dtype=np.float64) - else: - eval_coords = np.array([ - [0.5, 0.5, 0.5], - [0.25, 0.25, 0.25], - [0.75, 0.75, 0.75], - [0.1, 0.9, 0.1], - [0.9, 0.1, 0.9], - ], dtype=np.float64) - - print(f" Evaluating at {eval_coords.shape[0]} points") - request = pcms.EvaluationRequest.from_coordinates(eval_coords) - evaluator = factory.create_point_evaluator(request) - print(" Created point evaluator") - - eval_values = np.zeros((eval_coords.shape[0], num_components), - dtype=np.float64) - evaluator.evaluate(field, eval_values) - print(" Field evaluated successfully") - - for i in range(eval_coords.shape[0]): - coords_str = ", ".join([f"{eval_coords[i, j]:.3f}" for j in range(dim)]) - values_str = ", ".join( - [f"{eval_values[i, c]:.4f}" for c in range(num_components)] - ) - print(f" Point {i} ({coords_str}): values = [{values_str}]") - for c in range(num_components): - val = eval_values[i, c] - assert not np.isnan(val) - assert not np.isinf(val) - - -def test_tag_operations(world, dim): - """Exercise Omega_h mesh tag creation, mutation, and query helpers.""" - nx = 5 - ny = 5 if dim > 1 else 0 - nz = 5 if dim > 2 else 0 - - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, - nx, ny, nz, - False - ) - rng = np.random.default_rng(0) - - mesh.add_tag(0, "test_float", 1, dtype="float64") - mesh.add_tag(0, "test_int32", 1, dtype="int32") - mesh.add_tag(0, "test_int64", 1, dtype="int64") - - assert mesh.has_tag(0, "test_float") - assert mesh.has_tag(0, "test_int32") - assert mesh.has_tag(0, "test_int64") - - nverts = mesh.nverts() - temp_data = np.linspace(0.0, 100.0, nverts, dtype=np.float64) - ids_data = np.arange(nverts, dtype=np.int32) - velocity_data = rng.standard_normal(nverts * dim).astype(np.float64) - - mesh.add_tag(0, "temperature", 1, temp_data) - mesh.add_tag(0, "vertex_ids", 1, ids_data) - mesh.add_tag(0, "velocity", dim, velocity_data) - - np.testing.assert_allclose(mesh.get_tag(0, "temperature"), temp_data) - np.testing.assert_array_equal(mesh.get_tag(0, "vertex_ids"), ids_data) - np.testing.assert_allclose(mesh.get_tag(0, "velocity"), velocity_data) - - new_temp = np.ones(nverts, dtype=np.float64) * 50.0 - mesh.set_tag(0, "temperature", new_temp) - np.testing.assert_allclose(mesh.get_tag(0, "temperature"), new_temp) - - if dim == 3: - ncomps = 6 - stress_data = rng.standard_normal(mesh.nelems() * ncomps).astype( - np.float64) - mesh.add_tag(dim, "stress", ncomps, stress_data, - internal=False, - array_type=omega_h.ArrayType.SymmetricSquareMatrix) - assert mesh.has_tag(dim, "stress") - - internal_data = np.zeros(nverts, dtype=np.float64) - mesh.add_tag(0, "internal_temp", 1, internal_data, internal=True) - assert mesh.has_tag(0, "internal_temp") - - mesh.remove_tag(0, "test_float") - assert not mesh.has_tag(0, "test_float") - assert mesh.has_tag(0, "temperature") - assert mesh.ntags(0) > 0 - - nelems = mesh.nelems() - elem_quality = rng.random(nelems).astype(np.float64) - mesh.add_tag(dim, "quality", 1, elem_quality) - np.testing.assert_allclose(mesh.get_tag(dim, "quality"), elem_quality) - - if dim >= 2: - nedges = mesh.nedges() - edge_length = np.ones(nedges, dtype=np.float64) - mesh.add_tag(1, "edge_marker", 1, edge_length) - np.testing.assert_allclose(mesh.get_tag(1, "edge_marker"), edge_length) - - elem_verts = mesh.ask_elem_verts() - global_ids = mesh.globals(0) - verts_of_elems = mesh.ask_verts_of(dim) - owned_verts = mesh.owned(0) - - assert len(elem_verts) > 0 - assert len(global_ids) == nverts - assert len(verts_of_elems) > 0 - assert len(owned_verts) == nverts - assert np.sum(owned_verts) > 0 - assert isinstance(mesh.has_adj(0, dim), (bool, np.bool_)) - - -def test_entity_coordinates(world, dim): - """Exercise entity-coordinate and averaging helpers on Omega_h meshes.""" - nx = 4 - ny = 4 if dim > 1 else 0 - nz = 4 if dim > 2 else 0 - - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 1.0, - nx, ny, nz, - False - ) - - vertex_coords = mesh.coords() - assert len(vertex_coords) == mesh.nverts() * dim - - if dim >= 2: - edge_coords = omega_h.average_field(mesh, 1, dim, vertex_coords) - edge_coords_np = np.array(edge_coords).reshape(-1, dim) - assert edge_coords_np.shape[0] == mesh.nedges() - - if dim == 2: - elem_coords = omega_h.average_field(mesh, 2, dim, vertex_coords) - elem_coords_np = np.array(elem_coords).reshape(-1, dim) - assert elem_coords_np.shape[0] == mesh.nelems() - - if dim == 3: - face_coords = omega_h.average_field(mesh, 2, dim, vertex_coords) - face_coords_np = np.array(face_coords).reshape(-1, dim) - assert face_coords_np.shape[0] == mesh.nfaces() - - region_coords = omega_h.average_field(mesh, 3, dim, vertex_coords) - region_coords_np = np.array(region_coords).reshape(-1, dim) - assert region_coords_np.shape[0] == mesh.nregions() - - vertex_field = np.arange(mesh.nverts(), dtype=np.float64) - mesh.add_tag(0, "test_vertex_field", 1, vertex_field) - averaged_field = omega_h.average_field(mesh, dim, 1, vertex_field) - assert averaged_field.shape[0] == mesh.nents(dim) - - -def main(): - """Run all test cases""" - print("Testing Omega_h Field Python bindings...") - lib = omega_h.OmegaHLibrary() - world = lib.world() - print("Initialized Omega_h library and world") - - test_field_methods(world, 2, 1, 1) - test_field_methods(world, 2, 2, 1) - - print("\n" + "=" * 60) - print("Testing field evaluation...") - print("=" * 60) - test_field_evaluation(world, 2, 1, 1) - test_field_evaluation(world, 2, 2, 1) - - print("\n" + "=" * 60) - print("Testing field transfer...") - print("=" * 60) - test_field_transfer(world, 2, 1, 1) - test_field_transfer(world, 2, 2, 1) - - print("\n" + "=" * 60) - print("Testing tag operations...") - print("=" * 60) - test_tag_operations(world, 2) - test_tag_operations(world, 3) - - print("\n" + "=" * 60) - print("Testing entity coordinate computation...") - print("=" * 60) - test_entity_coordinates(world, 2) - test_entity_coordinates(world, 3) - - print("\n✓ All tests passed!") - del world - del lib - - -if __name__ == "__main__": - main() \ No newline at end of file diff --git a/src/pcms/pythonapi/test_uniform_grid_field.py b/src/pcms/pythonapi/test_uniform_grid_field.py deleted file mode 100644 index fe0745f3..00000000 --- a/src/pcms/pythonapi/test_uniform_grid_field.py +++ /dev/null @@ -1,344 +0,0 @@ -""" -Field-centric Python examples for uniform-grid-backed function spaces. -""" -import numpy as np -import pcms -import PyOmega_h as omega_h - - -def test_uniform_grid_field_creation(): - """Create a 2D field from a uniform-grid function space.""" - # Create a simple 2D structured grid. - grid = pcms.UniformGrid2D() - grid.bot_left = [0.0, 0.0] - grid.edge_length = [10.0, 10.0] - grid.divisions = [4, 4] - - # Build the function space from the grid, then create a Field from it. - factory = pcms.LagrangeFunctionSpace.from_uniform_grid( - grid, 1, pcms.CoordinateSystem.Cartesian - ) - field = factory.create_field() - - # Order-1 2D grids store data at vertices, so a 4x4 cell grid has 5x5 - # DOF holders. - expected_vertices = (grid.divisions[0] + 1) * (grid.divisions[1] + 1) - print(f"Grid cells: {grid.get_num_cells()}") - print(f"Num components: {field.get_num_components()}") - print(f"Num owned DOF holders: {field.get_num_dof_holders()}") - assert field.get_num_components() == 1 - assert field.get_num_dof_holders() == expected_vertices - print(f"Field created: {field is not None}") - - -def test_uniform_grid_field_data_operations(): - """Set and get flat DOF data through Field.""" - # A 2x2 cell grid has 3x3 vertex DOF holders. - grid = pcms.UniformGrid2D() - grid.bot_left = [0.0, 0.0] - grid.edge_length = [10.0, 10.0] - grid.divisions = [2, 2] - - factory = pcms.LagrangeFunctionSpace.from_uniform_grid( - grid, 1, pcms.CoordinateSystem.Cartesian - ) - field = factory.create_field() - - # Field data is set and retrieved as a flat 1D DOF array. - data = np.arange(field.get_num_dof_holders(), dtype=np.float64) - print(f"Number of vertices: {field.get_num_dof_holders()}") - print(f"Setting data: {data}") - field.set_dof_holder_data(data) - retrieved_data = field.get_dof_holder_data() - print(f"Retrieved data: {retrieved_data}") - np.testing.assert_allclose(retrieved_data, data) - print("Data successfully set and retrieved!") - - -def test_uniform_grid_field_coordinates_2d(): - """Expose 2D DOF-holder coordinates through Field.""" - # Coordinate queries should now flow through Field rather than layout - # objects. - grid = pcms.UniformGrid2D() - grid.bot_left = [0.0, 0.0] - grid.edge_length = [10.0, 10.0] - grid.divisions = [2, 3] - - factory = pcms.LagrangeFunctionSpace.from_uniform_grid( - grid, 1, pcms.CoordinateSystem.Cartesian - ) - field = factory.create_field() - - coords = field.get_dof_holder_coordinates() - expected_vertices = (grid.divisions[0] + 1) * (grid.divisions[1] + 1) - - assert isinstance(coords, np.ndarray) - assert coords.shape == (expected_vertices, 2) - np.testing.assert_allclose(coords[0], [0.0, 0.0]) - np.testing.assert_allclose(coords[-1], [10.0, 10.0]) - print("2D field coordinates verified!") - - -def test_uniform_grid_closest_cell(): - """UniformGrid helpers remain available for grid setup.""" - # Grid setup helpers remain public even though layout classes do not. - grid = pcms.UniformGrid2D() - grid.bot_left = [0.0, 0.0] - grid.edge_length = [10.0, 10.0] - grid.divisions = [4, 4] - - cell_id = grid.closest_cell_id(np.array([5.0, 5.0])) - cell_id_outside = grid.closest_cell_id(np.array([-1.0, -1.0])) - print(f"Point {[5.0, 5.0]} is in cell {cell_id}") - print(f"Point {[-1.0, -1.0]} (outside) maps to cell {cell_id_outside}") - assert cell_id >= 0 - assert cell_id_outside >= 0 - - -def test_uniform_grid_field_coordinates_3d(): - """Expose 3D DOF-holder coordinates and data through Field.""" - # 3D coordinate access and data operations follow the same pattern as 2D. - grid = pcms.UniformGrid3D() - grid.bot_left = [0.0, 0.0, 0.0] - grid.edge_length = [10.0, 10.0, 10.0] - grid.divisions = [2, 1, 3] - - factory = pcms.LagrangeFunctionSpace.from_uniform_grid( - grid, 1, pcms.CoordinateSystem.Cartesian - ) - field = factory.create_field() - - coords = field.get_dof_holder_coordinates() - expected_vertices = ((grid.divisions[0] + 1) * - (grid.divisions[1] + 1) * - (grid.divisions[2] + 1)) - - assert coords.shape == (expected_vertices, 3) - np.testing.assert_allclose(coords[0], [0.0, 0.0, 0.0]) - np.testing.assert_allclose(coords[-1], [10.0, 10.0, 10.0]) - print("3D field coordinates verified!") - - # Data set/get on a 3D grid field. - data = np.ones(expected_vertices, dtype=np.float64) * 42.0 - field.set_dof_holder_data(data) - np.testing.assert_allclose(field.get_dof_holder_data(), data) - print("3D field data set/get verified!") - - -def test_uniform_grid_field_evaluation(): - """Evaluate a uniform-grid field at explicit query points.""" - grid = pcms.UniformGrid2D() - grid.bot_left = [0.0, 0.0] - grid.edge_length = [10.0, 10.0] - grid.divisions = [2, 2] - - factory = pcms.LagrangeFunctionSpace.from_uniform_grid( - grid, 1, pcms.CoordinateSystem.Cartesian - ) - field = factory.create_field() - - # Set a linear field f(x,y) = x + y so the expected values are exact. - coords = field.get_dof_holder_coordinates() - data = np.array([coords[i, 0] + coords[i, 1] - for i in range(field.get_num_dof_holders())], - dtype=np.float64) - field.set_dof_holder_data(data) - - eval_coords = np.array([[2.5, 2.5], [7.5, 7.5]], dtype=np.float64) - request = pcms.EvaluationRequest.from_coordinates(eval_coords) - evaluator = factory.create_point_evaluator(request) - results = np.zeros((len(eval_coords), 1), dtype=np.float64) - evaluator.evaluate(field, results) - - print(f"Evaluation results: {results.flatten()}") - print(f"Expected (approximately): [5.0, 15.0]") - np.testing.assert_allclose(results.flatten(), [5.0, 15.0], atol=1e-10) - print("Uniform grid field evaluation verified!") - - -def test_uniform_grid_workflow(world): - """ - Test complete UniformGrid workflow with Omega_h mesh and field interpolation. - - This test: - 1. Creates an Omega_h mesh - 2. Creates a uniform grid from the mesh - 3. Creates a binary mask field - 4. Creates and initializes an Omega_h field with f(x,y) = x + 2*y - 5. Transfers the field to uniform grid via interpolation - 6. Verifies the transferred values - """ - # Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements - mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 0.0, - 4, 4, 0, - False - ) - - # Create a uniform grid that covers the mesh bounding box, then construct a - # binary mask field directly as a Field object. - grid = pcms.create_uniform_grid_from_mesh(mesh, [4, 4]) - print(f"Created uniform grid with {grid.get_num_cells()} cells") - mask_field = pcms.create_uniform_grid_binary_field(mesh, [4, 4]) - mask_data = mask_field.get_dof_holder_data() - print(f"Created mask field with {len(mask_data)} vertices") - - # Build a mesh-backed source field and initialize it with f(x,y)=x+2y. - omega_h_factory = pcms.LagrangeFunctionSpace.from_mesh(mesh, 1) - omega_h_field = omega_h_factory.create_field() - - coords = omega_h_field.get_dof_holder_coordinates() - omega_h_data = np.zeros(omega_h_field.get_num_dof_holders(), dtype=np.float64) - for i in range(omega_h_field.get_num_dof_holders()): - omega_h_data[i] = coords[i, 0] + 2.0 * coords[i, 1] - omega_h_field.set_dof_holder_data(omega_h_data) - print(f"Initialized Omega_h field with {omega_h_field.get_num_dof_holders()} nodes") - - ug_factory = pcms.LagrangeFunctionSpace.from_uniform_grid( - grid, 1, pcms.CoordinateSystem.Cartesian - ) - ug_field = ug_factory.create_field() - - # Transfer from the unstructured mesh field to the uniform-grid field. - interp = pcms.Interpolator(omega_h_factory, ug_factory) - interp.apply(omega_h_field, ug_field) - print("Field interpolation completed") - - ug_field_data = ug_field.get_dof_holder_data() - ug_coords = ug_field.get_dof_holder_coordinates() - - # The 4x4 cell grid should produce 5x5 vertex DOFs. - assert len(mask_data) == 25 - assert len(ug_field_data) == 25 - - # Verify interpolation at every target DOF holder. - num_errors = 0 - for vertex_id in range(len(ug_field_data)): - x = ug_coords[vertex_id, 0] - y = ug_coords[vertex_id, 1] - expected = x + 2.0 * y - actual = ug_field_data[vertex_id] - error = abs(actual - expected) - if error > 1e-10: - num_errors += 1 - print( - f"Vertex {vertex_id} at ({x:.2f}, {y:.2f}): " - f"expected {expected:.4f}, got {actual:.4f}, error {error:.2e}" - ) - assert error < 1e-10 - if num_errors == 0: - print("All uniform grid field values verified successfully!") - - -def test_omega_h_to_omega_h_transfer_workflow(world): - """ - Transfer a field from one Omega_h mesh to another Omega_h mesh. - - This preserves coverage for the workflow where both source and target use - mesh-backed function spaces rather than a uniform grid target. - """ - # Source mesh (coarser). - src_mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 0.0, - 4, 4, 0, - False - ) - - # Target mesh (finer). - tgt_mesh = omega_h.build_box( - world, - omega_h.Family.SIMPLEX, - 1.0, 1.0, 0.0, - 8, 8, 0, - False - ) - - # Create function spaces and fields on both meshes. - src_factory = pcms.LagrangeFunctionSpace.from_mesh(src_mesh, 1) - tgt_factory = pcms.LagrangeFunctionSpace.from_mesh(tgt_mesh, 1) - src_field = src_factory.create_field() - tgt_field = tgt_factory.create_field() - print("Created source and target Omega_h fields") - - # Initialize the source field with f(x,y)=x+2y. - src_coords = src_field.get_dof_holder_coordinates() - src_data = np.zeros(src_field.get_num_dof_holders()) - for i in range(len(src_data)): - src_data[i] = src_coords[i, 0] + 2.0 * src_coords[i, 1] - src_field.set_dof_holder_data(src_data) - print(f"Initialized source field with {len(src_data)} DOF values") - - # Transfer onto the target mesh and verify the target DOF values. - interp = pcms.Interpolator(src_factory, tgt_factory) - interp.apply(src_field, tgt_field) - print("Omega_h to Omega_h field interpolation completed") - - tgt_coords = tgt_field.get_dof_holder_coordinates() - tgt_data = tgt_field.get_dof_holder_data() - errors = 0 - for i in range(tgt_field.get_num_dof_holders()): - expected = tgt_coords[i, 0] + 2.0 * tgt_coords[i, 1] - error = abs(expected - tgt_data[i]) - if error > 1e-10: - errors += 1 - assert error < 1e-10 - print("Omega_h to Omega_h transfer verification completed") - - -def main(): - lib = omega_h.OmegaHLibrary() - world = lib.world() - - print("=" * 60) - print("Testing UniformGrid Field Creation") - print("=" * 60) - test_uniform_grid_field_creation() - - print("\n" + "=" * 60) - print("Testing UniformGrid Field Data Operations") - print("=" * 60) - test_uniform_grid_field_data_operations() - - print("\n" + "=" * 60) - print("Testing UniformGrid Field Coordinates (2D)") - print("=" * 60) - test_uniform_grid_field_coordinates_2d() - - print("\n" + "=" * 60) - print("Testing UniformGrid Field Coordinates (3D)") - print("=" * 60) - test_uniform_grid_field_coordinates_3d() - - print("\n" + "=" * 60) - print("Testing Closest Cell ID") - print("=" * 60) - test_uniform_grid_closest_cell() - - print("\n" + "=" * 60) - print("Testing UniformGrid Field Evaluation") - print("=" * 60) - test_uniform_grid_field_evaluation() - - print("\n" + "=" * 60) - print("Testing UniformGrid Workflow") - print("=" * 60) - test_uniform_grid_workflow(world) - - print("\n" + "=" * 60) - print("Testing Omega_h to Omega_h Transfer Workflow") - print("=" * 60) - test_omega_h_to_omega_h_transfer_workflow(world) - - print("\n" + "=" * 60) - print("All tests passed!") - print("=" * 60) - - del world - - -if __name__ == "__main__": - main() \ No newline at end of file diff --git a/test/python/__init__.py b/test/python/__init__.py new file mode 100644 index 00000000..ae78246e --- /dev/null +++ b/test/python/__init__.py @@ -0,0 +1 @@ +# tests/__init__.py diff --git a/test/python/conftest.py b/test/python/conftest.py new file mode 100644 index 00000000..8d936d98 --- /dev/null +++ b/test/python/conftest.py @@ -0,0 +1,27 @@ +""" +Shared pytest fixtures for pcms Python API tests. +""" + +import pytest +import PyOmega_h as omega_h + + +@pytest.fixture(scope="session") +def omega_h_lib(): + """Create Omega_h library instance. + + Session-scoped so the (potentially expensive) Omega_h initialization + happens exactly once per test run. + """ + lib = omega_h.OmegaHLibrary() + yield lib + + +@pytest.fixture(scope="session") +def world(omega_h_lib): + """Create Omega_h world communicator. + + Session-scoped; the world is obtained from the session-scoped library + and shared by all tests that need mesh-building capabilities. + """ + return omega_h_lib.world() diff --git a/test/python/test_field_copy.py b/test/python/test_field_copy.py new file mode 100644 index 00000000..6850b7d0 --- /dev/null +++ b/test/python/test_field_copy.py @@ -0,0 +1,60 @@ +"""Test copying Omega_h field data.""" + +import numpy as np +import pytest +import pcms +import PyOmega_h as omega_h + + +class TestFieldCopy: + """Tests for field data copying between Omega_h-backed fields.""" + + @pytest.mark.parametrize("dim, order, num_components", [ + (2, 1, 1), + (2, 2, 1), + ]) + def test_copy(self, world, dim, order, num_components): + """Test copying omega_h field data.""" + nx = 100 + ny = 100 if dim > 1 else 0 + nz = 100 if dim > 2 else 0 + + # Build mesh + mesh = omega_h.build_box( + world, + omega_h.Family.SIMPLEX, + 1.0, 1.0, 1.0, + nx, ny, nz, + False + ) + + # Create factory and layout + factory = pcms.LagrangeFunctionSpace.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian + ) + + # Create original field and set data + original = factory.create_field() + ndata = original.get_num_dof_holders() * original.get_num_components() + + # Create sequential array of IDs + ids = np.arange(ndata, dtype=np.float64) + original.set_dof_holder_data(ids) + + # Create copied field and copy data + copied = factory.create_field() + copier = pcms.Copy(factory, factory) + copier.apply(original, copied) + + # Get copied data + copied_array = copied.get_dof_holder_data() + + # Verify the copy + assert len(copied_array) == ndata, \ + f"Expected {ndata} elements, got {len(copied_array)}" + + # Check that all values match + differences = np.abs(ids - copied_array) + num_matches = np.sum(differences < 1e-12) + assert num_matches == ndata, \ + f"Only {num_matches}/{ndata} elements matched" diff --git a/test/python/test_mls_interpolation.py b/test/python/test_mls_interpolation.py new file mode 100644 index 00000000..a2ac218a --- /dev/null +++ b/test/python/test_mls_interpolation.py @@ -0,0 +1,89 @@ +""" +Tests for MLS-backed PolynomialReconstructionFunctionSpace. +""" + +import numpy as np +import pcms + + +def poly_value(x, y, degree): + """Evaluate a simple polynomial at (x, y).""" + if degree == 0: + return 3.0 + if degree == 1: + return x + y + if degree == 2: + return x**2 + y**2 + if degree == 3: + return x**3 + y**3 + raise ValueError(f"Unsupported polynomial degree: {degree}") + + +class TestMLSInterpolation: + """Tests for MLS polynomial reproduction.""" + + def test_polynomial_reproduction(self): + """ + Verify that MLS reproduces polynomials up to the configured degree. + + Source points are a 4x4 grid on [0,1]^2. Target points are six + arbitrary interior locations. For each interpolation degree d, MLS + must exactly reproduce any polynomial of degree <= d (up to the + configured tolerance). + """ + tolerance = 5e-4 + + grid_vals = np.linspace(0.0, 1.0, 4) + source_xy = np.array( + [(x, y) for x in grid_vals for y in grid_vals], + dtype=np.float64, + ) + target_xy = np.array( + [ + [0.11, 0.19], + [0.23, 0.77], + [0.51, 0.49], + [0.87, 0.26], + [0.69, 0.91], + [0.35, 0.62], + ], + dtype=np.float64, + ) + num_targets = target_xy.shape[0] + + for interp_degree in range(1, 4): + opts = pcms.MLSOptions() + opts.degree = interp_degree + opts.radius = 1.5 + opts.adapt_radius = False + opts.basis = pcms.RadialBasisFunction.NO_OP + + factory = pcms.PolynomialReconstructionFunctionSpace.from_coords( + source_xy, pcms.CoordinateSystem.Cartesian, opts + ) + field = factory.create_field() + request = pcms.EvaluationRequest.from_coordinates(target_xy) + evaluator = factory.create_point_evaluator(request) + results = np.zeros((num_targets, 1), dtype=np.float64) + + for func_degree in range(interp_degree, -1, -1): + source_values = np.array( + [poly_value(x, y, func_degree) for (x, y) in source_xy], + dtype=np.float64, + ) + exact_target_values = np.array( + [poly_value(x, y, func_degree) for (x, y) in target_xy], + dtype=np.float64, + ) + + field.set_dof_holder_data(source_values) + evaluator.evaluate(field, results) + approx_target_values = results[:, 0] + + max_abs_err = np.max( + np.abs(exact_target_values - approx_target_values) + ) + assert max_abs_err < tolerance, ( + f"MLS failed for interp_degree={interp_degree}, " + f"func_degree={func_degree}: max_abs_err={max_abs_err}" + ) diff --git a/test/python/test_omega_h_field.py b/test/python/test_omega_h_field.py new file mode 100644 index 00000000..f9c537c4 --- /dev/null +++ b/test/python/test_omega_h_field.py @@ -0,0 +1,140 @@ +""" +Field-centric tests for Omega_h-backed function spaces. +""" + +import numpy as np +import pytest +import pcms +import PyOmega_h as omega_h + +class TestOmegaHField: + """Tests for Omega_h-backed Field operations.""" + + @staticmethod + def _build_mesh(world, dim, nx=10): + """Build a box mesh of the requested dimension.""" + ny = nx if dim > 1 else 0 + nz = nx if dim > 2 else 0 + return omega_h.build_box( + world, omega_h.Family.SIMPLEX, 1.0, 1.0, 1.0, nx, ny, nz, False, + ) + + @pytest.mark.parametrize("dim, order, num_components", [ + (2, 1, 1), (2, 2, 1), + ]) + def test_field_methods(self, world, dim, order, num_components): + """Create an Omega_h-backed Field and exercise the public Field API.""" + mesh = self._build_mesh(world, dim) + factory = pcms.LagrangeFunctionSpace.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field() + + assert field.get_num_components() == num_components + assert field.get_num_dof_holders() > 0 + + coords = field.get_dof_holder_coordinates() + assert coords.shape[0] == field.get_num_dof_holders() + assert coords.shape[1] == dim + + ndata = field.get_num_dof_holders() * num_components + test_data = np.arange(ndata, dtype=np.float64) + expected_data = test_data.copy() + field.set_dof_holder_data(test_data) + test_data.fill(-999.0) + np.testing.assert_allclose(field.get_dof_holder_data(), expected_data) + + @pytest.mark.parametrize("dim, order, num_components", [ + (2, 1, 1), (2, 2, 1), + ]) + def test_field_transfer(self, world, dim, order, num_components): + """Transfer data between identical Omega_h-backed function spaces.""" + if num_components != 1: + pytest.skip("Multi-component transfer not yet supported") + + mesh = self._build_mesh(world, dim) + source_space = pcms.LagrangeFunctionSpace.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian + ) + target_space = pcms.LagrangeFunctionSpace.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian + ) + source = source_space.create_field() + target = target_space.create_field() + + coords = source.get_dof_holder_coordinates() + source_data = np.zeros(source.get_num_dof_holders(), dtype=np.float64) + for i in range(source.get_num_dof_holders()): + source_data[i] = np.sum(coords[i, :dim]) + source.set_dof_holder_data(source_data) + + interp = pcms.Interpolator(source_space, target_space) + interp.apply(source, target) + np.testing.assert_allclose( + target.get_dof_holder_data(), source_data, atol=1e-14 + ) + + @pytest.mark.parametrize("dim, order, num_components", [ + (2, 1, 1), (2, 2, 1), + ]) + def test_field_evaluation(self, world, dim, order, num_components): + """Evaluate an Omega_h-backed field at explicit query points.""" + mesh = self._build_mesh(world, dim) + factory = pcms.LagrangeFunctionSpace.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field() + + mesh_coords = mesh.coords() + num_verts = mesh.nverts() + num_owned = field.get_num_dof_holders() + ndata = num_owned * num_components + + def test_func(x, y, z, component): + return np.sin(5.0 * x * y) + float(component) + + test_data = np.zeros(ndata, dtype=np.float64) + for i in range(min(num_verts, num_owned)): + x = mesh_coords[i * dim + 0] if dim >= 1 else 0.0 + y = mesh_coords[i * dim + 1] if dim >= 2 else 0.0 + z = mesh_coords[i * dim + 2] if dim >= 3 else 0.0 + for c in range(num_components): + idx = i * num_components + c + test_data[idx] = test_func(x, y, z, c) + + if order == 2 and num_owned > num_verts: + for i in range(num_verts, num_owned): + for c in range(num_components): + idx = i * num_components + c + test_data[idx] = float(c) + 0.5 + + field.set_dof_holder_data(test_data) + + if dim == 2: + eval_coords = np.array([ + [0.5, 0.5], + [0.25, 0.25], + [0.75, 0.75], + [0.1, 0.9], + [0.9, 0.1], + ], dtype=np.float64) + else: + eval_coords = np.array([ + [0.5, 0.5, 0.5], + [0.25, 0.25, 0.25], + [0.75, 0.75, 0.75], + [0.1, 0.9, 0.1], + [0.9, 0.1, 0.9], + ], dtype=np.float64) + + request = pcms.EvaluationRequest.from_coordinates(eval_coords) + evaluator = factory.create_point_evaluator(request) + eval_values = np.zeros((eval_coords.shape[0], num_components), + dtype=np.float64) + evaluator.evaluate(field, eval_values) + + for i in range(eval_coords.shape[0]): + for c in range(num_components): + val = eval_values[i, c] + assert not np.isnan(val) + assert not np.isinf(val) diff --git a/test/python/test_uniform_grid_field.py b/test/python/test_uniform_grid_field.py new file mode 100644 index 00000000..4864cd9c --- /dev/null +++ b/test/python/test_uniform_grid_field.py @@ -0,0 +1,208 @@ +""" +Field-centric tests for uniform-grid-backed function spaces. +""" + +import numpy as np +import pytest +import pcms +import PyOmega_h as omega_h + +class TestUniformGridField: + """Tests for UniformGrid-backed Field operations (no Omega_h needed).""" + + def test_field_creation(self): + """Create a 2D field from a uniform-grid function space.""" + grid = pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [4, 4] + + factory = pcms.LagrangeFunctionSpace.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field() + expected = (grid.divisions[0] + 1) * (grid.divisions[1] + 1) + assert field.get_num_components() == 1 + assert field.get_num_dof_holders() == expected + + def test_data_operations(self): + """Set and get flat DOF data through Field.""" + grid = pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [2, 2] + + factory = pcms.LagrangeFunctionSpace.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field() + data = np.arange(field.get_num_dof_holders(), dtype=np.float64) + expected_data = data.copy() + field.set_dof_holder_data(data) + data.fill(-999.0) + np.testing.assert_allclose(field.get_dof_holder_data(), expected_data) + + def test_coordinates_2d(self): + """Expose 2D DOF-holder coordinates through Field.""" + grid = pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [2, 3] + + factory = pcms.LagrangeFunctionSpace.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field() + coords = field.get_dof_holder_coordinates() + expected = (grid.divisions[0] + 1) * (grid.divisions[1] + 1) + + assert isinstance(coords, np.ndarray) + assert coords.shape == (expected, 2) + np.testing.assert_allclose(coords[0], [0.0, 0.0]) + np.testing.assert_allclose(coords[-1], [10.0, 10.0]) + + def test_closest_cell(self): + """UniformGrid helpers remain available for grid setup.""" + grid = pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [4, 4] + + cell_id = grid.closest_cell_id(np.array([5.0, 5.0])) + cell_id_outside = grid.closest_cell_id(np.array([-1.0, -1.0])) + assert isinstance(cell_id, (int, np.integer)) + assert isinstance(cell_id_outside, (int, np.integer)) + + def test_coordinates_3d(self): + """Expose 3D DOF-holder coordinates and data through Field.""" + grid = pcms.UniformGrid3D() + grid.bot_left = [0.0, 0.0, 0.0] + grid.edge_length = [10.0, 10.0, 10.0] + grid.divisions = [2, 2, 2] + + factory = pcms.LagrangeFunctionSpace.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field() + coords = field.get_dof_holder_coordinates() + expected = ( + (grid.divisions[0] + 1) + * (grid.divisions[1] + 1) + * (grid.divisions[2] + 1) + ) + assert coords.shape == (expected, 3) + + data = np.arange(expected, dtype=np.float64) + expected_data = data.copy() + field.set_dof_holder_data(data) + data.fill(-999.0) + np.testing.assert_allclose(field.get_dof_holder_data(), expected_data) + + def test_field_evaluation(self): + """Evaluate a uniform-grid field at explicit query points.""" + grid = pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [1.0, 1.0] + grid.divisions = [10, 10] + + factory = pcms.LagrangeFunctionSpace.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field() + coords = field.get_dof_holder_coordinates() + num_dofs = field.get_num_dof_holders() + data = np.zeros(num_dofs, dtype=np.float64) + for i in range(num_dofs): + data[i] = coords[i, 0] + 2.0 * coords[i, 1] + field.set_dof_holder_data(data) + + query_pts = np.array( + [[0.1, 0.2], [0.5, 0.5], [0.9, 0.8]], dtype=np.float64 + ) + request = pcms.EvaluationRequest.from_coordinates(query_pts) + evaluator = factory.create_point_evaluator(request) + results = np.zeros((len(query_pts), 1), dtype=np.float64) + evaluator.evaluate(field, results) + + for i, pt in enumerate(query_pts): + expected = pt[0] + 2.0 * pt[1] + assert np.abs(results[i, 0] - expected) < 1e-6 + + +class TestUniformGridOmegaHWorkflow: + """Tests combining UniformGrid with Omega_h meshes (needs world).""" + + def test_uniform_grid_to_omega_h_workflow(self, world): + """ + Test complete UniformGrid workflow with Omega_h mesh and field + interpolation (Omega_h → UniformGrid). + """ + # Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements. + mesh = omega_h.build_box( + world, omega_h.Family.SIMPLEX, 1.0, 1.0, 0.0, 4, 4, 0, False, + ) + + grid = pcms.create_uniform_grid_from_mesh(mesh, [4, 4]) + mask_field = pcms.create_uniform_grid_binary_field(mesh, [4, 4]) + mask_data = mask_field.get_dof_holder_data() + + # Build a mesh-backed source field: f(x,y) = x + 2y. + omega_h_factory = pcms.LagrangeFunctionSpace.from_mesh(mesh, 1) + omega_h_field = omega_h_factory.create_field() + coords = omega_h_field.get_dof_holder_coordinates() + omega_h_data = np.zeros(omega_h_field.get_num_dof_holders(), + dtype=np.float64) + for i in range(omega_h_field.get_num_dof_holders()): + omega_h_data[i] = coords[i, 0] + 2.0 * coords[i, 1] + omega_h_field.set_dof_holder_data(omega_h_data) + + ug_factory = pcms.LagrangeFunctionSpace.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + ug_field = ug_factory.create_field() + + # Transfer from unstructured mesh to uniform grid. + interp = pcms.Interpolator(omega_h_factory, ug_factory) + interp.apply(omega_h_field, ug_field) + + ug_field_data = ug_field.get_dof_holder_data() + ug_coords = ug_field.get_dof_holder_coordinates() + + assert len(mask_data) == 25 + assert len(ug_field_data) == 25 + + for vertex_id in range(len(ug_field_data)): + x, y = ug_coords[vertex_id, 0], ug_coords[vertex_id, 1] + expected = x + 2.0 * y + actual = ug_field_data[vertex_id] + assert abs(actual - expected) < 1e-10 + + def test_omega_h_to_omega_h_transfer(self, world): + """Transfer a field from one Omega_h mesh to another Omega_h mesh.""" + src_mesh = omega_h.build_box( + world, omega_h.Family.SIMPLEX, 1.0, 1.0, 0.0, 4, 4, 0, False, + ) + tgt_mesh = omega_h.build_box( + world, omega_h.Family.SIMPLEX, 1.0, 1.0, 0.0, 8, 8, 0, False, + ) + + src_factory = pcms.LagrangeFunctionSpace.from_mesh(src_mesh, 1) + tgt_factory = pcms.LagrangeFunctionSpace.from_mesh(tgt_mesh, 1) + src_field = src_factory.create_field() + tgt_field = tgt_factory.create_field() + + src_coords = src_field.get_dof_holder_coordinates() + src_data = np.zeros(src_field.get_num_dof_holders()) + for i in range(len(src_data)): + src_data[i] = src_coords[i, 0] + 2.0 * src_coords[i, 1] + src_field.set_dof_holder_data(src_data) + + interp = pcms.Interpolator(src_factory, tgt_factory) + interp.apply(src_field, tgt_field) + + tgt_coords = tgt_field.get_dof_holder_coordinates() + tgt_data = tgt_field.get_dof_holder_data() + for i in range(tgt_field.get_num_dof_holders()): + expected = tgt_coords[i, 0] + 2.0 * tgt_coords[i, 1] + assert abs(expected - tgt_data[i]) < 1e-10 +