diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 000000000..f424cdd4b --- /dev/null +++ b/.devcontainer/devcontainer.json @@ -0,0 +1,37 @@ +{ + "name": "ANUGA", + "image": "mcr.microsoft.com/devcontainers/miniconda:latest", + + // onCreateCommand is cached by Codespaces prebuilds — slow conda install runs once. + "onCreateCommand": "conda env create -f .devcontainer/environment.yml", + + // postCreateCommand runs after the source tree is available — installs ANUGA itself and registers the Jupyter kernel. + "postCreateCommand": "conda run -n anuga_env pip install --no-build-isolation -e . && conda run -n anuga_env python -m ipykernel install --user --name anuga_env --display-name 'Python (anuga_env)'", + + // Activate the conda environment in every new terminal. + "postStartCommand": "conda init bash && echo 'conda activate anuga_env' >> ~/.bashrc", + + "customizations": { + "vscode": { + "extensions": [ + "ms-python.python", + "ms-python.pylance", + "ms-toolsai.jupyter" + ], + "settings": { + "python.defaultInterpreterPath": "/opt/conda/envs/anuga_env/bin/python", + "python.terminal.activateEnvironment": true, + "jupyter.kernels.filter": [ + { + "path": "/opt/conda/envs/anuga_env/bin/python", + "type": "pythonEnvironment" + } + ] + } + } + }, + + "remoteEnv": { + "CONDA_DEFAULT_ENV": "anuga_env" + } +} diff --git a/.devcontainer/environment.yml b/.devcontainer/environment.yml new file mode 100644 index 000000000..541fee9fa --- /dev/null +++ b/.devcontainer/environment.yml @@ -0,0 +1,39 @@ +name: anuga_env +channels: + - conda-forge +dependencies: + - python=3.12 + - pip + # Build tools + - compilers + - cython + - meson + - meson-python=0.17 + - ninja + - pybind11 + - pkg-config + # Runtime dependencies + - dill + - matplotlib + - meshpy + - netcdf4 + - numpy + - pyproj + - scipy + - xarray + # Optional but useful + - mpi4py + - pymetis + - rasterio + - fiona + - shapely + - pandas + - openpyxl + - cartopy + # Testing and notebooks + - pytest + - pytest-regressions + - jupyter + - ipykernel + - pip: + - pmw diff --git a/.github/workflows/conda-setup.yml b/.github/workflows/conda-setup.yml index 09824f3b8..7b9ba73ef 100644 --- a/.github/workflows/conda-setup.yml +++ b/.github/workflows/conda-setup.yml @@ -33,7 +33,7 @@ jobs: - uses: actions/checkout@v6 with: fetch-depth: 0 # fetch full history so git describe can find tags - - uses: conda-incubator/setup-miniconda@v3 + - uses: conda-incubator/setup-miniconda@v4 with: miniforge-version: latest # Installs Miniforge instead of Miniconda use-mamba: true # Recommended for faster installs diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 3bd37071f..045b64013 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -1,129 +1,184 @@ -# This workflow will upload a Python Package to PyPI when a release is created -# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python#publishing-to-package-registries +# Build wheels on every push / PR for early feedback. +# On a GitHub Release publish, build wheels then upload to PyPI. -# This workflow uses actions that are not certified by GitHub. -# They are provided by a third-party and are governed by -# separate terms of service, privacy policy, and support -# documentation. - -name: Python Publish PYPI +name: Build and Publish to PyPI on: + push: + branches: [main, develop] + pull_request: + branches: [main, develop] release: types: [published] - - permissions: contents: read jobs: - release-build: + # --------------------------------------------------------------------------- + # Build one wheel per (python-version, os) combination. + # Runs on every push/PR so broken wheels are caught before tagging. + # --------------------------------------------------------------------------- + build-wheels: + name: Build wheel (py${{ matrix.python-version }}, ${{ matrix.os }}) runs-on: ${{ matrix.os }} strategy: - fail-fast: true + fail-fast: false matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, macos-latest, windows-latest] python-version: ['3.10', '3.11', '3.12', '3.13', '3.14'] - #python-version: ['3.10', '3.11'] steps: - uses: actions/checkout@v6 with: - fetch-depth: 0 # fetch full history so git describe can find tags - - uses: conda-incubator/setup-miniconda@v3 + fetch-depth: 0 # git describe needs full history for version + + - uses: conda-incubator/setup-miniconda@v4 with: - miniforge-version: latest # Installs Miniforge instead of Miniconda - use-mamba: true # Recommended for faster installs + miniforge-version: latest + use-mamba: true python-version: ${{ matrix.python-version }} - auto-update-conda: true # Auto-update conda or mamba + auto-update-conda: true + conda-remove-defaults: true - - name: Install Our own environment + - name: Create conda environment with build deps shell: bash -el {0} run: | - conda create -n anuga_env -c conda-forge python=${{matrix.python-version}} + conda create -n anuga_env -c conda-forge \ + python=${{ matrix.python-version }} pip \ + meson-python meson ninja cython pybind11 numpy pkg-config - - name: Install gcc compilers on Windows + - name: Install compilers (Windows) if: runner.os == 'Windows' shell: bash -el {0} - run: | - conda install -c conda-forge -n anuga_env libpython gcc_win-64 gxx_win-64 + run: conda install -c conda-forge -n anuga_env libpython gcc_win-64 gxx_win-64 - # conda install -c conda-forge -n anuga_env compilers - # The compilers package on windows which uses clang. We run into - # problems when testing anuga, so we skip it for now. - - - name: Install clang with openmp compiler on macOS + - name: Install compilers (macOS) if: runner.os == 'macOS' shell: bash -el {0} - run: | - conda install -c conda-forge -n anuga_env compilers + run: conda install -c conda-forge -n anuga_env compilers - - name: Install gxx with openmp compiler on Linux + - name: Install compilers (Linux) if: runner.os == 'Linux' shell: bash -el {0} - run: | - conda install -c conda-forge -n anuga_env compilers + run: conda install -c conda-forge -n anuga_env compilers - - name: Build wheels + - name: Build wheel shell: bash -el {0} run: | conda activate anuga_env - python -m pip install build - python -m build --wheel --outdir dist-wheel + export PKG_CONFIG_PATH="$CONDA_PREFIX/lib/pkgconfig:$CONDA_PREFIX/share/pkgconfig" + pip wheel --no-build-isolation --no-deps -w dist-wheel . - - name: Repair wheels + - name: Repair wheel (bundle external libs; self-contained per platform) shell: bash -el {0} run: | conda activate anuga_env - python -m pip install repairwheel - repairwheel -o dist dist-wheel/*.whl + pip install repairwheel + # repairwheel dispatches per platform: auditwheel (Linux, manylinux), + # delocate (macOS), delvewheel (Windows). On macOS this copies + # libomp.dylib into the wheel and rewrites the @rpath so the wheel is + # self-contained — without it `import anuga` from the published wheel + # fails with "Library not loaded: @rpath/libomp.dylib" (the build links + # conda's libomp via an rpath into the conda env). delvewheel does the + # equivalent DLL bundling on Windows. Pass the conda lib dirs that + # exist on this OS so dependent libraries are found at repair time. + libdirs=() + for d in "$CONDA_PREFIX/lib" "$CONDA_PREFIX/Library/bin"; do + [ -d "$d" ] && libdirs+=(-l "$d") + done + repairwheel "${libdirs[@]}" -o dist dist-wheel/*.whl + + - name: Upload wheel artifact + uses: actions/upload-artifact@v7 + with: + name: wheel-py${{ matrix.python-version }}-${{ matrix.os }} + path: dist/*.whl + + # --------------------------------------------------------------------------- + # Build the source distribution (once, on Linux / Python 3.13). + # --------------------------------------------------------------------------- + build-sdist: + name: Build source distribution + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v6 + with: + fetch-depth: 0 - - name: Build source distribution - if: runner.os == 'linux' && matrix.python-version == '3.13' + - uses: conda-incubator/setup-miniconda@v4 + with: + miniforge-version: latest + use-mamba: true + auto-update-conda: true + conda-remove-defaults: true + + - name: Create conda environment + shell: bash -el {0} + run: | + conda create -n anuga_env -c conda-forge python=3.13 pip \ + meson-python meson ninja cython pybind11 numpy pkg-config + + - name: Build sdist shell: bash -el {0} run: | conda activate anuga_env + # PKG_CONFIG_PATH is inherited by the isolated build venv so meson + # finds conda's Python headers instead of the broken system ones. + export PKG_CONFIG_PATH="$CONDA_PREFIX/lib/pkgconfig:$CONDA_PREFIX/share/pkgconfig" + pip install build python -m build --sdist --outdir dist - - - name: Upload distributions + + - name: Upload sdist artifact uses: actions/upload-artifact@v7 with: - name: release-dists-${{ matrix.python-version }}-${{ matrix.os }} - path: dist/ - - - - - pypi-publish: - name: Upload release to PyPI + name: sdist + path: dist/*.tar.gz + + # --------------------------------------------------------------------------- + # Publish to PyPI — only runs when a GitHub Release is published. + # Requires all wheel builds and the sdist to succeed first. + # --------------------------------------------------------------------------- + publish: + name: Publish to PyPI + needs: [build-wheels, build-sdist] runs-on: ubuntu-latest - needs: - - release-build - + if: github.event_name == 'release' && github.event.action == 'published' + environment: name: anuga url: https://pypi.org/p/anuga permissions: - id-token: write # IMPORTANT: this permission is mandatory for trusted publishing + id-token: write # required for trusted publishing (OIDC) + contents: write # required to upload release assets steps: - - name: Retrieve release distributions + - name: Download all wheel artifacts uses: actions/download-artifact@v8 with: - pattern: release-dists-* + pattern: wheel-* path: dist merge-multiple: true - - run: ls ./dist # All wheels from all matrix jobs are now here + - name: Download sdist artifact + uses: actions/download-artifact@v8 + with: + name: sdist + path: dist + - name: List distributions to upload + run: ls -lh dist/ - - name: Publish release distributions to PyPI + - name: Publish to PyPI uses: pypa/gh-action-pypi-publish@release/v1 with: - #repository-url: https://test.pypi.org/legacy/ verify-metadata: false verbose: true packages-dir: dist/ + - name: Upload wheels and sdist as GitHub Release assets + env: + GH_TOKEN: ${{ github.token }} + run: gh release upload ${{ github.ref_name }} dist/* --clobber --repo ${{ github.repository }} diff --git a/.gitignore b/.gitignore index 19b5267d3..e6b3c478f 100644 --- a/.gitignore +++ b/.gitignore @@ -225,3 +225,4 @@ doc/source/generated ################################### anuga/__config__.py anuga/version.py +benchmarks/results/ diff --git a/README.rst b/README.rst index 043798491..1d616adaa 100644 --- a/README.rst +++ b/README.rst @@ -29,7 +29,11 @@ :target: https://anuga.readthedocs.io/en/latest/?badge=latest :alt: Documentation Status -|badge0| |badge1| |badge3| |badge4| |badge5| |badge6| |badge7| +.. |badge8| image:: https://github.com/codespaces/badge.svg + :target: https://codespaces.new/anuga-community/anuga_core + :alt: Open in GitHub Codespaces + +|badge0| |badge1| |badge3| |badge4| |badge5| |badge6| |badge7| |badge8| @@ -92,6 +96,17 @@ Once the conda-forge channel has been enabled, anuga can be installed with conda For more installation instructions, see https://anuga.readthedocs.io/en/latest/installation.html +GitHub Codespaces +----------------- + +Click the **Open in GitHub Codespaces** badge above to launch ANUGA in a fully configured +cloud environment — no local install required. The environment installs all dependencies +and compiles the C/Cython extensions automatically. + +Once the Codespace is ready, open any notebook and select the **Python (anuga_env)** kernel +from the kernel picker in the top-right corner. + + Documentation and Help ---------------------- diff --git a/anuga/__init__.py b/anuga/__init__.py index c8b7dde12..8d4845936 100644 --- a/anuga/__init__.py +++ b/anuga/__init__.py @@ -84,6 +84,7 @@ from anuga.geometry.polygon_function import Polygon_function from anuga.coordinate_transforms.lat_long_UTM_conversion import LLtoUTM, UTMtoLL +from anuga.coordinate_transforms.redfearn import epsg_to_ll, ll_to_epsg from anuga.abstract_2d_finite_volumes.pmesh2domain import \ pmesh_to_domain_instance, \ @@ -323,6 +324,8 @@ def get_args(): from anuga.extras import create_basic_mesh_from_regions from anuga.utilities import log as log +from anuga.utilities.log import set_logfile, TeeStream, file_only +from anuga.utilities.log import verbose as log_verbose from anuga.config import g from anuga.config import velocity_protection @@ -436,6 +439,8 @@ def get_args(): # Coordinate transforms 'LLtoUTM', 'UTMtoLL', + 'epsg_to_ll', + 'll_to_epsg', # Parallel API 'barrier', 'collect_value', @@ -491,6 +496,8 @@ def get_args(): 'memory_stats', 'print_memory_stats', 'log', + 'set_logfile', + 'TeeStream', 'parse_standard_args', 'parse_time', 'plot_utils', diff --git a/anuga/abstract_2d_finite_volumes/basic_mesh.py b/anuga/abstract_2d_finite_volumes/basic_mesh.py index b5c34632d..2b5be575f 100644 --- a/anuga/abstract_2d_finite_volumes/basic_mesh.py +++ b/anuga/abstract_2d_finite_volumes/basic_mesh.py @@ -146,6 +146,14 @@ def reorder(self, new_order, in_place=True): ------- Basic_mesh """ + # Force neighbour computation before any reordering. + # _triangle_neighbours (if set) contains pre-reorder indices; once + # the triangle numbering changes those indices become stale. Building + # _neighbours now ensures they are remapped correctly below and + # prevents _build_neighbours() from later reconstructing from the + # stale _triangle_neighbours cache. + _ = self.neighbours # triggers _build_neighbours() if not yet done + new_order = num.array(new_order, int) N = self.number_of_triangles inv_order = num.empty_like(new_order) diff --git a/anuga/abstract_2d_finite_volumes/gauge.py b/anuga/abstract_2d_finite_volumes/gauge.py index e22917461..2e9ea4075 100644 --- a/anuga/abstract_2d_finite_volumes/gauge.py +++ b/anuga/abstract_2d_finite_volumes/gauge.py @@ -608,7 +608,7 @@ def gauge_get_from_file(filename): if len(fields) > 2: elev.append(float(fields[elev_index])) loc = fields[name_index] - gaugelocation.append(loc.strip(r'\n')) + gaugelocation.append(loc.strip()) return gauges, gaugelocation, elev diff --git a/anuga/abstract_2d_finite_volumes/generic_domain.py b/anuga/abstract_2d_finite_volumes/generic_domain.py index d6d577feb..52128977a 100644 --- a/anuga/abstract_2d_finite_volumes/generic_domain.py +++ b/anuga/abstract_2d_finite_volumes/generic_domain.py @@ -576,6 +576,14 @@ def set_georeference(self, *args, **kwargs): self.mesh.set_georeference(*args, **kwargs) self.geo_reference = self.mesh.geo_reference + def get_epsg(self): + """Return the EPSG code of the domain's coordinate reference system.""" + return self.geo_reference.get_epsg() + + def set_epsg(self, epsg): + """Set the EPSG code of the domain's coordinate reference system.""" + self.geo_reference.epsg = int(epsg) + def build_boundary_dictionary(self, *args, **kwargs): self.mesh.build_boundary_dictionary(*args, **kwargs) diff --git a/anuga/abstract_2d_finite_volumes/quantity.py b/anuga/abstract_2d_finite_volumes/quantity.py index d58899b42..102d91db8 100644 --- a/anuga/abstract_2d_finite_volumes/quantity.py +++ b/anuga/abstract_2d_finite_volumes/quantity.py @@ -1731,41 +1731,15 @@ def set_values_from_lat_long_grid_file(self, if verbose: print(self.domain.geo_reference) - utm_zone = self.domain.geo_reference.get_zone() - utm_hemisphere = self.domain.geo_reference.get_hemisphere() - - northern = True - if utm_hemisphere == 'southern': - northern = False - - #import re - #utm_zone_number = re.findall(r'\d+', utm_zone)[0] - #utm_zone_letter = re.findall(r'[A-z]+', utm_zone)[0] - - #print(utm_zone) - #print(points) - - # we could use anuga's utmtoLL but it has not been vectorised so lets - # use this library, but we will have to download via pip - import utm - lat, long = utm.to_latlon(points[:,0], points[:,1], utm_zone, northern=northern) - - #print(lat) - #print(long) + from anuga.coordinate_transforms.redfearn import epsg_to_ll + epsg = self.domain.geo_reference.get_epsg() + lat, long = epsg_to_ll(points[:,0], points[:,1], epsg) lat = num.reshape(lat, (-1,1)) long = num.reshape(long, (-1,1)) points_ll = num.hstack((long,lat)) - # need to pull out the the utm zone number and letter - - - - #import utm - #points_ll = utm.to_latlon(easting = points[:,0], northing=points[:,1]) - - #print('points_ll', points_ll) from anuga.fit_interpolate.interpolate2d import interpolate2d diff --git a/anuga/abstract_2d_finite_volumes/tests/test_gauge.py b/anuga/abstract_2d_finite_volumes/tests/test_gauge.py index 54782ceda..060d0e0f4 100644 --- a/anuga/abstract_2d_finite_volumes/tests/test_gauge.py +++ b/anuga/abstract_2d_finite_volumes/tests/test_gauge.py @@ -486,28 +486,36 @@ def test_sww2csv_output_centroid_attribute(self): def test_sww2csv_multiple_files(self): """ - This is testing the sww2csv_gauges function, by creating multiple + This is testing the sww2csv_gauges function, by creating multiple sww files and then exporting the gauges and checking the results. """ + import shutil + tmpdir = tempfile.mkdtemp() + orig_dir = os.getcwd() + os.chdir(tmpdir) + try: + self._test_sww2csv_multiple_files_impl() + finally: + self.sww = None # files are inside tmpdir; tearDown must not touch them + os.chdir(orig_dir) + shutil.rmtree(tmpdir, ignore_errors=True) + + def _test_sww2csv_multiple_files_impl(self): timestep=2.0 domain = self.domain domain.set_starttime(0.) + domain.set_datadir('.') + # Create two sww files with timestep at end. These are to be # stored consecutively in the gauge csv files basename='datatest1' - domain.set_name(basename) + domain.set_name(basename) self._create_sww(stage=10.,timestep=timestep) - domain.set_name(basename+str(time.time())) + domain.set_name(basename+str(time.time())) domain.set_time(domain.get_time()+timestep) self._create_sww(stage=20.,timestep=timestep) - #points_file = tempfile.mktemp(".csv") - #file_id = open(points_file,"w") - - # test the function at these points - points = [[5.0,1.],[0.5,2.]] - # create a csv file containing our gauge points points_file = tempfile.mktemp(".csv") points_handle = open(points_file,"w") @@ -516,16 +524,12 @@ def test_sww2csv_multiple_files(self): point2, 0.5, 2.0\n") points_handle.close() - - sww2csv_gauges(basename+".sww", + sww2csv_gauges(basename+".sww", points_file, quantities=['stage', 'elevation'], use_cache=False, verbose=False) - point1_answers_array = [[0.0,1.0,-5.0], [2.0,10.0,-5.0],[4.0,10.0,-5.0], - [6.0,20.0,-5.0], [0.0,1.0,-5.0]] - point1_answers_array = [[0.0, 1.0, -3.0], [2.0, 10.0, -3.0], [4.0, 10.0, -3.0], [6.0, 20.0, -3.0]] @@ -538,46 +542,26 @@ def test_sww2csv_multiple_files(self): for i,row in enumerate(point1_reader): # note the 'hole' (element 1) below - skip the new 'hours' field line.append([float(row[0]),float(row[2]),float(row[3])]) - #print 'i', i - #print 'row',row - #print 'line',line[i],'point1',point1_answers_array[i] assert num.allclose(line[i], point1_answers_array[i]) + point1_handle.close() - #point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5],[4.0,10.0,-0.5], - # [6.0,20.0,-0.5], [0.0,1.0,-0.5]] point2_answers_array = [[0.0, 1.0, -2.416666666666667], [2.0, 10.000000000000002, -2.416666666666667], [4.0, 10.000000000000002, -2.416666666666667], [6.0, 20.000000000000004, -2.416666666666667]] - - - - point2_filename = 'gauge_point2.csv' + point2_filename = 'gauge_point2.csv' point2_handle = open(point2_filename) point2_reader = reader(point2_handle) next(point2_reader) - + line=[] for i,row in enumerate(point2_reader): # note the 'hole' (element 1) below - skip the new 'hours' field line.append([float(row[0]),float(row[2]),float(row[3])]) - #print 'line',line[i],'point2'#,point2_answers_array[i] assert num.allclose(line[i], point2_answers_array[i]) - - # clean up - point1_handle.close() point2_handle.close() - try: - os.remove(points_file) - os.remove(point1_filename) - os.remove(point2_filename) - #remove second swwfile not removed by tearDown - os.remove(basename+".sww") - except OSError: - pass - #------------------------------------------------------------- diff --git a/anuga/abstract_2d_finite_volumes/tests/test_neighbour_mesh_reorder.py b/anuga/abstract_2d_finite_volumes/tests/test_neighbour_mesh_reorder.py index a969a3da3..b6f6dda4e 100644 --- a/anuga/abstract_2d_finite_volumes/tests/test_neighbour_mesh_reorder.py +++ b/anuga/abstract_2d_finite_volumes/tests/test_neighbour_mesh_reorder.py @@ -198,6 +198,54 @@ def test_reorder_larger(self): compare_meshes(new_mesh, reorder_mesh) + def test_basic_mesh_reorder_neighbours_consistent(self): + """Basic_mesh.reorder() must produce correct neighbours even when + _neighbours has not been accessed before the reorder call. + + Regression test for a bug where the lazy-built _neighbours were + constructed from an unreordered _triangle_neighbours cache after + reorder(), causing ghost-layer BFS to follow wrong adjacency and + produce more ghost triangles than distribute() for the same mesh. + """ + import numpy as np + from anuga.abstract_2d_finite_volumes.basic_mesh import ( + rectangular_cross_basic_mesh) + + # Build a small mesh. _neighbours starts as None (not yet accessed). + bm = rectangular_cross_basic_mesh(5, 5, len1=5.0, len2=5.0) + assert bm._neighbours is None, \ + "Precondition: _neighbours should be None before first access" + + # Capture ground-truth neighbours BEFORE reordering. + nbrs_before = bm.neighbours.copy() # triggers _build_neighbours + + # Reset so we can test the lazy path. + bm2 = rectangular_cross_basic_mesh(5, 5, len1=5.0, len2=5.0) + assert bm2._neighbours is None + + # Apply a non-trivial permutation. + N = bm2.number_of_triangles + rng = np.random.default_rng(42) + new_order = rng.permutation(N) + + reordered = bm2.reorder(new_order, in_place=False) + + # After reorder the neighbours must be self-consistent: + # for each triangle i and each edge j, if reordered.neighbours[i,j] == k + # then triangle k must be adjacent to triangle i. + nbrs = reordered.neighbours + tris = reordered.triangles + for i in range(N): + for j in range(3): + k = nbrs[i, j] + if k < 0: + continue # boundary edge, skip + # Triangles i and k must share exactly 2 nodes. + shared = set(tris[i]) & set(tris[k]) + assert len(shared) == 2, ( + f"Triangle {i} and neighbour {k} share {len(shared)} " + f"nodes (expected 2) after reorder — neighbours are stale") + # def test_reorder_larger_16_16(self): # """Test larger mesh which failed in sequential_dist example diff --git a/anuga/abstract_2d_finite_volumes/tests/test_quantity.py b/anuga/abstract_2d_finite_volumes/tests/test_quantity.py index b8d20b970..5389f3c21 100644 --- a/anuga/abstract_2d_finite_volumes/tests/test_quantity.py +++ b/anuga/abstract_2d_finite_volumes/tests/test_quantity.py @@ -2308,12 +2308,12 @@ def test_set_values_from_ll_tif_file_north(self): lat_ll, lon_ll = 34.37, 150.90 lat_ur, lon_ur = 34.39, 150.92 - import utm - utm_east_ll, utm_north_ll, zone_ll, zone_letter_ll = utm.from_latlon(lat_ll, lon_ll) - utm_east_ur, utm_north_ur, zone_ur, zone_letter_ur = utm.from_latlon(lat_ur, lon_ur) + from anuga.coordinate_transforms.lat_long_UTM_conversion import LLtoUTM + zone_ll, utm_east_ll, utm_north_ll = LLtoUTM(lat_ll, lon_ll) + zone_ur, utm_east_ur, utm_north_ur = LLtoUTM(lat_ur, lon_ur) - hemisphere_ll = zone_letter_to_hemisphere(zone_letter_ll) - hemisphere_ur = zone_letter_to_hemisphere(zone_letter_ur) + hemisphere_ll = 'southern' if lat_ll < 0 else 'northern' + hemisphere_ur = 'southern' if lat_ur < 0 else 'northern' lat = numpy.linspace(lat_ll, lat_ur, 11) lon = numpy.linspace(lon_ll, lon_ur, 11) @@ -2417,12 +2417,12 @@ def test_set_values_from_ll_tif_file_north_indices(self): lat_ll, lon_ll = 34.37, 150.90 lat_ur, lon_ur = 34.39, 150.92 - import utm - utm_east_ll, utm_north_ll, zone_ll, zone_letter_ll = utm.from_latlon(lat_ll, lon_ll) - utm_east_ur, utm_north_ur, zone_ur, zone_letter_ur = utm.from_latlon(lat_ur, lon_ur) + from anuga.coordinate_transforms.lat_long_UTM_conversion import LLtoUTM + zone_ll, utm_east_ll, utm_north_ll = LLtoUTM(lat_ll, lon_ll) + zone_ur, utm_east_ur, utm_north_ur = LLtoUTM(lat_ur, lon_ur) - hemisphere_ll = zone_letter_to_hemisphere(zone_letter_ll) - hemisphere_ur = zone_letter_to_hemisphere(zone_letter_ur) + hemisphere_ll = 'southern' if lat_ll < 0 else 'northern' + hemisphere_ur = 'southern' if lat_ur < 0 else 'northern' lat = numpy.linspace(lat_ll, lat_ur, 11) lon = numpy.linspace(lon_ll, lon_ur, 11) @@ -2526,12 +2526,12 @@ def test_set_values_from_ll_tif_file_south(self): lat_ll, lon_ll = -34.39, 150.90 lat_ur, lon_ur = -34.37, 150.92 - import utm - utm_east_ll, utm_north_ll, zone_ll, zone_letter_ll = utm.from_latlon(lat_ll, lon_ll) - utm_east_ur, utm_north_ur, zone_ur, zone_letter_ur = utm.from_latlon(lat_ur, lon_ur) + from anuga.coordinate_transforms.lat_long_UTM_conversion import LLtoUTM + zone_ll, utm_east_ll, utm_north_ll = LLtoUTM(lat_ll, lon_ll) + zone_ur, utm_east_ur, utm_north_ur = LLtoUTM(lat_ur, lon_ur) - hemisphere_ll = zone_letter_to_hemisphere(zone_letter_ll) - hemisphere_ur = zone_letter_to_hemisphere(zone_letter_ur) + hemisphere_ll = 'southern' if lat_ll < 0 else 'northern' + hemisphere_ur = 'southern' if lat_ur < 0 else 'northern' lat = numpy.linspace(lat_ll, lat_ur, 11) lon = numpy.linspace(lon_ll, lon_ur, 11) @@ -2625,12 +2625,12 @@ def test_set_values_from_utm_tif_file(self): lat_ll, lon_ll = -34.39, 150.90 lat_ur, lon_ur = -34.37, 150.92 - import utm - utm_east_ll, utm_north_ll, zone_ll, zone_letter_ll = utm.from_latlon(lat_ll, lon_ll) - utm_east_ur, utm_north_ur, zone_ur, zone_letter_ur = utm.from_latlon(lat_ur, lon_ur) + from anuga.coordinate_transforms.lat_long_UTM_conversion import LLtoUTM + zone_ll, utm_east_ll, utm_north_ll = LLtoUTM(lat_ll, lon_ll) + zone_ur, utm_east_ur, utm_north_ur = LLtoUTM(lat_ur, lon_ur) - hemisphere_ll = zone_letter_to_hemisphere(zone_letter_ll) - hemisphere_ur = zone_letter_to_hemisphere(zone_letter_ur) + hemisphere_ll = 'southern' if lat_ll < 0 else 'northern' + hemisphere_ur = 'southern' if lat_ur < 0 else 'northern' easting = numpy.linspace(utm_east_ll, utm_east_ur, 11) northing = numpy.linspace(utm_north_ll, utm_north_ur, 11) diff --git a/anuga/abstract_2d_finite_volumes/tests/test_util.py b/anuga/abstract_2d_finite_volumes/tests/test_util.py index 5240d7f2d..2fb03c5e8 100644 --- a/anuga/abstract_2d_finite_volumes/tests/test_util.py +++ b/anuga/abstract_2d_finite_volumes/tests/test_util.py @@ -1598,7 +1598,7 @@ def test_calc_bearings(self): if 314 < angle < 316: v=1 assert v==1 - def test_calc_bearings_zero_vector(self): + def test_calc_bearings_zero_vector(self): from math import atan, degrees uh = 0 @@ -1606,7 +1606,29 @@ def test_calc_bearings_zero_vector(self): angle = calc_bearing(uh, vh) assert angle == NAN - + + def test_get_gauges_from_file(self): + """get_gauges_from_file must parse a gauge file. + + Regression test: the helper called ``gauge_get_from_file`` without + importing it (that function lives in ``gauge``), so every call raised + ``NameError``. It should return the (locations, names, elevations) + triple read from the file. + """ + handle, filename = tempfile.mkstemp('.txt') + os.close(handle) + with open(filename, 'w') as fid: + fid.write('easting, northing, name, elevation\n') + fid.write('308500, 6193000, gauge_a, 1.0\n') + fid.write('308700, 6193200, gauge_b, 2.5\n') + + gauges, gaugelocation, elev = get_gauges_from_file(filename) + os.remove(filename) + + assert gauges == [[308500.0, 6193000.0], [308700.0, 6193200.0]] + assert gaugelocation == ['gauge_a', 'gauge_b'] + assert elev == [1.0, 2.5] + #------------------------------------------------------------- if __name__ == "__main__": diff --git a/anuga/abstract_2d_finite_volumes/util.py b/anuga/abstract_2d_finite_volumes/util.py index fd5084a3a..4ee6b2995 100644 --- a/anuga/abstract_2d_finite_volumes/util.py +++ b/anuga/abstract_2d_finite_volumes/util.py @@ -131,6 +131,8 @@ def get_textual_float(value, format = '%.2f'): return format % float(value) def get_gauges_from_file(filename): + # Imported locally to avoid a circular import: gauge imports from util. + from anuga.abstract_2d_finite_volumes.gauge import gauge_get_from_file return gauge_get_from_file(filename) diff --git a/anuga/coordinate_transforms/lat_long_UTM_conversion.py b/anuga/coordinate_transforms/lat_long_UTM_conversion.py index 6a9c2eb81..34f76f0ba 100644 --- a/anuga/coordinate_transforms/lat_long_UTM_conversion.py +++ b/anuga/coordinate_transforms/lat_long_UTM_conversion.py @@ -57,77 +57,39 @@ #Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System #1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency -#def LLtoUTM(int ReferenceEllipsoid, const double Lat, const double Long, -# double &UTMNorthing, double &UTMEasting, char* UTMZone) - -def LLtoUTM( Lat, Long, ReferenceEllipsoid=23): - """ - converts lat/long to UTM coords. Equations from USGS Bulletin 1532 - East Longitudes are positive, West longitudes are negative. - North latitudes are positive, South latitudes are negative - Lat and Long are in decimal degrees - Written by Chuck Gantz- chuck.gantz@globalstar.com +def LLtoUTM(Lat, Long, ReferenceEllipsoid=23): + """Convert latitude/longitude to UTM coordinates (WGS84). + + Parameters + ---------- + Lat : float + Latitude in decimal degrees (positive north). + Long : float + Longitude in decimal degrees (positive east). + ReferenceEllipsoid : int, optional + Ignored — retained for backward compatibility. pyproj always uses WGS84. + + Returns + ------- + ZoneNumber : int + UTMEasting : float + UTMNorthing : float """ - a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius] - eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared] - k0 = 0.9996 - - #Make sure the longitude is between -180.00 .. 179.9 - LongTemp = (Long+180)-int((Long+180)/360)*360-180 # -180.00 .. 179.9 - - LatRad = Lat*_deg2rad - LongRad = LongTemp*_deg2rad - - ZoneNumber = int((LongTemp + 180)/6) + 1 - - if Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0: + # Compute UTM zone number, retaining special-zone handling for Norway/Svalbard. + LongTemp = (Long + 180) - int((Long + 180) / 360) * 360 - 180 + ZoneNumber = int((LongTemp + 180) / 6) + 1 + if 56.0 <= Lat < 64.0 and 3.0 <= LongTemp < 12.0: ZoneNumber = 32 + if 72.0 <= Lat < 84.0: + if LongTemp < 9.0: ZoneNumber = 31 + elif LongTemp < 21.0: ZoneNumber = 33 + elif LongTemp < 33.0: ZoneNumber = 35 + elif LongTemp < 42.0: ZoneNumber = 37 - # Special zones for Svalbard - if Lat >= 72.0 and Lat < 84.0: - if LongTemp >= 0.0 and LongTemp < 9.0:ZoneNumber = 31 - elif LongTemp >= 9.0 and LongTemp < 21.0: ZoneNumber = 33 - elif LongTemp >= 21.0 and LongTemp < 33.0: ZoneNumber = 35 - elif LongTemp >= 33.0 and LongTemp < 42.0: ZoneNumber = 37 - - LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 #+3 puts origin in middle of zone - LongOriginRad = LongOrigin * _deg2rad - - #compute the UTM Zone from the latitude and longitude - UTMZone = "%d%c" % (ZoneNumber, _UTMLetterDesignator(Lat)) - - eccPrimeSquared = (eccSquared)/(1-eccSquared) - N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad)) - T = tan(LatRad)*tan(LatRad) - C = eccPrimeSquared*cos(LatRad)*cos(LatRad) - A = cos(LatRad)*(LongRad-LongOriginRad) - - M = a*((1 - - eccSquared/4 - - 3*eccSquared*eccSquared/64 - - 5*eccSquared*eccSquared*eccSquared/256)*LatRad - - (3*eccSquared/8 - + 3*eccSquared*eccSquared/32 - + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad) - + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad) - - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad)) - - UTMEasting = (k0*N*(A+(1-T+C)*A*A*A/6 - + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120) - + 500000.0) - - UTMNorthing = (k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 - + (61 - -58*T - +T*T - +600*C - -330*eccPrimeSquared)*A*A*A*A*A*A/720))) - - if Lat < 0: - UTMNorthing = UTMNorthing + 10000000.0; #10000000 meter offset for southern hemisphere - #UTMZone was originally returned here. I don't know what the - #letter at the end was for. - return (ZoneNumber, UTMEasting, UTMNorthing) + from anuga.coordinate_transforms.redfearn import ll_to_epsg + epsg = 32700 + ZoneNumber if Lat < 0 else 32600 + ZoneNumber + easting, northing = ll_to_epsg(Lat, Long, epsg) + return (ZoneNumber, float(easting), float(northing)) def _UTMLetterDesignator(Lat): @@ -157,92 +119,32 @@ def _UTMLetterDesignator(Lat): elif -72 > Lat >= -80: return 'C' else: return 'Z' # if the Latitude is outside the UTM limits -#void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone, -# double& Lat, double& Long ) - def UTMtoLL(northing, easting, zone, isSouthernHemisphere=True, ReferenceEllipsoid=23): + """Convert UTM coordinates to latitude/longitude (WGS84). + + Parameters + ---------- + northing : float + UTM northing in metres. + easting : float + UTM easting in metres. + zone : int + UTM zone number (1–60). + isSouthernHemisphere : bool, optional + True (default) for southern hemisphere, False for northern. + ReferenceEllipsoid : int, optional + Ignored — retained for backward compatibility. pyproj always uses WGS84. + + Returns + ------- + lat : float + lon : float """ - converts UTM coords to lat/long. Equations from USGS Bulletin 1532 - East Longitudes are positive, West longitudes are negative. - North latitudes are positive, South latitudes are negative - Lat and Long are in decimal degrees. - Written by Chuck Gantz- chuck.gantz@globalstar.com - Converted to Python by Russ Nelson - - FIXME: This is set up to work for the Southern Hemisphere. - -Using -http://www.ga.gov.au/geodesy/datums/redfearn_geo_to_grid.jsp - - Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) -Zone: 36 -Easting: 511669.521 Northing: 19328195.112 -Latitude: 84 0 ' 0.00000 '' Longitude: 34 0 ' 0.00000 '' -Grid Convergence: 0 -59 ' 40.28 '' Point Scale: 0.99960166 - -____________ -Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) -Zone: 36 -Easting: 519384.803 Northing: 1118247.585 -Latitude: -80 0 ' 0.00000 '' Longitude: 34 0 ' 0.00000 '' -Grid Convergence: 0 59 ' 5.32 '' Point Scale: 0.99960459 -___________ -Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) -Zone: 36 -Easting: 611263.812 Northing: 10110547.106 -Latitude: 1 0 ' 0.00000 '' Longitude: 34 0 ' 0.00000 '' -Grid Convergence: 0 -1 ' 2.84 '' Point Scale: 0.99975325 -______________ -Site Name: GDA-MGA: (UTM with GRS80 ellipsoid) -Zone: 36 -Easting: 611263.812 Northing: 9889452.894 -Latitude: -1 0 ' 0.00000 '' Longitude: 34 0 ' 0.00000 '' -Grid Convergence: 0 1 ' 2.84 '' Point Scale: 0.99975325 - -So this uses a false northing of 10000000 in the both hemispheres. -ArcGIS used a false northing of 0 in the northern hem though. -Therefore it is difficult to actually know what hemisphere you are in. - """ - k0 = 0.9996 - a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius] - eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared] - e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared)) - - x = easting - 500000.0 #remove 500,000 meter offset for longitude - y = northing - - ZoneNumber = int(zone) - if isSouthernHemisphere: - y -= 10000000.0 # remove 10,000,000 meter offset used - # for southern hemisphere - - LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 # +3 puts origin in middle of zone - - eccPrimeSquared = (eccSquared)/(1-eccSquared) - - M = y/ k0 - mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256)) - - phi1Rad = (mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu) - + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu) - +(151*e1*e1*e1/96)*sin(6*mu)) - phi1 = phi1Rad*_rad2deg; - - N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad)) - T1 = tan(phi1Rad)*tan(phi1Rad) - C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad) - R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5) - D = x/(N1*k0) - - Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24 - +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720) - Lat = Lat * _rad2deg - - Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1) - *D*D*D*D*D/120.0)/cos(phi1Rad) - Long = LongOrigin + Long * _rad2deg - return (Lat, Long) + from anuga.coordinate_transforms.redfearn import epsg_to_ll + epsg = 32700 + int(zone) if isSouthernHemisphere else 32600 + int(zone) + lat, lon = epsg_to_ll(easting, northing, epsg) + return float(lat), float(lon) if __name__ == '__main__': (z, e, n) = LLtoUTM(-45.00, -75.00, 23) diff --git a/anuga/coordinate_transforms/redfearn.py b/anuga/coordinate_transforms/redfearn.py index 3f5f07be5..76cf66cf3 100644 --- a/anuga/coordinate_transforms/redfearn.py +++ b/anuga/coordinate_transforms/redfearn.py @@ -10,6 +10,66 @@ import numpy as num + +def epsg_to_ll(easting, northing, epsg): + """Convert projected coordinates to latitude/longitude using pyproj. + + Parameters + ---------- + easting : array-like + Easting coordinates in the CRS defined by epsg. + northing : array-like + Northing coordinates in the CRS defined by epsg. + epsg : int + EPSG code of the source CRS (e.g. 32755 for UTM zone 55S). + + Returns + ------- + lat : ndarray + lon : ndarray + """ + from pyproj import Transformer + t = Transformer.from_crs(epsg, 4326, always_xy=True) + e = num.asarray(easting, dtype=num.float64) + n = num.asarray(northing, dtype=num.float64) + # pyproj's fast scalar path calls float() on its inputs internally. + # In NumPy >= 2.0, float(array_with_ndim > 0) is deprecated. + # For single-element inputs use .item() to give pyproj a plain Python + # float; for multi-element arrays the array path is taken automatically. + if e.size == 1: + lon, lat = t.transform(e.item(), n.item()) + return num.asarray(lat).reshape(e.shape), num.asarray(lon).reshape(e.shape) + lon, lat = t.transform(e, n) + return num.asarray(lat), num.asarray(lon) + + +def ll_to_epsg(lat, lon, epsg): + """Convert latitude/longitude to projected coordinates using pyproj. + + Parameters + ---------- + lat : array-like + Latitudes in decimal degrees. + lon : array-like + Longitudes in decimal degrees. + epsg : int + EPSG code of the target CRS (e.g. 32755 for UTM zone 55S). + + Returns + ------- + easting : ndarray + northing : ndarray + """ + from pyproj import Transformer + t = Transformer.from_crs(4326, epsg, always_xy=True) + la = num.asarray(lat, dtype=num.float64) + lo = num.asarray(lon, dtype=num.float64) + if la.size == 1: + easting, northing = t.transform(lo.item(), la.item()) + return num.asarray(easting).reshape(la.shape), num.asarray(northing).reshape(la.shape) + easting, northing = t.transform(lo, la) + return num.asarray(easting), num.asarray(northing) + def degminsec2decimal_degrees(dd,mm,ss): assert abs(mm) == mm assert abs(ss) == ss diff --git a/anuga/file_conversion/tif2point_values.py b/anuga/file_conversion/tif2point_values.py index be8526dd8..27c3b0381 100644 --- a/anuga/file_conversion/tif2point_values.py +++ b/anuga/file_conversion/tif2point_values.py @@ -11,7 +11,6 @@ def tif2point_values(filename, zone=None, south=True, points=None, verbose=False import numpy as np import rasterio from pyproj import CRS, Transformer - from affine import Affine with rasterio.open(filename) as raster: ncols = raster.width @@ -45,7 +44,13 @@ def tif2point_values(filename, zone=None, south=True, points=None, verbose=False tif_georeference = CRS.from_epsg(4326) transformer = Transformer.from_crs(points_utm, tif_georeference) - points_lat, points_lon = transformer.transform(points[:, 0], points[:, 1]) + # pyproj dispatches to _transform_point (scalar path) when given a + # 1-element array, which fails with numpy >= 1.25. Pass plain Python + # lists so pyproj always uses the array path, then convert back. + _lat, _lon = transformer.transform( + points[:, 0].tolist(), points[:, 1].tolist()) + points_lat = np.asarray(_lat) + points_lon = np.asarray(_lon) ilocs = np.array(~affine_transform * (points_lon, points_lat)) diff --git a/anuga/geospatial_data/geospatial_data.py b/anuga/geospatial_data/geospatial_data.py index 197983f0e..c63bd22d5 100644 --- a/anuga/geospatial_data/geospatial_data.py +++ b/anuga/geospatial_data/geospatial_data.py @@ -17,7 +17,6 @@ import numpy as num from numpy.random import randint, seed -from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL from anuga.utilities.numerical_tools import ensure_numeric from anuga.coordinate_transforms.geo_reference import Geo_reference, \ TitleError, DEFAULT_ZONE, ensure_geo_reference, write_NetCDF_georeference @@ -352,15 +351,15 @@ def get_data_points(self, if as_lat_long is True: msg = "Points need a zone to be converted into lats and longs" assert self.geo_reference is not None, msg - zone = self.geo_reference.get_zone() assert self.geo_reference.get_zone() is not DEFAULT_ZONE, msg - lats_longs = [] - for point in self.get_data_points(True): - # UTMtoLL(northing, easting, zone, - lat_calced, long_calced = UTMtoLL(point[1], point[0], - zone, isSouthHemisphere) - lats_longs.append((lat_calced, long_calced)) # to hash - return lats_longs + from anuga.coordinate_transforms.redfearn import epsg_to_ll + pts = self.get_data_points(True) + epsg = self.geo_reference.get_epsg() + if epsg is None: + zone = self.geo_reference.get_zone() + epsg = 32700 + zone if isSouthHemisphere else 32600 + zone + lats, lons = epsg_to_ll(pts[:, 0], pts[:, 1], epsg) + return list(zip(lats.tolist(), lons.tolist())) if absolute is True and geo_reference is None: return self.geo_reference.get_absolute(self.data_points) diff --git a/anuga/scenario/prepare_data.py b/anuga/scenario/prepare_data.py index c4af586ae..2d6ea7109 100644 --- a/anuga/scenario/prepare_data.py +++ b/anuga/scenario/prepare_data.py @@ -17,6 +17,7 @@ import numpy import anuga from anuga.utilities import spatialInputUtil as su +from anuga.utilities.log import TeeStream from anuga.parallel import myid, barrier, send, receive, numprocs # Local modules @@ -25,23 +26,6 @@ from .parse_input_data import ProjectData -class Logger(object): - - """Makes it simple to get stdout to flush, and seems to work in - parallel - """ - - def __init__(self, logfile): - self.log = open(logfile, "a") - - def write(self, message): - self.log.write(message) - self.log.flush() - - def flush(self): - self.log.flush() - - class PrepareData(ProjectData): """Converts the information in the configuration file to inputs @@ -104,7 +88,7 @@ def define_output_directory_and_redirect_stdout(self, print('OUTPUT_DIRECTORY: ' + str(self.output_dir)) - # Send stdout to a file inside the output directory + # Tee stdout to a file inside the output directory if output_log is not None: if make_directories: stdout_file = self.output_dir + '/' + output_log @@ -113,11 +97,7 @@ def define_output_directory_and_redirect_stdout(self, if myid == 0: print('Redirecting output now to ' + stdout_file) - sys.stdout = Logger(stdout_file) - barrier() - - if myid != 0: - sys.stdout = Logger(stdout_file) + sys.stdout = TeeStream(stdout_file) return diff --git a/anuga/scenario/setup_mesh.py b/anuga/scenario/setup_mesh.py index a234390ce..51229260e 100644 --- a/anuga/scenario/setup_mesh.py +++ b/anuga/scenario/setup_mesh.py @@ -18,6 +18,7 @@ import anuga from anuga.parallel import myid, numprocs, barrier from anuga.utilities import spatialInputUtil as su +import anuga.utilities.log as log from anuga.parallel.parallel_api import pypar_available if pypar_available: from anuga import sequential_distribute_load @@ -25,7 +26,6 @@ from anuga.parallel.sequential_distribute import \ sequential_distribute_load_pickle_file -verbose = True def build_mesh(project): """ @@ -40,44 +40,43 @@ def build_mesh(project): su.combine_breakLines_and_riverWalls_for_mesh(project.breaklines, project.riverwalls) - # Make the mesh - - anuga.create_pmesh_from_regions( - project.bounding_polygon, - boundary_tags=project.boundary_tags, - maximum_triangle_area=project.default_res, - filename=project.meshname, - interior_regions=project.interior_regions, - use_cache=False, - verbose=verbose, - breaklines=mesh_breaklines, - regionPtArea=project.region_point_areas, - ) + # Make the mesh — verbose output captured to file only + with log.file_only(): + anuga.create_pmesh_from_regions( + project.bounding_polygon, + boundary_tags=project.boundary_tags, + maximum_triangle_area=project.default_res, + filename=project.meshname, + interior_regions=project.interior_regions, + use_cache=False, + verbose=True, + breaklines=mesh_breaklines, + regionPtArea=project.region_point_areas, + ) # Make the domain using the mesh domain = anuga.create_domain_from_file(project.meshname) - # Print some stats about mesh and domain - - print('Number of triangles = ', len(domain)) - print('The extent is ', domain.get_extent()) - print(domain.statistics()) + # Key mesh stats go to screen (info level) + log.info('Number of triangles = %d' % len(domain)) + log.info('The extent is %s' % str(domain.get_extent())) - # Print info on the smallest triangles + # Detailed stats go to file only (verbose level) + log.verbose(domain.statistics()) small_areas = domain.areas.argsort() - print('') - print('LOCATIONS OF TRIANGLES WITH SMALLEST AREAS') + log.verbose('') + log.verbose('LOCATIONS OF TRIANGLES WITH SMALLEST AREAS') for i in range(10): j = small_areas[i] x = domain.centroid_coordinates[j, 0] \ + domain.geo_reference.xllcorner y = domain.centroid_coordinates[j, 1] \ + domain.geo_reference.yllcorner - print(' Area ' + str(domain.areas[j]) + ' location: ' \ - + str(round(x, 1)) + ',' + str(round(y, 1))) - print('') + log.verbose(' Area %s location: %s,%s' + % (domain.areas[j], round(x, 1), round(y, 1))) + log.verbose('') return domain @@ -112,64 +111,59 @@ def setup_mesh(project, setup_initial_conditions=None): else: if myid == 0: - if verbose: - print('Hello from processor ', myid) + log.verbose('Hello from processor %d' % myid) pickle_name = 'domain' + '_P%g_%g.pickle' % (1, 0) pickle_name = join(project.partition_dir, pickle_name) if os.path.exists(pickle_name): - if verbose: - print('Saved domain seems to already exist') + log.verbose('Saved domain seems to already exist') else: - if verbose: - print('CREATING PARTITIONED DOMAIN') + log.info('Creating partitioned domain') domain = build_mesh(project) if setup_initial_conditions is not None: setup_initial_conditions.setup_initial_conditions( domain, project) - if verbose: - print('Saving Domain') - sequential_distribute_dump(domain, 1, - partition_dir=project.partition_dir, - verbose=verbose) + log.verbose('Saving domain') + with log.file_only(): + sequential_distribute_dump(domain, 1, + partition_dir=project.partition_dir, + verbose=True) par_pickle_name = 'domain' + '_P%g_%g.pickle' % (numprocs, 0) par_pickle_name = join(project.partition_dir, par_pickle_name) if os.path.exists(par_pickle_name): - if verbose: - print('Saved partitioned domain seems to already exist') + log.verbose('Saved partitioned domain seems to already exist') else: - if verbose: - print('Load in saved sequential pickled domain') - domain = sequential_distribute_load_pickle_file( - pickle_name, np=1, verbose=verbose) + log.verbose('Load in saved sequential pickled domain') + with log.file_only(): + domain = sequential_distribute_load_pickle_file( + pickle_name, np=1, verbose=True) - if verbose: - print('Dump partitioned domains') - sequential_distribute_dump( - domain, numprocs, - partition_dir=project.partition_dir, verbose=verbose) + log.verbose('Dump partitioned domains') + with log.file_only(): + sequential_distribute_dump( + domain, numprocs, + partition_dir=project.partition_dir, verbose=True) domain = None gc.collect() else: domain = None - if verbose: - print('Hello from processor ', myid) + log.verbose('Hello from processor %d' % myid) barrier() - if myid == 0: - print('LOADING PARTITIONED DOMAIN') + log.info('Loading partitioned domain') - domain = sequential_distribute_load( - filename=join(project.partition_dir, 'domain'), - verbose=verbose) + with log.file_only(): + domain = sequential_distribute_load( + filename=join(project.partition_dir, 'domain'), + verbose=True) # ######################################################################### # Set output directories diff --git a/anuga/utilities/log.py b/anuga/utilities/log.py index dc75e18df..7aea8bf0c 100644 --- a/anuga/utilities/log.py +++ b/anuga/utilities/log.py @@ -1,83 +1,200 @@ #!/usr/bin/env python """ -A simple logging module that logs to the console and a logfile, and has a -configurable threshold loglevel for each of console and logfile output. +A simple logging module that logs to the console and a logfile. + +Basic usage (print + log both go to terminal and file): -Use it this way: import anuga.utilities.log as log - # configure my logging - log.console_logging_level = log.INFO - log.log_logging_level = log.DEBUG - log.log_filename = './my.log' + log.set_logfile('./my.log') # activates tee to file - # log away! log.debug('A message at DEBUG level') log.info('Another message, INFO level') + print('This also goes to both terminal and file') -This class uses the 'borg' pattern - there is never more than one instance -of log data. See the URL for the basic idea used here: modules *are* -singletons! +Or via the public API: - + import anuga + anuga.set_logfile('./my.log') -Until the first call to log() the user is free to play with the module data -to configure the logging. +Level defaults when a logfile is active: + console: INFO (info/warning/error/critical visible on terminal) + file: DEBUG (everything recorded in the file) +This module uses the 'borg' pattern — modules are singletons. """ import os import sys import traceback import logging +from contextlib import contextmanager -DefaultConsoleLogLevel = logging.CRITICAL -DefaultFileLogLevel = logging.INFO -TimingDelimiter ='#@# ' +DefaultConsoleLogLevel = logging.INFO +DefaultFileLogLevel = logging.DEBUG +TimingDelimiter = '#@# ' ################################################################################ -# Module variables - only one copy of these, ever. -# -# The console logging level is set to a high level, like CRITICAL. The logfile -# logging is set lower, between DEBUG and CRITICAL. The idea is to log least to -# the console, but ensure that everything that goes to the console *will* also -# appear in the log file. There is code to ensure log <= console levels. -# -# If console logging level is set to CRITICAL+1 then nothing will print on the -# console. +# TeeStream — write to both terminal and a log file simultaneously ################################################################################ -# flag variable to determine if logging set up or not -_setup = False +class TeeStream: + """Tee sys.stdout to a file: every write goes to both terminal and file. -# logging level for the console -console_logging_level = DefaultConsoleLogLevel + Usage: + sys.stdout = TeeStream('run.log') + # From now on, print() and sys.stdout.write() go to both places. + sys.stdout.close() # when done (optional) + """ -# logging level for the logfile -log_logging_level = DefaultFileLogLevel + def __init__(self, logfile_path, mode='a'): + self._terminal = sys.__stdout__ + self._log = open(logfile_path, mode, encoding='utf-8') -# The default name of the file to log to. + def write(self, message): + self._terminal.write(message) + self._terminal.flush() + self._log.write(message) + self._log.flush() -def current_datetime(): - from datetime import datetime - return datetime.now().strftime("%Y%m%d_%H%M%S%z") + def flush(self): + self._terminal.flush() + self._log.flush() + + def close(self): + self._log.close() + + # Proxy attribute reads to the underlying terminal so code that inspects + # sys.stdout (e.g. checks for .encoding) still works. + def __getattr__(self, name): + return getattr(self._terminal, name) + + +class _FileOnlyStream: + """Write to the log file only — used by the file_only() context manager.""" + + def __init__(self, log_fh): + self._log = log_fh + def write(self, message): + self._log.write(message) + self._log.flush() -log_filename = os.path.join('.', f'anuga_{current_datetime()}.log') + def flush(self): + self._log.flush() + + def __getattr__(self, name): + return getattr(self._log, name) + + +@contextmanager +def file_only(): + """Context manager: send all print() output to the log file only. + + Terminal output is suppressed for the duration of the block. + Requires set_logfile() to have been called first; if no logfile + is active, output is suppressed entirely. + + Typical use — capture verbose internal output without cluttering + the terminal:: + + with log.file_only(): + anuga.create_pmesh_from_regions(..., verbose=True, ...) + """ + original = sys.stdout + if isinstance(sys.stdout, TeeStream): + sys.stdout = _FileOnlyStream(sys.stdout._log) + else: + # No logfile active — discard output + import io + sys.stdout = io.StringIO() + try: + yield + finally: + sys.stdout = original + + +################################################################################ +# Module variables — only one copy, ever. +################################################################################ + +# flag: has logging been set up yet? +_setup = False + +# logging level for the console handler +console_logging_level = DefaultConsoleLogLevel + +# logging level for the file handler +log_logging_level = DefaultFileLogLevel + +# Path to the log file. None = file logging disabled (no file created). +log_filename = None # set module variables so users don't have to do 'import logging'. CRITICAL = logging.CRITICAL -ERROR = logging.ERROR -WARNING = logging.WARNING -INFO = logging.INFO -DEBUG = logging.DEBUG -NOTSET = logging.NOTSET +ERROR = logging.ERROR +WARNING = logging.WARNING +INFO = logging.INFO +DEBUG = logging.DEBUG +NOTSET = logging.NOTSET + -# set _new_python to True if python version 2.5 or later -_new_python = (sys.version_info[0]*10 + sys.version_info[1] >= 25) # 2.5.x.x +################################################################################ +# set_logfile — the main entry point for enabling file+tee logging +################################################################################ +VERBOSE = logging.DEBUG # level used by log.verbose() — file only by default + + +def set_logfile(path, + console_level=DefaultConsoleLogLevel, + file_level=DefaultFileLogLevel, + verbose_to_screen=False): + """Enable logging to *path*, tee-ing all print() output as well. + + After this call: + - sys.stdout is replaced with a TeeStream so every print() goes to + both the terminal and *path*. + - log.info() writes to both terminal and file. + - log.verbose() / log.debug() write to the file only (unless + verbose_to_screen=True). + - The previous log file (if any) is closed. + + Parameters + ---------- + path : str + File path for the log file. + console_level : int + Logging level for console output (default INFO). + log.verbose() and log.debug() are below this threshold and go + to the file only. + file_level : int + Logging level for file output (default DEBUG — everything). + verbose_to_screen : bool + If True, lower the console threshold to DEBUG so that + log.verbose() output also appears on the terminal. Useful + when debugging without needing a clean screen. + """ + if verbose_to_screen: + console_level = logging.DEBUG + global log_filename, console_logging_level, log_logging_level, _setup + + # Close any existing TeeStream + if isinstance(sys.stdout, TeeStream): + sys.stdout.close() + + log_filename = path + console_logging_level = console_level + log_logging_level = file_level + _setup = False # force re-initialisation on next log() call + + # Tee stdout so print() goes to both terminal and file + sys.stdout = TeeStream(path) + + # Trigger logging setup now + log('Logfile opened: ' + path, INFO) ################################################################################ @@ -90,58 +207,54 @@ def log(msg, level=None): msg: The message string to log. level: The logging level to log with (defaults to console level). - The first call to this method (by anybody) initializes logging and - then logs the message. Subsequent calls just log the message. + The first call to this method initialises the logging.FileHandler if a + log_filename has been configured. ''' global _setup, log_logging_level - fname = '' # default to no frame name if it cannot be found + + fname = '' lnum = 0 - # have we been setup? if not _setup: - # sanity check the logging levels, require console >= file - if log_logging_level > console_logging_level: - log_logging_level = console_logging_level - - # setup the file logging system - if _new_python: + # File logging: only if a filename has been configured + if log_filename is not None: fmt = '%(asctime)s %(levelname)-8s %(mname)25s:%(lnum)-4d|%(message)s' - else: - fmt = '%(asctime)s %(levelname)-8s|%(message)s' - logging.basicConfig(level=log_logging_level, format=fmt, - filename=log_filename, filemode='w') - - # define a console handler which writes to sys.stdout - console = logging.StreamHandler(sys.stdout) - console.setLevel(console_logging_level) - formatter = logging.Formatter('%(message)s') - console.setFormatter(formatter) - logging.getLogger('').addHandler(console) - - # catch exceptions - sys.excepthook = log_exception_hook + file_handler = logging.FileHandler(log_filename, mode='a') + file_handler.setLevel(log_logging_level) + file_handler.setFormatter(logging.Formatter(fmt)) + + root = logging.getLogger('') + root.setLevel(min(log_logging_level, console_logging_level)) + + # Remove any pre-existing handlers to avoid duplicates on re-init + for h in root.handlers[:]: + root.removeHandler(h) - # tell the world how we are set up - start_msg = ("Logfile is '%s' with logging level of %s, " - "console logging level is %s" - % (log_filename, - logging.getLevelName(log_logging_level), - logging.getLevelName(console_logging_level))) - if _new_python: - logging.log(logging.INFO, start_msg, - extra={'mname': __name__, 'lnum': 0}) + root.addHandler(file_handler) + + console = logging.StreamHandler(sys.__stdout__) + console.setLevel(console_logging_level) + console.setFormatter(logging.Formatter('%(message)s')) + root.addHandler(console) else: - logging.log(logging.INFO, start_msg) + # No file configured: just console at console_logging_level + root = logging.getLogger('') + root.setLevel(console_logging_level) + for h in root.handlers[:]: + root.removeHandler(h) + console = logging.StreamHandler(sys.__stdout__) + console.setLevel(console_logging_level) + console.setFormatter(logging.Formatter('%(message)s')) + root.addHandler(console) - # mark module as *setup* + sys.excepthook = log_exception_hook _setup = True - # if logging level not supplied, assume console level if level is None: level = console_logging_level - # get caller information - look back for first module != + # get caller information frames = traceback.extract_stack() frames.reverse() @@ -158,63 +271,48 @@ def log(msg, level=None): if fname != mod_name: break - # why are we here? ... Oh yes! Log the message! - if _new_python: - logging.log(level, msg, extra={'mname': fname, 'lnum': lnum}) - else: - logging.log(level, msg) + logging.log(level, msg, extra={'mname': fname, 'lnum': lnum}) def log_exception_hook(type, value, tb): - '''Hook function to process uncaught exceptions. - - type: Type of exception. - value: The exception data. - tb: Traceback object. - - This has the same interface as sys.excepthook(). - ''' - + '''Hook function to process uncaught exceptions.''' msg = '\n' + ''.join(traceback.format_exception(type, value, tb)) critical(msg) ################################################################################ -# Shortcut routines to make for simpler user code. +# Shortcut routines ################################################################################ -def debug(msg=''): - '''Shortcut for log(DEBUG, msg).''' +def verbose(msg=''): + """Log a verbose/internal message — goes to file only (not screen). + Use this instead of print() inside ANUGA code that has a verbose flag. + Output appears on screen only when set_logfile(..., verbose_to_screen=True). + """ log(msg, logging.DEBUG) +def debug(msg=''): + """Log a DEBUG-level message (file only by default).""" + log(msg, logging.DEBUG) def info(msg=''): - '''Shortcut for log(INFO, msg).''' - + """Log an INFO-level message (terminal and file).""" log(msg, logging.INFO) - def warning(msg=''): - '''Shortcut for log(WARNING, msg).''' - + """Log a WARNING-level message (terminal and file).""" log(msg, logging.WARNING) - def error(msg=''): - '''Shortcut for log(ERROR, msg).''' - + """Log an ERROR-level message (terminal and file).""" log(msg, logging.ERROR) - def critical(msg=''): - '''Shortcut for log(CRITICAL, msg).''' - + """Log a CRITICAL-level message (terminal and file).""" log(msg, logging.CRITICAL) def timingInfo(msg=''): - '''Shortcut for log(timingDelimiter, msg).''' - log(TimingDelimiter + msg, logging.INFO) @@ -228,57 +326,44 @@ def resource_usage(level=logging.INFO): _proc_status = '/proc/%d/status' % os.getpid() def _VmB(VmKey): - '''Get number of virtual bytes used.''' - - # get pseudo file /proc//status try: t = open(_proc_status) v = t.read() t.close() except IOError: return 0.0 - - # get VmKey line, eg: 'VmRSS: 999 kB\n ... i = v.index(VmKey) v = v[i:].split(None, 3) if len(v) < 3: return 0.0 - - # convert Vm value to bytes return float(v[1]) * _scale[v[2]] def memory(since=0.0): - '''Get virtual memory usage in bytes.''' - return _VmB('VmSize:') - since def resident(since=0.0): - '''Get resident memory usage in bytes.''' - return _VmB('VmRSS:') - since def stacksize(since=0.0): - '''Get stack size in bytes.''' - return _VmB('VmStk:') - since msg = ('Resource usage: memory=%.1fMB resident=%.1fMB stacksize=%.1fMB' - % ((memory() / _scale['MB']), - (resident() / _scale['MB']), - (stacksize() / _scale['MB']))) + % (memory() / _scale['MB'], + resident() / _scale['MB'], + stacksize() / _scale['MB'])) log(msg, level) else: - # Windows code from: http://code.activestate.com/recipes/511491/ try: import ctypes import winreg except ImportError: - log(level, 'Windows resource usage not available') + log('Windows resource usage not available', level) return kernel32 = ctypes.windll.kernel32 c_ulong = ctypes.c_ulong c_ulonglong = ctypes.c_ulonglong + class MEMORYSTATUSEX(ctypes.Structure): _fields_ = [('dwLength', c_ulong), ('dwMemoryLoad', c_ulong), @@ -288,18 +373,22 @@ class MEMORYSTATUSEX(ctypes.Structure): ('ullAvailPageFile', c_ulonglong), ('ullTotalVirtual', c_ulonglong), ('ullAvailVirtual', c_ulonglong), - ('ullAvailExtendedVirtual', c_ulonglong) - ] + ('ullAvailExtendedVirtual', c_ulonglong)] memoryStatusEx = MEMORYSTATUSEX() memoryStatusEx.dwLength = ctypes.sizeof(MEMORYSTATUSEX) kernel32.GlobalMemoryStatusEx(ctypes.byref(memoryStatusEx)) msg = ('Resource usage: total memory=%.1fMB free memory=%.1fMB' - % ((memoryStatusEx.ullTotalPhys / _scale['MB']), - (memoryStatusEx.ullAvailPhys / _scale['MB']))) + % (memoryStatusEx.ullTotalPhys / _scale['MB'], + memoryStatusEx.ullAvailPhys / _scale['MB'])) log(msg, level) + +def current_datetime(): + from datetime import datetime + return datetime.now().strftime("%Y%m%d_%H%M%S%z") + def CurrentDateTime(): from datetime import datetime return datetime.now().strftime("%Y-%m-%d_%H:%M:%S") @@ -309,8 +398,8 @@ def TimeStamp(): return datetime.now().strftime('%Y%m%d_%H%M%S') -def resource_usage_timing(level=logging.INFO, prefix =""): - '''Log memory usage at given log level.''' +def resource_usage_timing(level=logging.INFO, prefix=''): + '''Log memory usage with timing info.''' _scale = {'KB': 1024, 'MB': 1024*1024, 'GB': 1024*1024*1024, 'kB': 1024, 'mB': 1024*1024, 'gB': 1024*1024*1024} @@ -319,61 +408,43 @@ def resource_usage_timing(level=logging.INFO, prefix =""): _proc_status = '/proc/%d/status' % os.getpid() def _VmB(VmKey): - '''Get number of virtual bytes used.''' - - # get pseudo file /proc//status try: t = open(_proc_status) v = t.read() t.close() except IOError: return 0.0 - - # get VmKey line, eg: 'VmRSS: 999 kB\n ... i = v.index(VmKey) v = v[i:].split(None, 3) if len(v) < 3: return 0.0 - - # convert Vm value to bytes return float(v[1]) * _scale[v[2]] - def memory(since=0.0): - '''Get virtual memory usage in bytes.''' - - return _VmB('VmSize:') - since - - def resident(since=0.0): - '''Get resident memory usage in bytes.''' - - return _VmB('VmRSS:') - since - - def stacksize(since=0.0): - '''Get stack size in bytes.''' - - return _VmB('VmStk:') - since + memory = lambda since=0.0: _VmB('VmSize:') - since + resident = lambda since=0.0: _VmB('VmRSS:') - since + stacksize= lambda since=0.0: _VmB('VmStk:') - since msg = ('Resource usage: memory=%.1fMB resident=%.1fMB stacksize=%.1fMB' - % ((memory() / _scale['MB']), - (resident() / _scale['MB']), - (stacksize() / _scale['MB']))) + % (memory() / _scale['MB'], + resident() / _scale['MB'], + stacksize() / _scale['MB'])) log(msg, level) timingInfo('sys_platform, ' + sys.platform) - timingInfo(prefix + 'memory, ' + str((memory() / _scale['MB']))) - timingInfo(prefix + 'resident, ' + str((resident() / _scale['MB']))) - timingInfo(prefix + 'stacksize, ' + str((stacksize() / _scale['MB']))) + timingInfo(prefix + 'memory, ' + str(memory() / _scale['MB'])) + timingInfo(prefix + 'resident, ' + str(resident() / _scale['MB'])) + timingInfo(prefix + 'stacksize, ' + str(stacksize() / _scale['MB'])) else: - # Windows code from: http://code.activestate.com/recipes/511491/ try: import ctypes import winreg except ImportError: - log(level, 'Windows resource usage not available') + log('Windows resource usage not available', level) return kernel32 = ctypes.windll.kernel32 c_ulong = ctypes.c_ulong c_ulonglong = ctypes.c_ulonglong + class MEMORYSTATUSEX(ctypes.Structure): _fields_ = [('dwLength', c_ulong), ('dwMemoryLoad', c_ulong), @@ -383,42 +454,26 @@ class MEMORYSTATUSEX(ctypes.Structure): ('ullAvailPageFile', c_ulonglong), ('ullTotalVirtual', c_ulonglong), ('ullAvailVirtual', c_ulonglong), - ('ullAvailExtendedVirtual', c_ulonglong) - ] + ('ullAvailExtendedVirtual', c_ulonglong)] memoryStatusEx = MEMORYSTATUSEX() memoryStatusEx.dwLength = ctypes.sizeof(MEMORYSTATUSEX) kernel32.GlobalMemoryStatusEx(ctypes.byref(memoryStatusEx)) msg = ('Resource usage: total memory=%.1fMB free memory=%.1fMB' - % ((memoryStatusEx.ullTotalPhys / _scale['MB']), - (memoryStatusEx.ullAvailPhys / _scale['MB']))) + % (memoryStatusEx.ullTotalPhys / _scale['MB'], + memoryStatusEx.ullAvailPhys / _scale['MB'])) log(msg, level) timingInfo('sys_platform, ' + sys.platform) - timingInfo(prefix + 'total_memory, ' + str((memoryStatusEx.ullTotalPhys / _scale['MB']))) - timingInfo(prefix + 'free_memory, ' + str((memoryStatusEx.ullAvailPhys / _scale['MB']))) + timingInfo(prefix + 'total_memory, ' + str(memoryStatusEx.ullTotalPhys / _scale['MB'])) + timingInfo(prefix + 'free_memory, ' + str(memoryStatusEx.ullAvailPhys / _scale['MB'])) ################################################################################ if __name__ == '__main__': + set_logfile('/tmp/anuga_test.log') critical('#' * 80) warning('Test of logging...') - log('CRITICAL+1', CRITICAL+1) - log('CRITICAL', CRITICAL) - log('CRITICAL-1', CRITICAL-1) - log('CRITICAL-2', CRITICAL-2) - log('default - CRITICAL?') - - def test_it(num=100): - if num > 0: - test_it(num-1) - else: - resource_usage() - - import numpy as num - - a = num.zeros((1000,1000), float) - - info('sys.version_info=%s, _new_python=%s' - % (str(sys.version_info), str(_new_python))) - test_it() + info('An info message') + debug('A debug message (file only if console level is INFO)') + print('This print() goes to both terminal and /tmp/anuga_test.log') diff --git a/archive/cupy_cuda/README.md b/archive/cupy_cuda/README.md new file mode 100644 index 000000000..c5b320506 --- /dev/null +++ b/archive/cupy_cuda/README.md @@ -0,0 +1,22 @@ +# CuPy/CUDA archive + +These files implement an earlier GPU acceleration approach using CuPy as a +NumPy drop-in replacement, alongside a raw CUDA C kernel (`cuda_anuga.cu`). + +They were moved out of `anuga/shallow_water/` to reduce clutter. The active +GPU/OpenMP offloading work (for SC26) lives in the main `develop` branch under +`anuga/shallow_water/sw_domain_openmp_ext.pyx` and the associated C headers. + +## Contents + +| File | Description | +|------|-------------| +| `sw_domain_cupy.py` | Domain subclass using CuPy arrays on the GPU | +| `sw_domain_cuda.py` | Domain subclass driving the raw CUDA kernel | +| `cuda_anuga.cu` | Raw CUDA C kernel for flux/update computations | +| `tests/` | CuPy-specific test scripts and pytest file | + +## Status + +Archived 2026-04-05. Not wired into meson, not imported by `anuga/__init__.py`. +Preserved in case the CuPy approach is revisited. diff --git a/anuga/shallow_water/cuda_anuga.cu b/archive/cupy_cuda/cuda_anuga.cu similarity index 100% rename from anuga/shallow_water/cuda_anuga.cu rename to archive/cupy_cuda/cuda_anuga.cu diff --git a/anuga/shallow_water/sw_domain_cuda.py b/archive/cupy_cuda/sw_domain_cuda.py similarity index 100% rename from anuga/shallow_water/sw_domain_cuda.py rename to archive/cupy_cuda/sw_domain_cuda.py diff --git a/anuga/shallow_water/sw_domain_cupy.py b/archive/cupy_cuda/sw_domain_cupy.py similarity index 100% rename from anuga/shallow_water/sw_domain_cupy.py rename to archive/cupy_cuda/sw_domain_cupy.py diff --git a/anuga/shallow_water/cupy_tests/run_cupy_cft.py b/archive/cupy_cuda/tests/run_cupy_cft.py similarity index 100% rename from anuga/shallow_water/cupy_tests/run_cupy_cft.py rename to archive/cupy_cuda/tests/run_cupy_cft.py diff --git a/anuga/shallow_water/cupy_tests/run_cupy_extrapolation.py b/archive/cupy_cuda/tests/run_cupy_extrapolation.py similarity index 100% rename from anuga/shallow_water/cupy_tests/run_cupy_extrapolation.py rename to archive/cupy_cuda/tests/run_cupy_extrapolation.py diff --git a/anuga/shallow_water/cupy_tests/run_cupy_flux_distribute.py b/archive/cupy_cuda/tests/run_cupy_flux_distribute.py similarity index 100% rename from anuga/shallow_water/cupy_tests/run_cupy_flux_distribute.py rename to archive/cupy_cuda/tests/run_cupy_flux_distribute.py diff --git a/anuga/shallow_water/cupy_tests/run_cupy_protect_negative.py b/archive/cupy_cuda/tests/run_cupy_protect_negative.py similarity index 100% rename from anuga/shallow_water/cupy_tests/run_cupy_protect_negative.py rename to archive/cupy_cuda/tests/run_cupy_protect_negative.py diff --git a/anuga/shallow_water/cupy_tests/run_cupy_update_conserved_quantities.py b/archive/cupy_cuda/tests/run_cupy_update_conserved_quantities.py similarity index 100% rename from anuga/shallow_water/cupy_tests/run_cupy_update_conserved_quantities.py rename to archive/cupy_cuda/tests/run_cupy_update_conserved_quantities.py diff --git a/anuga/shallow_water/cupy_tests/run_sw_domain_cupy.py b/archive/cupy_cuda/tests/run_sw_domain_cupy.py similarity index 100% rename from anuga/shallow_water/cupy_tests/run_sw_domain_cupy.py rename to archive/cupy_cuda/tests/run_sw_domain_cupy.py diff --git a/anuga/shallow_water/cupy_tests/test_sw_domain_cupy.py b/archive/cupy_cuda/tests/test_sw_domain_cupy.py similarity index 100% rename from anuga/shallow_water/cupy_tests/test_sw_domain_cupy.py rename to archive/cupy_cuda/tests/test_sw_domain_cupy.py diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 000000000..8529a8d7d --- /dev/null +++ b/benchmarks/README.md @@ -0,0 +1,112 @@ +# ANUGA Benchmarks + +Tools to track wall-time and memory performance across commits. + +## Quick start + +```bash +# Activate your env +conda activate anuga_env_3.14 + +# Run small + medium scenarios, all modes (takes ~40 s) +python benchmarks/run_benchmarks.py + +# Quick sanity check — small only (~5 s) +python benchmarks/run_benchmarks.py --sizes small + +# List saved result files +python benchmarks/compare_benchmarks.py --list + +# Compare two results +python benchmarks/compare_benchmarks.py benchmarks/results/before.json \ + benchmarks/results/after.json +``` + +## Scenarios + +| Name | Triangles | finaltime | Typical wall (mode 0) | +|--------|----------:|----------:|----------------------:| +| small | 10 000 | 200 s | ~1 s | +| medium | 90 000 | 100 s | ~16 s | +| large | 360 000 | 50 s | ~90 s | + +`large` is not run by default. Use `--sizes small,medium,large` to include it. + +## Modes + +| Mode | Description | +|------|-------------| +| 0 | Python Euler (default) | +| 1 | Python RK2 | +| 2 | C RK2 / GPU (CPU_ONLY_MODE if no GPU present) | + +## Metrics + +| Metric | Meaning | +|-------------|---------| +| `cells/s` | n_triangles × n_steps / wall_time — primary performance figure | +| `setup MB` | RSS after domain creation (before any evolve) | +| `peak MB` | Peak RSS sampled at 100 ms intervals during evolve | +| `MB/Ktri` | peak_MB / (n_triangles / 1000) — memory per 1 000 triangles | + +## Workflow: before/after comparison + +```bash +# 1. Capture baseline on current commit +python benchmarks/run_benchmarks.py --output /tmp/before.json + +# 2. Make your code changes ... + +# 3. Capture new result +python benchmarks/run_benchmarks.py --output /tmp/after.json + +# 4. Compare +python benchmarks/compare_benchmarks.py /tmp/before.json /tmp/after.json + +# Show only rows that changed by more than 5% +python benchmarks/compare_benchmarks.py /tmp/before.json /tmp/after.json --threshold 5 +``` + +## Result files + +Results are saved to `benchmarks/results/__.json`. +The `results/` directory is git-ignored so committed baselines don't pollute +the repo. Copy important baselines elsewhere if you want to keep them. + +--- + +## Parallel distribution benchmarks (MPI) + +`distribute_benchmarks.py` compares four approaches for distributing a mesh +across MPI ranks. Requires `mpi4py`. + +```bash +# Run all four methods on a 500×500 mesh (~1M tris) with 8 ranks +mpirun -np 8 python benchmarks/distribute_benchmarks.py --size 500 + +# Morton scheme, 3 repetitions for stable medians +mpirun -np 8 python benchmarks/distribute_benchmarks.py --size 500 --scheme morton --reps 3 +``` + +| Method | Description | +|--------|-------------| +| `distribute()` | Traditional: full Domain on rank 0, then distribute | +| `distribute_collaborative()` | Shared-memory cooperative version | +| `distribute_basic_mesh()` | Mesh-first: only Basic_mesh on rank 0, quantities set locally after | +| `dump()+load()` | Rank 0 partitions and writes files; all ranks read | + +### Grid sweep (multiple np × scheme combinations) + +```bash +# Run the full np=[10,20,30] × scheme=[metis,morton,hilbert] grid +python benchmarks/run_benchmark_grid.py --size 1000 + +# Custom grid +python benchmarks/run_benchmark_grid.py --size 500 --np 4,8,16 --schemes metis,morton + +# Dry run (print commands without executing) +python benchmarks/run_benchmark_grid.py --dry-run +``` + +Results are saved to `benchmarks/results/dist/bench_np_.txt` and +a consolidated `summary.txt`. diff --git a/benchmarks/compare_benchmarks.py b/benchmarks/compare_benchmarks.py new file mode 100644 index 000000000..7a94f6a4f --- /dev/null +++ b/benchmarks/compare_benchmarks.py @@ -0,0 +1,212 @@ +#!/usr/bin/env python3 +""" +ANUGA Benchmark Comparison Tool +--------------------------------- +Compares two JSON results files produced by run_benchmarks.py and prints +a side-by-side delta table. + +Usage +----- + python benchmarks/compare_benchmarks.py before.json after.json + + # Show only speed changes > 5% + python benchmarks/compare_benchmarks.py before.json after.json --threshold 5 + + # List available result files + python benchmarks/compare_benchmarks.py --list +""" + +import argparse +import json +import os +import sys + + +def _load(path): + with open(path) as fh: + return json.load(fh) + + +def _pct(old, new): + """Return signed % change from old to new, or None if old==0.""" + if old == 0: + return None + return (new - old) / old * 100.0 + + +def _fmt_pct(pct): + """Format a signed % change.""" + if pct is None: + return ' n/a' + return f'{pct:+.1f}%' + + +def _arrow(pct, invert=False): + """▲/▼/= based on direction and whether lower-is-better.""" + if pct is None or abs(pct) < 0.5: + return ' ' + improved = (pct < 0) if invert else (pct > 0) + return '↑' if improved else '↓' + + +def compare(before, after, threshold=0.0): + b_results = {r['name']: r for r in before['results']} + a_results = {r['name']: r for r in after['results']} + + all_names = sorted(set(b_results) | set(a_results)) + + print(f"\nBenchmark comparison") + print(f" Before : {before['git']['commit']} branch={before['git']['branch']} " + f"t={before['timestamp']}") + print(f" After : {after['git']['commit']} branch={after['git']['branch']} " + f"t={after['timestamp']}") + print() + + # Header + col = '{:<28} {:>7} {:>4} {:>10} {:>8} {:>10} {:>8} {:>10} {:>8}' + hdr = col.format( + 'Scenario', 'tris', 'mode', + 'wall(s)', 'Δwall', + 'cells/s', 'Δspeed', + 'peak MB', 'Δmem', + ) + rule = '─' * len(hdr) + print(hdr) + print(rule) + + any_printed = False + for name in all_names: + b = b_results.get(name) + a = a_results.get(name) + + if b is None: + print(f' {name:<28} (only in after)') + any_printed = True + continue + if a is None: + print(f' {name:<28} (only in before)') + any_printed = True + continue + + d_wall = _pct(b['wall_time_s'], a['wall_time_s']) + d_speed = _pct(b['cells_per_s'], a['cells_per_s']) + d_mem = _pct(b['peak_rss_mb'], a['peak_rss_mb']) + + # Skip if all changes are below threshold + def _mag(v): + return abs(v) if v is not None else 0.0 + + if threshold > 0 and max(_mag(d_wall), _mag(d_speed), _mag(d_mem)) < threshold: + continue + + w_arrow = _arrow(d_wall, invert=True) + s_arrow = _arrow(d_speed, invert=False) + m_arrow = _arrow(d_mem, invert=True) + + print(col.format( + name, + f"{b['n_triangles']:,}", b['mode'], + f"{a['wall_time_s']:.2f}", + f"{_fmt_pct(d_wall)}{w_arrow}", + f"{a['cells_per_s']:,.0f}", + f"{_fmt_pct(d_speed)}{s_arrow}", + f"{a['peak_rss_mb']:.1f}", + f"{_fmt_pct(d_mem)}{m_arrow}", + )) + any_printed = True + + if not any_printed: + print(f' (no changes exceed {threshold}% threshold)') + + print(rule) + print() + + # Summary statistics + common = [(b_results[n], a_results[n]) + for n in all_names if n in b_results and n in a_results] + if common: + speed_pcts = [_pct(b['cells_per_s'], a['cells_per_s']) + for b, a in common if b['cells_per_s'] > 0] + mem_pcts = [_pct(b['peak_rss_mb'], a['peak_rss_mb']) + for b, a in common if b['peak_rss_mb'] > 0] + + def _mean(lst): + lst = [v for v in lst if v is not None] + return sum(lst) / len(lst) if lst else None + + ms = _mean(speed_pcts) + mm = _mean(mem_pcts) + + if ms is not None: + dir_s = 'faster' if ms > 0 else 'slower' + print(f' Average speed change : {ms:+.1f}% ({dir_s})') + if mm is not None: + dir_m = 'less memory' if mm < 0 else 'more memory' + print(f' Average memory change: {mm:+.1f}% ({dir_m})') + print() + + +def list_results(): + results_dir = os.path.join(os.path.dirname(__file__), 'results') + if not os.path.isdir(results_dir): + print('No results/ directory found. Run run_benchmarks.py first.') + return + + files = sorted( + f for f in os.listdir(results_dir) if f.endswith('.json') + ) + if not files: + print('No result files in benchmarks/results/') + return + + print(f'Result files in {results_dir}:') + print() + col = ' {:<50} {:>10} {:>10} {:>6}' + print(col.format('File', 'Commit', 'Branch', 'Runs')) + print(' ' + '─' * 80) + for fname in files: + path = os.path.join(results_dir, fname) + try: + data = _load(path) + commit = data.get('git', {}).get('commit', '?') + branch = data.get('git', {}).get('branch', '?') + n = len(data.get('results', [])) + print(col.format(fname, commit, branch, n)) + except Exception as exc: + print(f' {fname:<50} (error: {exc})') + print() + + +def main(): + parser = argparse.ArgumentParser( + description='Compare two ANUGA benchmark result files.', + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=__doc__, + ) + parser.add_argument('before', nargs='?', help='Before JSON file') + parser.add_argument('after', nargs='?', help='After JSON file') + parser.add_argument( + '--threshold', type=float, default=0.0, + help='Only show rows where any metric changes by at least this %% (default: 0)', + ) + parser.add_argument( + '--list', action='store_true', + help='List available result files in benchmarks/results/', + ) + args = parser.parse_args() + + if args.list: + list_results() + return 0 + + if not args.before or not args.after: + parser.error('Provide both before and after JSON files, or use --list') + + before = _load(args.before) + after = _load(args.after) + compare(before, after, threshold=args.threshold) + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/benchmarks/distribute_benchmarks.py b/benchmarks/distribute_benchmarks.py new file mode 100644 index 000000000..6b46cd433 --- /dev/null +++ b/benchmarks/distribute_benchmarks.py @@ -0,0 +1,574 @@ +#!/usr/bin/env python3 +"""Benchmark ANUGA parallel mesh distribution methods. + +Compares four approaches for distributing a mesh across MPI ranks: + + distribute() -- traditional: full Domain on rank 0, then distribute + distribute_collaborative() -- shared-memory cooperative version of distribute() + distribute_basic_mesh() -- mesh-first: only Basic_mesh on rank 0 (no quantities), + quantities set locally after distribution + dump()+load() -- rank 0 partitions and writes files; all ranks read + +Run with: + + mpirun -np 8 python benchmarks/distribute_benchmarks.py --size 500 + mpirun -np 8 python benchmarks/distribute_benchmarks.py --size 500 --scheme morton + mpirun -np 8 python benchmarks/distribute_benchmarks.py --size 500 --reps 3 + +Options +------- +--size M Grid size M: produces 4*M*M triangles (default 500) +--reps R Repetitions for median timing (default 1) +--scheme S Partition scheme: metis | morton | hilbert (default metis) +--interval T Memory ticker sample interval in seconds (default 1.0) +--no-evolve Skip the correctness evolve check after timing +""" + +import argparse +import gc +import shutil +import statistics +import tempfile +import threading +import time + +from mpi4py import MPI + +import anuga +from anuga.parallel.parallel_api import distribute, distribute_collaborative +from anuga.parallel.parallel_api import distribute_basic_mesh +from anuga.parallel.sequential_distribute import ( + sequential_distribute_dump, sequential_distribute_load) + +comm = MPI.COMM_WORLD +myid = comm.Get_rank() +nproc = comm.Get_size() + + +# --------------------------------------------------------------------------- +# Memory helpers +# --------------------------------------------------------------------------- + +def _rss_mb(): + """Current RSS in MiB from /proc/self/status.""" + try: + with open('/proc/self/status') as fh: + for line in fh: + if line.startswith('VmRSS:'): + return int(line.split()[1]) / 1024.0 + except OSError: + pass + return 0.0 + + +def _pss_mb(): + """Proportional Set Size in MiB (shared pages counted proportionally). + + Summing PSS across all ranks gives the true physical-memory footprint + of the job, unlike VmRSS which double-counts pages shared via MPI.Win. + """ + try: + with open('/proc/self/smaps_rollup') as fh: + for line in fh: + if line.startswith('Pss:'): + return int(line.split()[1]) / 1024.0 + except OSError: + pass + try: + total = 0 + with open('/proc/self/smaps') as fh: + for line in fh: + if line.startswith('Pss:'): + total += int(line.split()[1]) + return total / 1024.0 + except OSError: + pass + return _rss_mb() + + +class MemoryMonitor: + """Background thread sampling this process's RSS/PSS every `interval` s. + + Use as a context manager:: + + with MemoryMonitor(label='distribute()', interval=1.0) as mon: + fn(domain) + peak_pss_mib = mon.peak_pss_mb + """ + + def __init__(self, label='', interval=1.0, print_ticker=False): + self.label = label + self.interval = interval + self.print_ticker = print_ticker + self.peak_rss_mb = 0.0 + self.peak_pss_mb = 0.0 + self._stop = threading.Event() + self._thread = threading.Thread(target=self._run, daemon=True) + + def __enter__(self): + self.peak_rss_mb = _rss_mb() + self.peak_pss_mb = _pss_mb() + self._t0 = time.time() + self._stop.clear() + self._thread.start() + return self + + def __exit__(self, *_): + self._stop.set() + self._thread.join() + + def _run(self): + while not self._stop.wait(self.interval): + rss = _rss_mb() + pss = _pss_mb() + if rss > self.peak_rss_mb: + self.peak_rss_mb = rss + if pss > self.peak_pss_mb: + self.peak_pss_mb = pss + if self.print_ticker: + elapsed = time.time() - self._t0 + print(f' [{self.label}] t={elapsed:6.1f}s ' + f'rank-0 RSS={rss:,.0f} MiB PSS={pss:,.0f} MiB', + flush=True) + + +# --------------------------------------------------------------------------- +# Shared-memory diagnostic +# --------------------------------------------------------------------------- + +def check_shmem(): + """Test whether MPI.Win.Allocate_shared works. Returns (ok, node_size, reason).""" + import numpy as np + node_comm = comm.Split_type(MPI.COMM_TYPE_SHARED) + node_rank = node_comm.Get_rank() + node_size = node_comm.Get_size() + try: + win = MPI.Win.Allocate_shared( + 8 if node_rank == 0 else 0, 8, MPI.INFO_NULL, node_comm) + buf, _ = win.Shared_query(0) + arr = np.ndarray((1,), dtype=np.float64, buffer=buf) + if node_rank == 0: + arr[0] = 42.0 + node_comm.Barrier() + ok = float(arr[0]) == 42.0 + win.Free() + node_comm.Free() + return ok, node_size, ('Win.Allocate_shared OK' if ok + else 'shared query returned wrong value') + except Exception as e: + try: + node_comm.Free() + except Exception: + pass + return False, node_size, f'exception: {e}' + + +# --------------------------------------------------------------------------- +# Domain / mesh factories +# --------------------------------------------------------------------------- + +def make_domain(grid_size): + """Build a full Domain on rank 0 with a simple initial condition.""" + pts, verts, bnd = anuga.rectangular_cross(grid_size, grid_size, + len1=float(grid_size), + len2=float(grid_size)) + domain = anuga.Domain(pts, verts, bnd) + domain.set_quantity('elevation', 0.0) + domain.set_quantity('stage', lambda x, y: 1.0 + 0.01 * x) + domain.set_quantity('friction', 0.03) + return domain + + +def make_basic_mesh(grid_size): + """Build a Basic_mesh on rank 0 only (no quantities).""" + if myid != 0: + return None + from anuga.abstract_2d_finite_volumes.basic_mesh import ( + rectangular_cross_basic_mesh) + return rectangular_cross_basic_mesh(grid_size, grid_size, + len1=float(grid_size), + len2=float(grid_size)) + + +# --------------------------------------------------------------------------- +# Ghost-triangle statistics +# --------------------------------------------------------------------------- + +def collect_ghost_stats(pd): + """Return (ghost_sum, ghost_max, ghost_min) reduced to rank 0.""" + import numpy as np + n_ghost = int(np.sum(pd.tri_full_flag == 0)) + g_sum = comm.reduce(n_ghost, op=MPI.SUM, root=0) + g_max = comm.reduce(n_ghost, op=MPI.MAX, root=0) + g_min = comm.reduce(n_ghost, op=MPI.MIN, root=0) + return (g_sum or 0), (g_max or 0), (g_min or 0) + + +# --------------------------------------------------------------------------- +# Timed benchmarks +# --------------------------------------------------------------------------- + +def time_distribute(fn, grid_size, reps, ticker_interval, parameters=None): + """Time distribute() or distribute_collaborative(). + + Returns (times_s, pss_sum_mbs, rss_max_mbs, ghost_stats) where values + are only valid on rank 0. + """ + times = [] + pss_sum_mbs = [] + rss_max_mbs = [] + ghost_stats = (0, 0, 0) + label = fn.__name__ + + for _ in range(reps): + gc.collect() + comm.Barrier() + domain = make_domain(grid_size) + comm.Barrier() + + mon = MemoryMonitor(label=label, interval=ticker_interval, + print_ticker=(myid == 0)) + with mon: + t0 = MPI.Wtime() + pd = fn(domain, parameters=parameters) + t1 = MPI.Wtime() + + elapsed = t1 - t0 + peak_rss = max(mon.peak_rss_mb, _rss_mb()) + peak_pss = max(mon.peak_pss_mb, _pss_mb()) + ghost_stats = collect_ghost_stats(pd) + + wall = comm.reduce(elapsed, op=MPI.MAX, root=0) + pss_sum = comm.reduce(peak_pss, op=MPI.SUM, root=0) + rss_max = comm.reduce(peak_rss, op=MPI.MAX, root=0) + + if myid == 0: + times.append(wall) + pss_sum_mbs.append(pss_sum) + rss_max_mbs.append(rss_max) + + del pd, domain + + return times, pss_sum_mbs, rss_max_mbs, ghost_stats + + +def time_distribute_basic_mesh(grid_size, reps, ticker_interval, parameters=None): + """Time distribute_basic_mesh() -- mesh-first workflow. + + Rank 0 builds only a Basic_mesh (no quantities). All ranks call + distribute_basic_mesh(), then set quantities locally on the result. + + Returns (times_s, pss_sum_mbs, rss_max_mbs, ghost_stats). + """ + times = [] + pss_sum_mbs = [] + rss_max_mbs = [] + ghost_stats = (0, 0, 0) + + for _ in range(reps): + gc.collect() + comm.Barrier() + bm = make_basic_mesh(grid_size) + comm.Barrier() + + mon = MemoryMonitor(label='distribute_basic_mesh', interval=ticker_interval, + print_ticker=(myid == 0)) + with mon: + t0 = MPI.Wtime() + pd = distribute_basic_mesh(bm, parameters=parameters) + t1 = MPI.Wtime() + + # Set quantities post-distribution (part of the workflow) + pd.set_quantity('elevation', 0.0) + pd.set_quantity('stage', lambda x, y: 1.0 + 0.01 * x) + pd.set_quantity('friction', 0.03) + + elapsed = t1 - t0 + peak_rss = max(mon.peak_rss_mb, _rss_mb()) + peak_pss = max(mon.peak_pss_mb, _pss_mb()) + ghost_stats = collect_ghost_stats(pd) + + wall = comm.reduce(elapsed, op=MPI.MAX, root=0) + pss_sum = comm.reduce(peak_pss, op=MPI.SUM, root=0) + rss_max = comm.reduce(peak_rss, op=MPI.MAX, root=0) + + if myid == 0: + times.append(wall) + pss_sum_mbs.append(pss_sum) + rss_max_mbs.append(rss_max) + + del pd, bm + gc.collect() + comm.Barrier() + + return times, pss_sum_mbs, rss_max_mbs, ghost_stats + + +def time_dump_load(grid_size, reps, ticker_interval, parameters=None): + """Time sequential_distribute_dump (rank 0) + load (all ranks). + + Returns (times_dump, times_load, pss_sum_mbs, rss_max_mbs, ghost_stats). + """ + times_dump = [] + times_load = [] + pss_sum_mbs = [] + rss_max_mbs = [] + ghost_stats = (0, 0, 0) + + for _ in range(reps): + gc.collect() + comm.Barrier() + + if myid == 0: + tmpdir = tempfile.mkdtemp(prefix='anuga_bench_') + else: + tmpdir = None + tmpdir = comm.bcast(tmpdir, root=0) + + domain = make_domain(grid_size) + domain_name = domain.get_name() + comm.Barrier() + + # Dump phase (rank 0 only) + if myid == 0: + mon_dump = MemoryMonitor(label='dump', interval=ticker_interval, + print_ticker=True) + with mon_dump: + t0 = MPI.Wtime() + sequential_distribute_dump(domain, nproc, + partition_dir=tmpdir, + parameters=parameters) + t1 = MPI.Wtime() + dump_elapsed = t1 - t0 + else: + dump_elapsed = 0.0 + + del domain + comm.Barrier() + + # Load phase (all ranks) + mon_load = MemoryMonitor(label='load', interval=ticker_interval, + print_ticker=(myid == 0)) + with mon_load: + t0 = MPI.Wtime() + pd = sequential_distribute_load(filename=domain_name, + partition_dir=tmpdir) + t1 = MPI.Wtime() + + load_elapsed = t1 - t0 + peak_rss = max(mon_load.peak_rss_mb, _rss_mb()) + peak_pss = max(mon_load.peak_pss_mb, _pss_mb()) + ghost_stats = collect_ghost_stats(pd) + + dump_wall = comm.bcast(dump_elapsed, root=0) + load_wall = comm.reduce(load_elapsed, op=MPI.MAX, root=0) + pss_sum = comm.reduce(peak_pss, op=MPI.SUM, root=0) + rss_max = comm.reduce(peak_rss, op=MPI.MAX, root=0) + + if myid == 0: + times_dump.append(dump_wall) + times_load.append(load_wall) + pss_sum_mbs.append(pss_sum) + rss_max_mbs.append(rss_max) + + del pd + + if myid == 0: + shutil.rmtree(tmpdir, ignore_errors=True) + comm.Barrier() + + return times_dump, times_load, pss_sum_mbs, rss_max_mbs, ghost_stats + + +def run_evolve_check(grid_size, scheme): + """Distribute a small mesh and run one evolve step to verify correctness.""" + import numpy as np + + small = min(grid_size, 50) + if myid == 0: + print(f'\n Evolve check ({small}x{small} mesh, distribute_basic_mesh) ...') + + bm = make_basic_mesh(small) + pd = distribute_basic_mesh(bm, parameters={'partition_scheme': scheme}) + + pd.set_quantity('elevation', 0.0) + pd.set_quantity('stage', lambda x, y: np.where( + (x > small * 0.4) & (x < small * 0.6), 1.0, 0.0)) + pd.set_quantity('friction', 0.03) + Br = anuga.Reflective_boundary(pd) + pd.set_boundary({t: Br for t in pd.get_boundary_tags()}) + + t0 = time.perf_counter() + for _ in pd.evolve(yieldstep=0.01, finaltime=0.01): + pass + dt = time.perf_counter() - t0 + + if myid == 0: + print(f' Evolve check passed in {dt:.2f}s') + + +# --------------------------------------------------------------------------- +# Formatting helpers +# --------------------------------------------------------------------------- + +def fmt_time(times): + if not times: + return '--' + med = statistics.median(times) + if len(times) == 1: + return f'{med:.3f}s' + return f'{med:.3f}s (min {min(times):.3f} max {max(times):.3f})' + + +def fmt_mem(mbs): + if not mbs: + return '--' + med = statistics.median(mbs) + return f'{med/1024:.2f} GiB' if med >= 1024 else f'{med:.0f} MiB' + + +# --------------------------------------------------------------------------- +# Argument parsing +# --------------------------------------------------------------------------- + +def parse_args(): + p = argparse.ArgumentParser( + description='Benchmark ANUGA parallel mesh distribution methods.', + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=__doc__, + ) + p.add_argument('--size', type=int, default=500) + p.add_argument('--reps', type=int, default=1) + p.add_argument('--scheme', type=str, default='metis', + choices=['metis', 'morton', 'hilbert']) + p.add_argument('--interval', type=float, default=1.0) + p.add_argument('--no-evolve', action='store_true') + return p.parse_args() + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + args = parse_args() + M = args.size + ntri = 4 * M * M + params = {'partition_scheme': args.scheme} + + shmem_ok, node_size, shmem_reason = check_shmem() + + if myid == 0: + shmem_status = ('yes' if shmem_ok + else f'NO -- falling back to Bcast ({shmem_reason})') + W = 66 + print('=' * W) + print(' ANUGA distribute benchmark') + print(f' Mesh: rectangular_cross {M}x{M} ({ntri:,} triangles)') + print(f' MPI ranks: {nproc} Repetitions: {args.reps}') + print(f' Ranks/node: {node_size}') + print(f' Scheme: {args.scheme}') + print(f' Shared mem: {shmem_status}') + print('=' * W) + + # Run all four methods + times_std, pss_std, rss_std, ghost_std = time_distribute( + distribute, M, args.reps, args.interval, params) + + gc.collect() + comm.Barrier() + + times_col, pss_col, rss_col, ghost_col = time_distribute( + distribute_collaborative, M, args.reps, args.interval, params) + + gc.collect() + comm.Barrier() + + times_bm, pss_bm, rss_bm, ghost_bm = time_distribute_basic_mesh( + M, args.reps, args.interval, params) + + gc.collect() + comm.Barrier() + + times_dump, times_load, pss_dl, rss_dl, ghost_dl = time_dump_load( + M, args.reps, args.interval, params) + + if not args.no_evolve: + run_evolve_check(M, args.scheme) + + if myid == 0: + times_dl = [d + load_t for d, load_t in zip(times_dump, times_load)] + + W = 32 + C = 22 + + def row(label, *vals): + cells = ' '.join(f'{v:<{C}}' for v in vals) + print(f' {label:<{W}} {cells}') + + sep = '-' * (W + 2 + 4 * (C + 2)) + print(f'\n{sep}') + row('', 'distribute()', 'collaborative()', 'distribute_basic_mesh()', 'dump()+load()') + print(sep) + row('Wall time (median)', + fmt_time(times_std), fmt_time(times_col), + fmt_time(times_bm), fmt_time(times_dl)) + row('Peak PSS sum (physical total)', + fmt_mem(pss_std), fmt_mem(pss_col), + fmt_mem(pss_bm), fmt_mem(pss_dl)) + row('Peak RSS max single rank', + fmt_mem(rss_std), fmt_mem(rss_col), + fmt_mem(rss_bm), fmt_mem(rss_dl)) + print(sep) + print(' Note: PSS sums shared pages proportionally.') + if times_dump: + med_dump = statistics.median(times_dump) + med_load = statistics.median(times_load) + print(f' dump()+load() breakdown: ' + f'dump (serial) {med_dump:.3f}s + ' + f'load (parallel) {med_load:.3f}s') + + # Speedup summary + print() + candidates = [ + ('distribute()', statistics.median(times_std) if times_std else None), + ('collaborative()', statistics.median(times_col) if times_col else None), + ('distribute_basic_mesh()', statistics.median(times_bm) if times_bm else None), + ('dump()+load()', statistics.median(times_dl) if times_dl else None), + ] + valid = [(n, t) for n, t in candidates if t is not None] + best_name, best_time = min(valid, key=lambda x: x[1]) + print(f' Fastest: {best_name} ({best_time:.3f}s)') + for name, t in valid: + if name != best_name and t > 0: + print(f' vs {name}: {t / best_time:.2f}x slower') + + # Ghost triangle stats + print() + print(f' Partition quality ({args.scheme}, {nproc} ranks, {ntri:,} triangles):') + for label, gs in [('distribute()', ghost_std), + ('collaborative()', ghost_col), + ('distribute_basic_mesh()', ghost_bm), + ('dump()+load()', ghost_dl)]: + g_sum, g_max, g_min = gs + if g_sum > 0: + pct = 100.0 * g_sum / ntri + avg = g_sum / nproc + print(f' {label:<28} ghost={g_sum:,} ({pct:.1f}%) ' + f'avg={avg:,.0f}/rank min={g_min:,} max={g_max:,}') + else: + print(f' {label:<28} (no ghost data)') + print() + + +if __name__ == '__main__': + try: + main() + except Exception as e: + import sys + import traceback + print(f'\n[rank {myid}] ERROR: {e}', flush=True) + traceback.print_exc() + comm.Abort(1) + sys.exit(1) + MPI.Finalize() diff --git a/benchmarks/run_benchmark_grid.py b/benchmarks/run_benchmark_grid.py new file mode 100644 index 000000000..cb8ad7730 --- /dev/null +++ b/benchmarks/run_benchmark_grid.py @@ -0,0 +1,325 @@ +#!/usr/bin/env python3 +"""Run distribute_benchmarks.py over a grid of (np, scheme) values and +print a consolidated summary table. + +Each (np, scheme) combination is run as a subprocess and its stdout saved to +/bench_np_.txt. Existing files are reused unless +--force is given, so partial grids can be resumed cheaply. + +Usage +----- + python benchmarks/run_benchmark_grid.py [options] + +Options +------- +--size M Mesh grid size M (4*M*M triangles, default 1000) +--reps R Repetitions per function (default 1) +--interval S Ticker sample interval in seconds (default 5.0) +--np LIST Comma-separated process counts (default 10,20,30) +--schemes LIST Comma-separated schemes (default metis,morton,hilbert) +--outdir DIR Directory for per-run output files (default benchmarks/results/dist) +--mpirun CMD MPI launcher (default mpirun) +--script PATH Path to distribute_benchmarks.py (default: auto-detect) +--force Re-run even if output file already exists +--dry-run Print commands without executing them +--no-summary Skip summary; just run and save +""" + +import argparse +import re +import subprocess +from datetime import datetime +from pathlib import Path + + +# --------------------------------------------------------------------------- +# Argument parsing +# --------------------------------------------------------------------------- + +def parse_args(): + p = argparse.ArgumentParser(add_help=True) + p.add_argument('--size', type=int, default=1000) + p.add_argument('--reps', type=int, default=1) + p.add_argument('--interval', type=float, default=5.0) + p.add_argument('--np', type=str, default='10,20,30') + p.add_argument('--schemes', type=str, default='metis,morton,hilbert') + p.add_argument('--outdir', type=str, default=None) + p.add_argument('--mpirun', type=str, default='mpirun') + p.add_argument('--script', type=str, default=None) + p.add_argument('--force', action='store_true') + p.add_argument('--dry-run', action='store_true') + p.add_argument('--no-summary', action='store_true') + a = p.parse_args() + a.np_values = [int(x) for x in a.np.split(',')] + a.scheme_list = [s.strip() for s in a.schemes.split(',')] + if a.outdir is None: + here = Path(__file__).resolve().parent + a.outdir = str(here / 'results' / 'dist') + return a + + +def find_benchmark_script(explicit=None): + if explicit: + return explicit + here = Path(__file__).resolve().parent + candidate = here / 'distribute_benchmarks.py' + if candidate.exists(): + return str(candidate) + raise FileNotFoundError( + 'Cannot find distribute_benchmarks.py; use --script to specify it.') + + +# --------------------------------------------------------------------------- +# Run one (np, scheme) combination +# --------------------------------------------------------------------------- + +def output_path(outdir, np_val, scheme): + return Path(outdir) / f'bench_np{np_val}_{scheme}.txt' + + +def run_one(mpirun, script, np_val, scheme, size, reps, interval, + outfile, dry_run=False): + cmd = [mpirun, '-np', str(np_val), 'python', script, + '--size', str(size), + '--reps', str(reps), + '--interval', str(interval), + '--scheme', scheme, + '--no-evolve'] + print(f' $ {" ".join(cmd)}') + print(f' -> {outfile}') + if dry_run: + return True + try: + result = subprocess.run( + cmd, capture_output=True, text=True, timeout=3600) + text = result.stdout + if result.returncode != 0: + text += f'\n[STDERR]\n{result.stderr}' + Path(outfile).write_text(text) + if result.returncode != 0: + print(f' [FAILED -- exit code {result.returncode}]') + return False + return True + except subprocess.TimeoutExpired: + print(' [TIMED OUT]') + return False + except Exception as e: + print(f' [ERROR: {e}]') + return False + + +# --------------------------------------------------------------------------- +# Parse one run's stdout into a metrics dict +# --------------------------------------------------------------------------- + +def _first_float(pattern, text): + m = re.search(pattern, text) + return float(m.group(1)) if m else None + + +def _first_int(pattern, text): + m = re.search(pattern, text) + return int(m.group(1).replace(',', '')) if m else None + + +def parse_output(text): + """Extract benchmark metrics from one run's stdout.""" + r = {} + + r['triangles'] = _first_int(r'(\d[\d,]+) triangles', text) + r['ranks'] = _first_int(r'MPI ranks:\s+(\d+)', text) + r['scheme'] = (re.search(r'Scheme:\s+(\S+)', text) or + type('', (), {'group': lambda s, n: 'metis'})()).group(1) + + # Wall time row: 4 columns separated by whitespace + wall = re.search( + r'Wall time[^\d]+([\d.]+)s\s+([\d.]+)s\s+([\d.]+)s\s+([\d.]+)s', text) + if wall: + r['wall_dist'] = float(wall.group(1)) + r['wall_col'] = float(wall.group(2)) + r['wall_bm'] = float(wall.group(3)) + r['wall_dl'] = float(wall.group(4)) + + # PSS row: 4 columns + pss = re.search( + r'Peak PSS[^\d]+([\d.]+) GiB\s+([\d.]+) GiB\s+([\d.]+) GiB\s+([\d.]+) GiB', + text) + if pss: + r['pss_dist'] = float(pss.group(1)) + r['pss_col'] = float(pss.group(2)) + r['pss_bm'] = float(pss.group(3)) + r['pss_dl'] = float(pss.group(4)) + + # RSS row: 4 columns + rss = re.search( + r'Peak RSS[^\d]+([\d.]+) GiB\s+([\d.]+) GiB\s+([\d.]+) GiB\s+([\d.]+) GiB', + text) + if rss: + r['rss_dist'] = float(rss.group(1)) + r['rss_col'] = float(rss.group(2)) + r['rss_bm'] = float(rss.group(3)) + r['rss_dl'] = float(rss.group(4)) + + # dump+load breakdown + dl = re.search(r'dump \(serial\) ([\d.]+)s.*?load \(parallel\) ([\d.]+)s', text) + if dl: + r['dump_time'] = float(dl.group(1)) + r['load_time'] = float(dl.group(2)) + + # Ghost % (from distribute() column) + g_tot = re.search(r'distribute\(\)\s+ghost=[\d,]+\s+\(([\d.]+)%\)', text) + if g_tot: + r['ghost_pct'] = float(g_tot.group(1)) + + return r + + +# --------------------------------------------------------------------------- +# Summary printer +# --------------------------------------------------------------------------- + +def _t(val): + return f'{val:6.2f}s' if val is not None else ' -- ' + + +def _g(val): + return f'{val:6.2f}' if val is not None else ' -- ' + + +def _pct(val): + return f'{val:5.1f}%' if val is not None else ' -- ' + + +def print_summary(results, args, ntri, outdir): + lines = [] + + def out(*a): + s = ' '.join(str(x) for x in a) + lines.append(s) + print(s) + + W = 76 + out('=' * W) + out(' ANUGA distribute() benchmark grid summary') + out(f' Mesh: synthetic {args.size}x{args.size} ({ntri:,} triangles)') + out(f' np values: {args.np_values}') + out(f' Schemes: {args.scheme_list}') + out(f' Reps: {args.reps}') + out(f' Generated: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}') + out('=' * W) + + # Wall time table + out() + out(' Wall time (seconds)') + hdr = (f' {"np":>4} {"scheme":<8} ' + f'{"distribute()":>13} {"collaborative()":>15} ' + f'{"dist_basic_mesh()":>18} {"dump+load":>10} ' + f'{"dump":>7} {"load":>6}') + out(hdr) + out(' ' + '-' * (len(hdr) - 2)) + for np_val in args.np_values: + for scheme in args.scheme_list: + r = results.get((np_val, scheme), {}) + out(f' {np_val:>4} {scheme:<8} ' + f'{_t(r.get("wall_dist")):>13} ' + f'{_t(r.get("wall_col")):>15} ' + f'{_t(r.get("wall_bm")):>18} ' + f'{_t(r.get("wall_dl")):>10} ' + f'{_t(r.get("dump_time")):>7} ' + f'{_t(r.get("load_time")):>6}') + + # Memory table + out() + out(' Peak PSS -- physical memory total across all ranks (GiB)') + hdr2 = (f' {"np":>4} {"scheme":<8} ' + f'{"distribute()":>13} {"collaborative()":>15} ' + f'{"dist_basic_mesh()":>18} {"dump+load":>10}') + out(hdr2) + out(' ' + '-' * (len(hdr2) - 2)) + for np_val in args.np_values: + for scheme in args.scheme_list: + r = results.get((np_val, scheme), {}) + out(f' {np_val:>4} {scheme:<8} ' + f'{_g(r.get("pss_dist")):>13} ' + f'{_g(r.get("pss_col")):>15} ' + f'{_g(r.get("pss_bm")):>18} ' + f'{_g(r.get("pss_dl")):>10}') + + # Ghost triangle table + if any('ghost_pct' in r for r in results.values()): + out() + out(' Partition quality -- ghost triangles (distribute() column)') + out(f' {"np":>4} {"scheme":<8} {"ghost %":>8}') + for np_val in args.np_values: + for scheme in args.scheme_list: + r = results.get((np_val, scheme), {}) + out(f' {np_val:>4} {scheme:<8} {_pct(r.get("ghost_pct")):>8}') + + out() + out('=' * W) + + summary_path = Path(outdir) / 'summary.txt' + summary_path.write_text('\n'.join(lines) + '\n') + print(f'\n Summary saved to {summary_path}') + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + args = parse_args() + script = find_benchmark_script(args.script) + outdir = Path(args.outdir) + outdir.mkdir(parents=True, exist_ok=True) + + print(f'\nBenchmark grid: np={args.np_values} schemes={args.scheme_list}') + print(f'Output directory: {outdir.resolve()}') + print(f'Script: {script}\n') + + ntri = 4 * args.size * args.size + failed = [] + + for np_val in args.np_values: + for scheme in args.scheme_list: + outfile = output_path(outdir, np_val, scheme) + if outfile.exists() and not args.force: + print(f'[SKIP] np={np_val} scheme={scheme} ' + f'(file exists; use --force to re-run)') + continue + print(f'\n[RUN] np={np_val} scheme={scheme}') + ok = run_one(args.mpirun, script, np_val, scheme, + args.size, args.reps, args.interval, + outfile, dry_run=args.dry_run) + if not ok: + failed.append((np_val, scheme)) + + if args.dry_run or args.no_summary: + return + + results = {} + for np_val in args.np_values: + for scheme in args.scheme_list: + outfile = output_path(outdir, np_val, scheme) + if not outfile.exists(): + continue + text = outfile.read_text() + r = parse_output(text) + if r.get('wall_dist') is not None: + results[(np_val, scheme)] = r + else: + print(f' [WARN] Could not parse {outfile}') + + if not results: + print('\nNo parseable results found.') + return + + if failed: + print(f'\nWARNING: {len(failed)} run(s) failed: {failed}') + + print() + print_summary(results, args, ntri, outdir) + + +if __name__ == '__main__': + main() diff --git a/benchmarks/run_benchmarks.py b/benchmarks/run_benchmarks.py new file mode 100644 index 000000000..b4a4b4225 --- /dev/null +++ b/benchmarks/run_benchmarks.py @@ -0,0 +1,345 @@ +#!/usr/bin/env python3 +""" +ANUGA Performance Benchmark Suite +---------------------------------- +Measures wall time and peak memory for simulations of increasing size and +for each multiprocessor_mode (0=Python Euler, 1=Python RK2, 2=C RK2/GPU). + +Results are written to a JSON file for later comparison with +compare_benchmarks.py. + +Usage +----- + # Run all scenarios (small + medium), save to results/ + python benchmarks/run_benchmarks.py + + # Quick sanity check (small only) + python benchmarks/run_benchmarks.py --sizes small + + # All sizes including large (~360K tris, slow) + python benchmarks/run_benchmarks.py --sizes small,medium,large + + # Only mode 0 and 2 + python benchmarks/run_benchmarks.py --modes 0,2 + + # Custom output file + python benchmarks/run_benchmarks.py --output my_results.json + +Metrics +------- +- wall_time_s : wall-clock seconds for the evolve loop only +- n_steps : number of internal (CFL) timesteps taken +- cells_per_s : n_triangles * n_steps / wall_time (the primary figure) +- setup_rss_mb : RSS after domain creation (before any evolve) +- peak_rss_mb : peak RSS sampled during evolve (100 ms polling) +- mb_per_ktri : peak_rss_mb / (n_triangles / 1000) — memory efficiency +""" + +import argparse +import json +import os +import platform +import subprocess +import sys +import tempfile +import threading +import time +from datetime import datetime + + +# --------------------------------------------------------------------------- +# Memory sampler (background thread) +# --------------------------------------------------------------------------- + +class _MemSampler: + """Poll process RSS every `interval_s` seconds in a daemon thread.""" + + def __init__(self, interval_s=0.1): + self._interval = interval_s + self.peak_mb = 0.0 + self._running = False + self._thread = None + + def start(self): + self._running = True + self._thread = threading.Thread(target=self._run, daemon=True) + self._thread.start() + + def stop(self): + self._running = False + if self._thread: + self._thread.join(timeout=2.0) + + def current_mb(self): + try: + import psutil + return psutil.Process(os.getpid()).memory_info().rss / 1024 ** 2 + except ImportError: + pass + # fallback: /proc/self/status on Linux + try: + with open('/proc/self/status') as fh: + for line in fh: + if line.startswith('VmRSS:'): + return int(line.split()[1]) / 1024 + except OSError: + pass + return 0.0 + + def _run(self): + while self._running: + rss = self.current_mb() + if rss > self.peak_mb: + self.peak_mb = rss + time.sleep(self._interval) + + def reset_peak(self): + self.peak_mb = self.current_mb() + + +# --------------------------------------------------------------------------- +# Domain factory +# --------------------------------------------------------------------------- + +def _create_domain(nx, ny, mode, tmpdir): + """ + Create a rectangular dam-break domain with nx*ny cells per side. + + The domain has no file output (store=False) and uses reflective boundaries. + Initial condition: left half stage=2m, right half stage=0.5m. + """ + import numpy as np + import anuga + from anuga import rectangular_cross_domain, Reflective_boundary + + domain = rectangular_cross_domain(nx, ny, len1=1000.0, len2=1000.0) + domain.set_flow_algorithm('DE0') + domain.set_low_froude(0) + domain.set_name('bench') + domain.set_datadir(tmpdir) + domain.store = False + + domain.set_quantity('elevation', 0.0) + domain.set_quantity('stage', lambda x, y: np.where(x < 500.0, 2.0, 0.5)) + domain.set_quantity('xmomentum', 0.0) + domain.set_quantity('ymomentum', 0.0) + domain.set_boundary({t: Reflective_boundary(domain) + for t in domain.get_boundary_tags()}) + + if mode >= 1: + domain.set_multiprocessor_mode(mode) + + return domain + + +# --------------------------------------------------------------------------- +# Scenario definitions +# --------------------------------------------------------------------------- + +SCENARIOS = { + 'small': dict(nx=50, ny=50, finaltime=200.0, yieldstep=50.0), + 'medium': dict(nx=150, ny=150, finaltime=100.0, yieldstep=25.0), + 'large': dict(nx=300, ny=300, finaltime=50.0, yieldstep=12.5), +} + + +# --------------------------------------------------------------------------- +# Run one scenario +# --------------------------------------------------------------------------- + +def run_one(size, mode, omp_threads, sampler): + """ + Run a single benchmark scenario and return a result dict. + + Parameters + ---------- + size : str + One of 'small', 'medium', 'large'. + mode : int + multiprocessor_mode (0, 1, or 2). + omp_threads : int + Value of OMP_NUM_THREADS (informational only — caller must set env). + sampler : _MemSampler + Running memory sampler; peak is reset just before evolve starts. + + Returns + ------- + dict + Benchmark result record. + """ + cfg = SCENARIOS[size] + tmpdir = tempfile.mkdtemp() + + try: + domain = _create_domain(cfg['nx'], cfg['ny'], mode, tmpdir) + n_tris = domain.number_of_triangles + + # Memory snapshot after full domain setup (before any evolve) + setup_rss_mb = sampler.current_mb() + sampler.reset_peak() + + t0 = time.perf_counter() + for _ in domain.evolve(yieldstep=cfg['yieldstep'], + finaltime=cfg['finaltime']): + pass + wall_time_s = time.perf_counter() - t0 + + n_steps = domain.number_of_steps + peak_rss_mb = sampler.peak_mb + cells_per_s = (n_tris * n_steps / wall_time_s) if wall_time_s > 0 else 0.0 + + finally: + import shutil + shutil.rmtree(tmpdir, ignore_errors=True) + + return { + 'name': f'{size}_mode{mode}_t{omp_threads}', + 'size': size, + 'n_triangles': n_tris, + 'mode': mode, + 'omp_threads': omp_threads, + 'finaltime': cfg['finaltime'], + 'n_steps': n_steps, + 'wall_time_s': round(wall_time_s, 3), + 'cells_per_s': round(cells_per_s, 0), + 'setup_rss_mb': round(setup_rss_mb, 1), + 'peak_rss_mb': round(peak_rss_mb, 1), + 'mb_per_ktri': round(peak_rss_mb / (n_tris / 1000), 3) if n_tris else 0.0, + } + + +# --------------------------------------------------------------------------- +# Metadata helpers +# --------------------------------------------------------------------------- + +def _git_info(): + def _run(args): + try: + return subprocess.check_output( + args, cwd=os.path.dirname(__file__), + stderr=subprocess.DEVNULL).decode().strip() + except Exception: + return 'unknown' + + return { + 'commit': _run(['git', 'rev-parse', '--short', 'HEAD']), + 'branch': _run(['git', 'rev-parse', '--abbrev-ref', 'HEAD']), + 'commit_long': _run(['git', 'rev-parse', 'HEAD']), + } + + +def _env_info(): + omp = os.environ.get('OMP_NUM_THREADS', 'unset') + try: + import anuga + anuga_version = anuga.__version__ + except Exception: + anuga_version = 'unknown' + + return { + 'python_version': sys.version.split()[0], + 'platform': platform.system(), + 'hostname': platform.node().split('.')[0], + 'omp_num_threads_env': omp, + 'anuga_version': anuga_version, + } + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + parser = argparse.ArgumentParser( + description='ANUGA performance benchmark suite.', + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=__doc__, + ) + parser.add_argument( + '--sizes', default='small,medium', + help='Comma-separated sizes to run: small,medium,large (default: small,medium)', + ) + parser.add_argument( + '--modes', default='0,1,2', + help='Comma-separated multiprocessor_modes to test (default: 0,1,2)', + ) + parser.add_argument( + '--output', default=None, + help='Output JSON path. Default: benchmarks/results/__.json', + ) + args = parser.parse_args() + + sizes = [s.strip() for s in args.sizes.split(',')] + modes = [int(m.strip()) for m in args.modes.split(',')] + + for s in sizes: + if s not in SCENARIOS: + parser.error(f'Unknown size {s!r}. Choose from: {list(SCENARIOS)}') + + # --- output path --- + git = _git_info() + env = _env_info() + timestamp = datetime.now().strftime('%Y%m%d_%H%M%S') + omp_threads = int(os.environ.get('OMP_NUM_THREADS', 1)) + + if args.output: + outpath = args.output + else: + outdir = os.path.join(os.path.dirname(__file__), 'results') + os.makedirs(outdir, exist_ok=True) + fname = f"{git['branch'].replace('/', '_')}_{git['commit']}_{timestamp}.json" + outpath = os.path.join(outdir, fname) + + print(f'ANUGA benchmark commit={git["commit"]} branch={git["branch"]}') + print(f'Python {env["python_version"]} OMP_NUM_THREADS={omp_threads}') + print(f'Sizes: {sizes} Modes: {modes}') + print(f'Output: {outpath}') + print() + + sampler = _MemSampler(interval_s=0.1) + sampler.start() + + results = [] + header = f"{'Scenario':<28} {'tris':>8} {'mode':>4} {'thrd':>4} {'steps':>6} {'wall(s)':>8} {'cells/s':>10} {'setup MB':>9} {'peak MB':>8} {'MB/Ktri':>8}" + rule = '-' * len(header) + print(header) + print(rule) + + for size in sizes: + for mode in modes: + name = f'{size}_mode{mode}_t{omp_threads}' + sys.stdout.write(f' Running {name} ... ') + sys.stdout.flush() + try: + rec = run_one(size, mode, omp_threads, sampler) + results.append(rec) + print( + f"\r {rec['name']:<28} {rec['n_triangles']:>8} {rec['mode']:>4} " + f"{rec['omp_threads']:>4} {rec['n_steps']:>6} " + f"{rec['wall_time_s']:>8.2f} {rec['cells_per_s']:>10,.0f} " + f"{rec['setup_rss_mb']:>9.1f} {rec['peak_rss_mb']:>8.1f} " + f"{rec['mb_per_ktri']:>8.3f}" + ) + except Exception as exc: + print(f'\r {name:<28} FAILED: {exc}') + + sampler.stop() + print(rule) + print() + + payload = { + 'timestamp': timestamp, + 'git': git, + 'env': env, + 'omp_threads': omp_threads, + 'results': results, + } + with open(outpath, 'w') as fh: + json.dump(payload, fh, indent=2) + print(f'Saved {len(results)} results → {outpath}') + + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/docs/source/setup_anuga_script/index.rst b/docs/source/setup_anuga_script/index.rst index 9d65e5a82..5ba953d77 100644 --- a/docs/source/setup_anuga_script/index.rst +++ b/docs/source/setup_anuga_script/index.rst @@ -18,6 +18,7 @@ Script Structure operators evolve checkpointing + logging .. only:: html diff --git a/docs/source/setup_anuga_script/logging.rst b/docs/source/setup_anuga_script/logging.rst new file mode 100644 index 000000000..a11342cff --- /dev/null +++ b/docs/source/setup_anuga_script/logging.rst @@ -0,0 +1,156 @@ + +.. currentmodule:: anuga + +Logging +======= + +By default ANUGA prints progress messages to the terminal (stdout). +Calling :func:`set_logfile` activates file logging: from that point +every ``print()`` call and every ``log.*()`` call is written to both the +terminal **and** the named log file simultaneously. + +Quick start +----------- + +:: + + import anuga + + anuga.set_logfile('my_run.log') + + # All subsequent output goes to terminal AND my_run.log + print('Setting up domain ...') + + domain = anuga.rectangular_cross_domain(10, 5) + ... + + for t in domain.evolve(yieldstep=1.0, finaltime=10.0): + domain.print_timestepping_statistics() + +After this call the file ``my_run.log`` contains a complete record of the +run including timestep statistics, warnings, and any other printed output. + + +Output destinations +------------------- + +ANUGA's output is split into three categories: + +.. list-table:: + :header-rows: 1 + :widths: 20 15 15 50 + + * - Call + - Terminal + - Log file + - When to use + * - ``print(msg)`` + - ✓ + - ✓ + - User script milestones and results + * - ``log.info(msg)`` + - ✓ + - ✓ + - Significant simulation events + * - ``log.warning(msg)`` + - ✓ + - ✓ + - Non-fatal issues that need attention + * - ``log.verbose(msg)`` + - ✗ + - ✓ + - Internal ANUGA chatter (mesh stats, solver steps) + * - ``log.debug(msg)`` + - ✗ + - ✓ + - Developer diagnostics + +``log.verbose()`` and ``log.debug()`` are below the default console +threshold (``INFO``) so they are silently dropped when no log file is +active. + + +Capturing verbose third-party output +------------------------------------- + +Some ANUGA functions (mesh construction, domain distribution) produce +detailed ``print()`` output when called with ``verbose=True``. Use the +:func:`file_only` context manager to redirect that output to the log +file without showing it on the terminal:: + + import anuga + import anuga.utilities.log as log + + anuga.set_logfile('my_run.log') + + with log.file_only(): + domain = anuga.create_domain_from_regions( + bounding_polygon, + boundary_tags=tags, + maximum_triangle_area=res, + verbose=True, # full output goes to file, not screen + ) + +Outside the ``with`` block normal tee behaviour resumes immediately. + + +Showing verbose output on the terminal +--------------------------------------- + +When debugging it can be useful to see all output on the terminal as +well. Pass ``verbose_to_screen=True`` to :func:`set_logfile`:: + + anuga.set_logfile('debug.log', verbose_to_screen=True) + +This lowers the console threshold from ``INFO`` to ``DEBUG`` so +``log.verbose()`` and ``log.debug()`` messages also appear on screen. + + +Using log levels directly +-------------------------- + +The ``anuga.utilities.log`` module is available as ``anuga.log`` and +exposes the standard Python logging levels:: + + import anuga.utilities.log as log + + log.info('Mesh built successfully') + log.warning('Elevation data outside domain extent — using zero') + log.verbose('Triangle count: %d' % n) # file only + log.debug('CG solver iteration %d, residual %.2e' % (k, r)) + +Levels in order of decreasing severity: + +.. list-table:: + :header-rows: 1 + :widths: 20 80 + + * - Level + - Constant + * - ``CRITICAL`` + - ``log.CRITICAL`` (50) — fatal errors + * - ``ERROR`` + - ``log.ERROR`` (40) — recoverable errors + * - ``WARNING`` + - ``log.WARNING`` (30) — non-fatal issues + * - ``INFO`` + - ``log.INFO`` (20) — normal progress (default console threshold) + * - ``DEBUG`` / ``VERBOSE`` + - ``log.DEBUG`` (10) — detail (default file threshold) + + +API reference +------------- + +.. autofunction:: anuga.set_logfile + +.. autoclass:: anuga.utilities.log.TeeStream + :members: write, flush, close + +.. autofunction:: anuga.utilities.log.file_only + +.. autofunction:: anuga.utilities.log.verbose +.. autofunction:: anuga.utilities.log.info +.. autofunction:: anuga.utilities.log.warning +.. autofunction:: anuga.utilities.log.debug +.. autofunction:: anuga.utilities.log.critical diff --git a/examples/cairns_toml_excel/run_model.py b/examples/cairns_toml_excel/run_model.py index ef572e0fa..b1a0cfaa2 100644 --- a/examples/cairns_toml_excel/run_model.py +++ b/examples/cairns_toml_excel/run_model.py @@ -141,21 +141,12 @@ def progress(msg): print('Evolving') -import io -_logfile = sys.stdout # Logger or normal stdout barrier() for t in domain.evolve(yieldstep=project.yieldstep, finaltime=project.finaltime, outputstep=project.outputstep): if myid == 0: - buf = io.StringIO() - sys.stdout = buf domain.print_timestepping_statistics() - sys.stdout = _logfile - stats = buf.getvalue() - sys.__stdout__.write(stats) - sys.__stdout__.flush() - _logfile.write(stats) if project.report_mass_conservation_statistics: domain.report_water_volume_statistics() diff --git a/examples/cairns_toml_excel/setup/prepare_data.py b/examples/cairns_toml_excel/setup/prepare_data.py index 7b2a54786..109934f0f 100644 --- a/examples/cairns_toml_excel/setup/prepare_data.py +++ b/examples/cairns_toml_excel/setup/prepare_data.py @@ -17,6 +17,7 @@ import numpy import anuga from anuga.utilities import spatialInputUtil as su +from anuga.utilities.log import TeeStream from anuga.parallel import myid, barrier, send, receive, numprocs # Local modules @@ -25,23 +26,6 @@ from setup.parse_input_data import ProjectData -class Logger(object): - - """Makes it simple to get stdout to flush, and seems to work in - parallel - """ - - def __init__(self, logfile): - self.log = open(logfile, "a") - - def write(self, message): - self.log.write(message) - self.log.flush() - - def flush(self): - self.log.flush() - - class PrepareData(ProjectData): """Converts the information in the configuration file to inputs @@ -104,7 +88,7 @@ def define_output_directory_and_redirect_stdout(self, print('OUTPUT_DIRECTORY: ' + str(self.output_dir)) - # Send stdout to a file inside the output directory + # Tee stdout to a file inside the output directory if output_log is not None: if make_directories: stdout_file = self.output_dir + '/' + output_log @@ -113,11 +97,7 @@ def define_output_directory_and_redirect_stdout(self, if myid == 0: print('Redirecting output now to ' + stdout_file) - sys.stdout = Logger(stdout_file) - barrier() - - if myid != 0: - sys.stdout = Logger(stdout_file) + sys.stdout = TeeStream(stdout_file) return diff --git a/pyproject.toml b/pyproject.toml index b3b1c6d1e..cbb67221c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,7 +28,7 @@ dependencies = [ 'netCDF4>=1.6', 'scipy>=1.11', 'meshpy>=2022.1', - 'utm', + 'pyproj>=3.6', 'xarray', ] classifiers = [ @@ -56,8 +56,6 @@ parallel = [ 'pymetis>=2023.1', ] data = [ - 'pyproj>=3.6', - 'affine>=2.4', 'rasterio', 'fiona', 'shapely', diff --git a/scripts/anuga_run_toml.py b/scripts/anuga_run_toml.py index c525aec6f..8ce3d65cb 100755 --- a/scripts/anuga_run_toml.py +++ b/scripts/anuga_run_toml.py @@ -149,21 +149,12 @@ def progress(msg): progress('Evolving') -import io -_logfile = sys.stdout barrier() for t in domain.evolve(yieldstep=project.yieldstep, finaltime=project.finaltime, outputstep=project.outputstep): if myid == 0: - buf = io.StringIO() - sys.stdout = buf domain.print_timestepping_statistics() - sys.stdout = _logfile - stats = buf.getvalue() - sys.__stdout__.write(stats) - sys.__stdout__.flush() - _logfile.write(stats) if project.report_mass_conservation_statistics: domain.report_water_volume_statistics()