From 37473a40a924e545e9e7d81df5eabd7790a98df6 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sun, 5 Apr 2026 11:47:25 +1000 Subject: [PATCH 01/28] L1-L4 logging refactor: TeeStream, lazy log file, set_logfile() MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit L1: Replace Logger (file-only) with TeeStream in prepare_data.py so print() goes to both terminal and file. Remove the StringIO hack from anuga_run_toml.py and run_model.py — plain print() now works. L2: Change DefaultConsoleLogLevel from CRITICAL to INFO so log.info() is visible on the terminal by default. Remove the clamping rule that prevented file from logging more than console. L3: Default log_filename=None (lazy). No file is created on import; file logging only starts when set_logfile() is called. Eliminates spurious timestamped log files in test directories. L4: New set_logfile(path) function installs TeeStream on sys.stdout and configures the logging.FileHandler simultaneously — one call makes all print() and log.*() output go to both terminal and file. Exported as anuga.set_logfile() and anuga.TeeStream in public API. Co-Authored-By: Claude Sonnet 4.6 --- anuga/__init__.py | 3 + anuga/scenario/prepare_data.py | 26 +- anuga/utilities/log.py | 364 +++++++++--------- examples/cairns_toml_excel/run_model.py | 9 - .../cairns_toml_excel/setup/prepare_data.py | 26 +- scripts/anuga_run_toml.py | 9 - 6 files changed, 183 insertions(+), 254 deletions(-) diff --git a/anuga/__init__.py b/anuga/__init__.py index c8b7dde12..cda2cc4bb 100644 --- a/anuga/__init__.py +++ b/anuga/__init__.py @@ -323,6 +323,7 @@ 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 from anuga.config import g from anuga.config import velocity_protection @@ -491,6 +492,8 @@ def get_args(): 'memory_stats', 'print_memory_stats', 'log', + 'set_logfile', + 'TeeStream', 'parse_standard_args', 'parse_time', 'plot_utils', 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/utilities/log.py b/anuga/utilities/log.py index dc75e18df..2a35695db 100644 --- a/anuga/utilities/log.py +++ b/anuga/utilities/log.py @@ -1,30 +1,28 @@ #!/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 @@ -33,51 +31,112 @@ import logging -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) -log_filename = os.path.join('.', f'anuga_{current_datetime()}.log') + +################################################################################ +# 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 +################################################################################ +def set_logfile(path, + console_level=DefaultConsoleLogLevel, + file_level=DefaultFileLogLevel): + """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() / log.debug() etc. also write to *path* via the Python + logging module. + - 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). + file_level : int + Logging level for file output (default 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 +149,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 +213,35 @@ 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).''' - log(msg, logging.DEBUG) - def info(msg=''): - '''Shortcut for log(INFO, msg).''' - log(msg, logging.INFO) - def warning(msg=''): - '''Shortcut for log(WARNING, msg).''' - log(msg, logging.WARNING) - def error(msg=''): - '''Shortcut for log(ERROR, msg).''' - log(msg, logging.ERROR) - def critical(msg=''): - '''Shortcut for log(CRITICAL, msg).''' - log(msg, logging.CRITICAL) def timingInfo(msg=''): - '''Shortcut for log(timingDelimiter, msg).''' - log(TimingDelimiter + msg, logging.INFO) @@ -228,57 +255,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 +302,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 +327,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 +337,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 +383,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/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/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() From efaf1f0c131cfb2e6f8bd516ae07053cd111c616 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sun, 5 Apr 2026 12:02:45 +1000 Subject: [PATCH 02/28] Add log.verbose() and verbose_to_screen flag; quiet mesh construction output MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit log.verbose(msg): file-only output for internal ANUGA chatter. Same level as log.debug() but named to guide contributors: use this instead of print() inside verbose=True code paths. set_logfile(..., verbose_to_screen=True): drops console threshold to DEBUG so verbose output also appears on screen — useful for debugging. setup_mesh.py: replace verbose print() calls with log.verbose() / log.info() and pass verbose=False to create_pmesh_from_regions and sequential_distribute_* so mesh construction no longer floods the terminal. Key stats (triangle count, extent) still appear on screen via log.info(). Co-Authored-By: Claude Sonnet 4.6 --- anuga/__init__.py | 1 + anuga/scenario/setup_mesh.py | 63 +++++++++++++++--------------------- anuga/utilities/log.py | 29 ++++++++++++++--- 3 files changed, 52 insertions(+), 41 deletions(-) diff --git a/anuga/__init__.py b/anuga/__init__.py index cda2cc4bb..10d15de68 100644 --- a/anuga/__init__.py +++ b/anuga/__init__.py @@ -324,6 +324,7 @@ def get_args(): from anuga.utilities import log as log from anuga.utilities.log import set_logfile, TeeStream +from anuga.utilities.log import verbose as log_verbose from anuga.config import g from anuga.config import velocity_protection diff --git a/anuga/scenario/setup_mesh.py b/anuga/scenario/setup_mesh.py index a234390ce..4bc564ad2 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,8 +40,7 @@ def build_mesh(project): su.combine_breakLines_and_riverWalls_for_mesh(project.breaklines, project.riverwalls) - # Make the mesh - + # Make the mesh — verbose output goes to file only via log.verbose() anuga.create_pmesh_from_regions( project.bounding_polygon, boundary_tags=project.boundary_tags, @@ -49,7 +48,7 @@ def build_mesh(project): filename=project.meshname, interior_regions=project.interior_regions, use_cache=False, - verbose=verbose, + verbose=False, breaklines=mesh_breaklines, regionPtArea=project.region_point_areas, ) @@ -58,26 +57,25 @@ def build_mesh(project): 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 +110,55 @@ 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') + log.verbose('Saving domain') sequential_distribute_dump(domain, 1, partition_dir=project.partition_dir, - verbose=verbose) + verbose=False) 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') + log.verbose('Load in saved sequential pickled domain') domain = sequential_distribute_load_pickle_file( - pickle_name, np=1, verbose=verbose) + pickle_name, np=1, verbose=False) - if verbose: - print('Dump partitioned domains') + log.verbose('Dump partitioned domains') sequential_distribute_dump( domain, numprocs, - partition_dir=project.partition_dir, verbose=verbose) + partition_dir=project.partition_dir, verbose=False) 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) + verbose=False) # ######################################################################### # Set output directories diff --git a/anuga/utilities/log.py b/anuga/utilities/log.py index 2a35695db..974c033b5 100644 --- a/anuga/utilities/log.py +++ b/anuga/utilities/log.py @@ -100,16 +100,21 @@ def __getattr__(self, name): # 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): + 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() / log.debug() etc. also write to *path* via the Python - logging module. + - 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 @@ -118,9 +123,17 @@ def set_logfile(path, 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). + 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 @@ -226,6 +239,14 @@ def log_exception_hook(type, value, tb): # Shortcut routines ################################################################################ +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(msg, logging.DEBUG) From d251a0955adaf6f950401c6c9b8877895f626441 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sun, 5 Apr 2026 12:14:53 +1000 Subject: [PATCH 03/28] Add log.file_only() context manager; capture mesh verbose output to file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit _FileOnlyStream: file-only counterpart to TeeStream — writes to the log file handle without touching the terminal. file_only(): context manager that temporarily installs _FileOnlyStream as sys.stdout so any print() inside the block goes to the log file only. Restores original stdout on exit. If no logfile is active, output is discarded (StringIO). Exported as anuga.file_only(). setup_mesh.py: wrap create_pmesh_from_regions and sequential_distribute_* calls with file_only() and restore verbose=True so the full mesh construction output is captured in the log file but not shown on screen. Co-Authored-By: Claude Sonnet 4.6 --- anuga/__init__.py | 2 +- anuga/scenario/setup_mesh.py | 51 ++++++++++++++++++++---------------- anuga/utilities/log.py | 45 +++++++++++++++++++++++++++++++ 3 files changed, 74 insertions(+), 24 deletions(-) diff --git a/anuga/__init__.py b/anuga/__init__.py index 10d15de68..b07070e35 100644 --- a/anuga/__init__.py +++ b/anuga/__init__.py @@ -323,7 +323,7 @@ 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 +from anuga.utilities.log import set_logfile, TeeStream, file_only from anuga.utilities.log import verbose as log_verbose from anuga.config import g diff --git a/anuga/scenario/setup_mesh.py b/anuga/scenario/setup_mesh.py index 4bc564ad2..51229260e 100644 --- a/anuga/scenario/setup_mesh.py +++ b/anuga/scenario/setup_mesh.py @@ -40,18 +40,19 @@ def build_mesh(project): su.combine_breakLines_and_riverWalls_for_mesh(project.breaklines, project.riverwalls) - # Make the mesh — verbose output goes to file only via log.verbose() - 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=False, - 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 @@ -126,9 +127,10 @@ def setup_mesh(project, setup_initial_conditions=None): domain, project) log.verbose('Saving domain') - sequential_distribute_dump(domain, 1, - partition_dir=project.partition_dir, - verbose=False) + 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) @@ -137,13 +139,15 @@ def setup_mesh(project, setup_initial_conditions=None): log.verbose('Saved partitioned domain seems to already exist') else: log.verbose('Load in saved sequential pickled domain') - domain = sequential_distribute_load_pickle_file( - pickle_name, np=1, verbose=False) + with log.file_only(): + domain = sequential_distribute_load_pickle_file( + pickle_name, np=1, verbose=True) log.verbose('Dump partitioned domains') - sequential_distribute_dump( - domain, numprocs, - partition_dir=project.partition_dir, verbose=False) + with log.file_only(): + sequential_distribute_dump( + domain, numprocs, + partition_dir=project.partition_dir, verbose=True) domain = None gc.collect() @@ -156,9 +160,10 @@ def setup_mesh(project, setup_initial_conditions=None): log.info('Loading partitioned domain') - domain = sequential_distribute_load( - filename=join(project.partition_dir, 'domain'), - verbose=False) + 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 974c033b5..b75640f02 100644 --- a/anuga/utilities/log.py +++ b/anuga/utilities/log.py @@ -29,6 +29,7 @@ import sys import traceback import logging +from contextlib import contextmanager DefaultConsoleLogLevel = logging.INFO @@ -71,6 +72,50 @@ 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() + + 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. ################################################################################ From a4a434a9afccf04031455d5a64a605b347556467 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sun, 5 Apr 2026 12:26:22 +1000 Subject: [PATCH 04/28] Add logging documentation page and docstrings New docs/source/setup_anuga_script/logging.rst covering: - Quick start with set_logfile() - Output destination table (print/info/verbose/debug) - file_only() context manager for capturing verbose internals - verbose_to_screen flag for debugging - Log level reference table - API autofunction blocks Add short docstrings to log.py shortcut functions so autofunction has content to render. Co-Authored-By: Claude Sonnet 4.6 --- anuga/utilities/log.py | 5 + docs/source/setup_anuga_script/index.rst | 1 + docs/source/setup_anuga_script/logging.rst | 156 +++++++++++++++++++++ 3 files changed, 162 insertions(+) create mode 100644 docs/source/setup_anuga_script/logging.rst diff --git a/anuga/utilities/log.py b/anuga/utilities/log.py index b75640f02..7aea8bf0c 100644 --- a/anuga/utilities/log.py +++ b/anuga/utilities/log.py @@ -293,18 +293,23 @@ def verbose(msg=''): log(msg, logging.DEBUG) def debug(msg=''): + """Log a DEBUG-level message (file only by default).""" log(msg, logging.DEBUG) def info(msg=''): + """Log an INFO-level message (terminal and file).""" log(msg, logging.INFO) def warning(msg=''): + """Log a WARNING-level message (terminal and file).""" log(msg, logging.WARNING) def error(msg=''): + """Log an ERROR-level message (terminal and file).""" log(msg, logging.ERROR) def critical(msg=''): + """Log a CRITICAL-level message (terminal and file).""" log(msg, logging.CRITICAL) def timingInfo(msg=''): 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 From 246dbca0cb70f3c011afd6366b0d26f748dae3a5 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sun, 5 Apr 2026 12:31:37 +1000 Subject: [PATCH 05/28] Fix pyproj DeprecationWarning in tif2point_values for single-point queries When indices=[N] selects a single centroid, points[:,0] is a 1-element ndarray. pyproj dispatches to _transform_point (scalar path) which fails with numpy >= 1.25 when given a non-0d array. Fix: pass .tolist() to transformer.transform() so pyproj always takes the array code path, then convert the results back to numpy arrays for the affine transform multiplication. Co-Authored-By: Claude Sonnet 4.6 --- anuga/file_conversion/tif2point_values.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/anuga/file_conversion/tif2point_values.py b/anuga/file_conversion/tif2point_values.py index be8526dd8..421d27afe 100644 --- a/anuga/file_conversion/tif2point_values.py +++ b/anuga/file_conversion/tif2point_values.py @@ -45,7 +45,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)) From cb1f1eb38b9e268cf0d0d2de76e3898a4ec9f076 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sun, 5 Apr 2026 12:40:10 +1000 Subject: [PATCH 06/28] Archive CuPy/CUDA files out of anuga/shallow_water/ Move sw_domain_cupy.py, sw_domain_cuda.py, cuda_anuga.cu and cupy_tests/ to archive/cupy_cuda/ to reduce clutter. These are not wired into meson or the public API; the active GPU work is the OpenMP-offloading path in sw_domain_openmp_ext.pyx. A README.md explains the contents and archive date. Co-Authored-By: Claude Sonnet 4.6 --- archive/cupy_cuda/README.md | 22 +++++++++++++++++++ .../cupy_cuda}/cuda_anuga.cu | 0 .../cupy_cuda}/sw_domain_cuda.py | 0 .../cupy_cuda}/sw_domain_cupy.py | 0 .../cupy_cuda/tests}/run_cupy_cft.py | 0 .../tests}/run_cupy_extrapolation.py | 0 .../tests}/run_cupy_flux_distribute.py | 0 .../tests}/run_cupy_protect_negative.py | 0 .../run_cupy_update_conserved_quantities.py | 0 .../cupy_cuda/tests}/run_sw_domain_cupy.py | 0 .../cupy_cuda/tests}/test_sw_domain_cupy.py | 0 11 files changed, 22 insertions(+) create mode 100644 archive/cupy_cuda/README.md rename {anuga/shallow_water => archive/cupy_cuda}/cuda_anuga.cu (100%) rename {anuga/shallow_water => archive/cupy_cuda}/sw_domain_cuda.py (100%) rename {anuga/shallow_water => archive/cupy_cuda}/sw_domain_cupy.py (100%) rename {anuga/shallow_water/cupy_tests => archive/cupy_cuda/tests}/run_cupy_cft.py (100%) rename {anuga/shallow_water/cupy_tests => archive/cupy_cuda/tests}/run_cupy_extrapolation.py (100%) rename {anuga/shallow_water/cupy_tests => archive/cupy_cuda/tests}/run_cupy_flux_distribute.py (100%) rename {anuga/shallow_water/cupy_tests => archive/cupy_cuda/tests}/run_cupy_protect_negative.py (100%) rename {anuga/shallow_water/cupy_tests => archive/cupy_cuda/tests}/run_cupy_update_conserved_quantities.py (100%) rename {anuga/shallow_water/cupy_tests => archive/cupy_cuda/tests}/run_sw_domain_cupy.py (100%) rename {anuga/shallow_water/cupy_tests => archive/cupy_cuda/tests}/test_sw_domain_cupy.py (100%) 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 From 9b12820c142645a5f1978193845f02665e2d852f Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sun, 5 Apr 2026 12:49:34 +1000 Subject: [PATCH 07/28] Fix test_sww2csv_multiple_files: run in tempdir to avoid stale SWW glob get_all_swwfiles() globs for 'datatest1*.sww' in CWD. Leftover datatest1.sww files from prior runs were being picked up, producing more rows than the 4 expected and causing IndexError. Fix: chdir into a fresh tempdir for the test body, set domain datadir to '.', and clean up the whole tmpdir on exit. Nullify self.sww before restoring CWD so tearDown doesn't try to remove an already-deleted file. Co-Authored-By: Claude Sonnet 4.6 --- .../tests/test_gauge.py | 58 +++++++------------ 1 file changed, 21 insertions(+), 37 deletions(-) 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 - #------------------------------------------------------------- From cd6bb2116430d5d5bcf05499054abb2b5639e56a Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Tue, 7 Apr 2026 18:25:35 +1000 Subject: [PATCH 08/28] Add performance benchmark suite (run_benchmarks.py + compare_benchmarks.py) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Measures wall time and peak RSS for small/medium/large dam-break scenarios across multiprocessor_mode 0/1/2. Results saved as JSON for comparison across commits. compare_benchmarks.py shows ±% deltas for speed and memory. benchmarks/results/ is gitignored. Co-Authored-By: Claude Sonnet 4.6 --- .gitignore | 1 + benchmarks/README.md | 74 +++++++ benchmarks/compare_benchmarks.py | 212 +++++++++++++++++++ benchmarks/run_benchmarks.py | 345 +++++++++++++++++++++++++++++++ 4 files changed, 632 insertions(+) create mode 100644 benchmarks/README.md create mode 100644 benchmarks/compare_benchmarks.py create mode 100644 benchmarks/run_benchmarks.py 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/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 000000000..1df721fb2 --- /dev/null +++ b/benchmarks/README.md @@ -0,0 +1,74 @@ +# 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. 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/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()) From d76cf3e111f5354f1f4495081a12b58c56e23380 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Tue, 7 Apr 2026 19:27:19 +1000 Subject: [PATCH 09/28] Add distribute_benchmarks.py: unified MPI distribution benchmark (4 methods) Merges scripts/benchmark_distribute.py and scripts/benchmark_distribute_mesh.py into benchmarks/distribute_benchmarks.py covering all four approaches: distribute(), collaborative(), distribute_basic_mesh(), dump()+load() Also moves run_benchmark_grid.py to benchmarks/ with updated paths and 4-column table parsing for the new distribute_basic_mesh() column. Co-Authored-By: Claude Sonnet 4.6 --- benchmarks/README.md | 38 ++ benchmarks/distribute_benchmarks.py | 574 ++++++++++++++++++++++++++++ benchmarks/run_benchmark_grid.py | 325 ++++++++++++++++ 3 files changed, 937 insertions(+) create mode 100644 benchmarks/distribute_benchmarks.py create mode 100644 benchmarks/run_benchmark_grid.py diff --git a/benchmarks/README.md b/benchmarks/README.md index 1df721fb2..8529a8d7d 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -72,3 +72,41 @@ python benchmarks/compare_benchmarks.py /tmp/before.json /tmp/after.json --thres 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/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() From b01a5e6314bc9b62cf89253af00df8bf42d01ee4 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Tue, 7 Apr 2026 19:52:55 +1000 Subject: [PATCH 10/28] Fix Basic_mesh.reorder() producing stale neighbours after triangle reorder MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When _neighbours had not yet been accessed (None), reorder() skipped the index-remapping step. Later, accessing mesh.neighbours triggered _build_neighbours() which copied _triangle_neighbours — a pre-reorder array whose indices no longer matched the new triangle numbering. This caused distribute_basic_mesh() to compute ~59% more ghost triangles than distribute() for the same mesh and scheme (measured: 9,002 vs 5,655 on a 1M-triangle metis/4-rank run). Fix: trigger _build_neighbours() at the start of reorder() so _neighbours is always populated and correctly remapped. Adds a regression test that verifies each neighbour pair shares exactly 2 nodes after a random reordering of a Basic_mesh. Co-Authored-By: Claude Sonnet 4.6 --- .../abstract_2d_finite_volumes/basic_mesh.py | 8 ++++ .../tests/test_neighbour_mesh_reorder.py | 48 +++++++++++++++++++ 2 files changed, 56 insertions(+) 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/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 From afa74dabc55e8379405da4dda8e2405a1bfae24b Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 27 Apr 2026 19:58:52 +0000 Subject: [PATCH 11/28] Bump conda-incubator/setup-miniconda in the github-actions group Bumps the github-actions group with 1 update: [conda-incubator/setup-miniconda](https://github.com/conda-incubator/setup-miniconda). Updates `conda-incubator/setup-miniconda` from 3 to 4 - [Release notes](https://github.com/conda-incubator/setup-miniconda/releases) - [Changelog](https://github.com/conda-incubator/setup-miniconda/blob/main/CHANGELOG.md) - [Commits](https://github.com/conda-incubator/setup-miniconda/compare/v3...v4) --- updated-dependencies: - dependency-name: conda-incubator/setup-miniconda dependency-version: '4' dependency-type: direct:production update-type: version-update:semver-major dependency-group: github-actions ... Signed-off-by: dependabot[bot] --- .github/workflows/conda-setup.yml | 2 +- .github/workflows/python-publish-pypi.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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..533f356a6 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -31,7 +31,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 From 5ba1c3c608fae4386a4b0a12070c92ff03cee7ed Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Mon, 11 May 2026 19:47:25 +1000 Subject: [PATCH 12/28] fix(ci): add pip to conda create in PyPI publish workflow Bare `conda create python=X` does not include pip, causing `python -m pip install build` to fail with 'No module named pip'. Co-Authored-By: Claude Sonnet 4.6 --- .github/workflows/python-publish-pypi.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 533f356a6..1fffa5f80 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -41,7 +41,7 @@ jobs: - name: Install Our own environment 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 - name: Install gcc compilers on Windows if: runner.os == 'Windows' From a8368fa233199e4aca0ac430022faee3800dd10c Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Mon, 11 May 2026 20:02:20 +1000 Subject: [PATCH 13/28] fix(ci): use --no-isolation for wheel builds in PyPI workflow Build isolation creates a fresh venv that picks up the system pkg-config, which resolves to Ubuntu's system Python headers. On ubuntu-latest these headers require a split pyconfig.h that is not present, causing the Cython sanity check to fail with 'cannot compile programs'. Pre-install all build deps (meson-python, meson, ninja, cython, pybind11, numpy) into the conda env and pass --no-isolation to python -m build, so meson finds the conda env's Python headers. Mirrors the approach used in the conda-setup.yml test workflow. Also fix runner.os == 'linux' -> 'Linux' so the sdist step actually runs. Co-Authored-By: Claude Sonnet 4.6 --- .github/workflows/python-publish-pypi.yml | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 1fffa5f80..9129c1444 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -41,7 +41,8 @@ jobs: - name: Install Our own environment shell: bash -el {0} run: | - conda create -n anuga_env -c conda-forge python=${{matrix.python-version}} pip + conda create -n anuga_env -c conda-forge python=${{matrix.python-version}} pip \ + meson-python meson ninja cython pybind11 numpy - name: Install gcc compilers on Windows if: runner.os == 'Windows' @@ -69,8 +70,8 @@ jobs: shell: bash -el {0} run: | conda activate anuga_env - python -m pip install build - python -m build --wheel --outdir dist-wheel + pip install build + python -m build --wheel --no-isolation --outdir dist-wheel - name: Repair wheels shell: bash -el {0} @@ -80,11 +81,11 @@ jobs: repairwheel -o dist dist-wheel/*.whl - name: Build source distribution - if: runner.os == 'linux' && matrix.python-version == '3.13' + if: runner.os == 'Linux' && matrix.python-version == '3.13' shell: bash -el {0} run: | conda activate anuga_env - python -m build --sdist --outdir dist + python -m build --sdist --no-isolation --outdir dist - name: Upload distributions uses: actions/upload-artifact@v7 From 89ededfc17ae64e47b4303eee4b8455e7a6b6851 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Mon, 11 May 2026 20:18:35 +1000 Subject: [PATCH 14/28] fix(ci): use pip wheel instead of python -m build for wheel builds conda's ninja package installs the binary but does not register a Python distribution, so python -m build --no-isolation fails its package check with 'Missing dependencies: ninja'. pip wheel --no-build-isolation uses whatever is on PATH (including conda-installed meson/ninja/cython) without a Python-package scan. Co-Authored-By: Claude Sonnet 4.6 --- .github/workflows/python-publish-pypi.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 9129c1444..719598c17 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -70,8 +70,7 @@ jobs: shell: bash -el {0} run: | conda activate anuga_env - pip install build - python -m build --wheel --no-isolation --outdir dist-wheel + pip wheel --no-build-isolation --no-deps -w dist-wheel . - name: Repair wheels shell: bash -el {0} @@ -85,6 +84,7 @@ jobs: shell: bash -el {0} run: | conda activate anuga_env + pip install build python -m build --sdist --no-isolation --outdir dist - name: Upload distributions From c4391e289e17621b1192a05e95499fdd015926a3 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Mon, 11 May 2026 20:22:55 +1000 Subject: [PATCH 15/28] =?UTF-8?q?refactor(ci):=20restructure=20PyPI=20work?= =?UTF-8?q?flow=20=E2=80=94=20build=20on=20push,=20publish=20on=20release?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Previously the workflow only triggered on release:published, so there was no way to verify wheels built successfully before tagging. New structure: - build-wheels: runs on every push/PR to main/develop (early feedback) and on release events. Builds one wheel per (python, os) combination using pip wheel --no-build-isolation to avoid the conda-ninja/pkg-config mismatch. Repairs Linux wheels with repairwheel. - build-sdist: builds source distribution once (Linux/py3.13). - publish: only runs on release:published, requires both build jobs to pass, then downloads all artifacts and uploads via trusted publishing. Co-Authored-By: Claude Sonnet 4.6 --- .github/workflows/python-publish-pypi.yml | 154 ++++++++++++++-------- 1 file changed, 97 insertions(+), 57 deletions(-) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 719598c17..f955e98e8 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -1,130 +1,170 @@ -# 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 + 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}} pip \ + conda create -n anuga_env -c conda-forge \ + python=${{ matrix.python-version }} pip \ meson-python meson ninja cython pybind11 numpy - - 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 - - # 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. + run: conda install -c conda-forge -n anuga_env libpython gcc_win-64 gxx_win-64 - - 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 pip wheel --no-build-isolation --no-deps -w dist-wheel . - - name: Repair wheels + - name: Repair wheel (Linux — make manylinux-compatible) + if: runner.os == 'Linux' shell: bash -el {0} run: | conda activate anuga_env - python -m pip install repairwheel + pip install repairwheel repairwheel -o dist dist-wheel/*.whl - - name: Build source distribution - if: runner.os == 'Linux' && matrix.python-version == '3.13' + - name: Stage wheel (non-Linux) + if: runner.os != 'Linux' shell: bash -el {0} run: | - conda activate anuga_env - pip install build - python -m build --sdist --no-isolation --outdir dist - - - name: Upload distributions + mkdir -p dist + cp dist-wheel/*.whl dist/ + + - name: Upload wheel artifact uses: actions/upload-artifact@v7 with: - name: release-dists-${{ matrix.python-version }}-${{ matrix.os }} - path: dist/ + 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 + - 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 + - name: Build sdist + shell: bash -el {0} + run: | + conda activate anuga_env + pip install build + python -m build --sdist --no-isolation --outdir dist - pypi-publish: - name: Upload release to PyPI + - name: Upload sdist artifact + uses: actions/upload-artifact@v7 + with: + 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) 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/ - From af933f43a29a05fd6c9a10485fe82741ee7942d1 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Mon, 11 May 2026 20:37:16 +1000 Subject: [PATCH 16/28] fix(ci): add pkg-config and PKG_CONFIG_PATH to fix Linux wheel builds meson falls back to /usr/bin/pkg-config (system) when no pkg-config is specified in the native file. The system pkg-config finds Ubuntu's Python 3.12 headers which have a split pyconfig.h that doesn't exist, causing the Cython sanity check to fail. Two-pronged fix: - Add pkg-config to conda create so conda's pkg-config (aware of the conda env) shadows the system /usr/bin/pkg-config on PATH. - Export PKG_CONFIG_PATH=$CONDA_PREFIX/lib/pkgconfig:... before the build so meson finds the conda Python's .pc files first. Co-Authored-By: Claude Sonnet 4.6 --- .github/workflows/python-publish-pypi.yml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index f955e98e8..4f5d25814 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -46,7 +46,7 @@ jobs: run: | conda create -n anuga_env -c conda-forge \ python=${{ matrix.python-version }} pip \ - meson-python meson ninja cython pybind11 numpy + meson-python meson ninja cython pybind11 numpy pkg-config - name: Install compilers (Windows) if: runner.os == 'Windows' @@ -67,6 +67,7 @@ jobs: shell: bash -el {0} run: | conda activate anuga_env + export PKG_CONFIG_PATH="$CONDA_PREFIX/lib/pkgconfig:$CONDA_PREFIX/share/pkgconfig" pip wheel --no-build-isolation --no-deps -w dist-wheel . - name: Repair wheel (Linux — make manylinux-compatible) @@ -113,12 +114,13 @@ jobs: shell: bash -el {0} run: | conda create -n anuga_env -c conda-forge python=3.13 pip \ - meson-python meson ninja cython pybind11 numpy + meson-python meson ninja cython pybind11 numpy pkg-config - name: Build sdist shell: bash -el {0} run: | conda activate anuga_env + export PKG_CONFIG_PATH="$CONDA_PREFIX/lib/pkgconfig:$CONDA_PREFIX/share/pkgconfig" pip install build python -m build --sdist --no-isolation --outdir dist From e6aaf2a2e5d6b2650729a51c5047a8fdefc94149 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Mon, 11 May 2026 20:40:26 +1000 Subject: [PATCH 17/28] fix(ci): drop --no-isolation for sdist build to avoid ninja dep check python -m build --no-isolation checks that all build-system.requires entries are importable Python packages. conda's ninja binary is not a Python package, so the check fails with 'Missing dependencies: ninja'. For the sdist, drop --no-isolation so the build isolation venv installs PyPI's ninja (which does register as a package). PKG_CONFIG_PATH is still exported and is inherited by the isolated venv subprocess, so meson finds the conda env's Python headers rather than the broken system ones. Co-Authored-By: Claude Sonnet 4.6 --- .github/workflows/python-publish-pypi.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 4f5d25814..5ff2e018d 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -120,9 +120,11 @@ jobs: 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 --no-isolation --outdir dist + python -m build --sdist --outdir dist - name: Upload sdist artifact uses: actions/upload-artifact@v7 From b89ae4df5d357bbd19f00dd585f925b3b52730be Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Mon, 11 May 2026 21:05:18 +1000 Subject: [PATCH 18/28] feat(ci): upload wheels and sdist as GitHub Release assets After publishing to PyPI, upload all dist/* files to the GitHub Release using gh release upload. Requires contents:write permission on the publish job. Co-Authored-By: Claude Sonnet 4.6 --- .github/workflows/python-publish-pypi.yml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 5ff2e018d..293ce2d9a 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -148,6 +148,7 @@ jobs: permissions: id-token: write # required for trusted publishing (OIDC) + contents: write # required to upload release assets steps: - name: Download all wheel artifacts @@ -172,3 +173,8 @@ jobs: 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 }} From 6d982aa5d05987ca5ea1c251f9253d5db240dfbb Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 1 Apr 2026 17:06:19 +1100 Subject: [PATCH 19/28] Replace utm dependency with pyproj via epsg_to_ll/ll_to_epsg - Add vectorised epsg_to_ll() and ll_to_epsg() to redfearn.py using pyproj - Rewrite UTMtoLL and LLtoUTM as thin wrappers (public API preserved) - Replace utm.to_latlon() in quantity.py with epsg_to_ll via geo_reference.get_epsg() - Replace scalar UTMtoLL loop in geospatial_data.py with vectorised epsg_to_ll - Replace utm.from_latlon() in test_quantity.py with LLtoUTM - Remove unused `from affine import Affine` in tif2point_values.py - Move pyproj>=3.6 from optional [data] to required dependencies - Remove utm and affine from pyproject.toml Co-Authored-By: Claude Sonnet 4.6 --- anuga/__init__.py | 3 + anuga/abstract_2d_finite_volumes/quantity.py | 32 +-- .../tests/test_quantity.py | 40 ++-- .../lat_long_UTM_conversion.py | 204 +++++------------- anuga/coordinate_transforms/redfearn.py | 46 ++++ anuga/file_conversion/tif2point_values.py | 1 - anuga/geospatial_data/geospatial_data.py | 14 +- pyproject.toml | 4 +- 8 files changed, 131 insertions(+), 213 deletions(-) diff --git a/anuga/__init__.py b/anuga/__init__.py index b07070e35..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, \ @@ -438,6 +439,8 @@ def get_args(): # Coordinate transforms 'LLtoUTM', 'UTMtoLL', + 'epsg_to_ll', + 'll_to_epsg', # Parallel API 'barrier', 'collect_value', 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_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/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..9b261df06 100644 --- a/anuga/coordinate_transforms/redfearn.py +++ b/anuga/coordinate_transforms/redfearn.py @@ -10,6 +10,52 @@ 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) + lon, lat = t.transform(easting, northing) + 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) + easting, northing = t.transform(lon, lat) + 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 421d27afe..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 diff --git a/anuga/geospatial_data/geospatial_data.py b/anuga/geospatial_data/geospatial_data.py index 197983f0e..b1b9b0dba 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,12 @@ 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() + 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/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', From e0e8e9245a33a2d1da9efd5ff2b2f0c2fc020cca Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 13 May 2026 09:25:25 +1000 Subject: [PATCH 20/28] feat: add GitHub Codespaces devcontainer MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds .devcontainer/devcontainer.json and a dedicated environment.yml so users can open ANUGA in a fully configured cloud environment by clicking 'Code → Open with Codespaces' on GitHub — no local install needed. Setup: - onCreateCommand: creates the conda env (cached by Codespaces prebuilds) - postCreateCommand: pip install --no-build-isolation -e . (needs source) - Activates anuga_env in every new terminal via postStartCommand - Default Python interpreter and Jupyter kernel point to anuga_env - Extensions: ms-python.python, ms-python.pylance, ms-toolsai.jupyter The devcontainer environment.yml is kept separate from the main environments/ files so it can evolve independently and stays clean (no utm, correct pyproj as hard dep, jupyter/ipykernel included). Co-Authored-By: Claude Sonnet 4.6 --- .devcontainer/devcontainer.json | 37 +++++++++++++++++++++++++++++++ .devcontainer/environment.yml | 39 +++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+) create mode 100644 .devcontainer/devcontainer.json create mode 100644 .devcontainer/environment.yml diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 000000000..f627bdd17 --- /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. + "postCreateCommand": "conda run -n anuga_env pip install --no-build-isolation -e .", + + // 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 From b81523a9dff6203a2c5713c6efd0c20c23ddb65a Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 13 May 2026 13:36:06 +1000 Subject: [PATCH 21/28] fix(devcontainer): register anuga_env Jupyter kernel on creation Co-Authored-By: Claude Sonnet 4.6 --- .devcontainer/devcontainer.json | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index f627bdd17..f424cdd4b 100644 --- a/.devcontainer/devcontainer.json +++ b/.devcontainer/devcontainer.json @@ -5,8 +5,8 @@ // 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. - "postCreateCommand": "conda run -n anuga_env pip install --no-build-isolation -e .", + // 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", From 8e32c448e9dee288bb2bcbdf9738edabf42b0dcc Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 13 May 2026 13:57:47 +1000 Subject: [PATCH 22/28] docs: add GitHub Codespaces badge and section to README Co-Authored-By: Claude Sonnet 4.6 --- README.rst | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) 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 ---------------------- From a6e22a13e5859aa19fcef2fe09394dc9a6c17630 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 13 May 2026 14:57:28 +1000 Subject: [PATCH 23/28] feat: add get_epsg/set_epsg convenience methods to Domain Co-Authored-By: Claude Sonnet 4.6 --- anuga/abstract_2d_finite_volumes/generic_domain.py | 8 ++++++++ 1 file changed, 8 insertions(+) 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) From 7c24376651943feef39666f5faebb5ea2cacdf68 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 13 May 2026 18:46:50 +1000 Subject: [PATCH 24/28] =?UTF-8?q?fix:=20fallback=20zone=E2=86=92EPSG=20in?= =?UTF-8?q?=20get=5Fdata=5Fpoints=20when=20hemisphere=20is=20undefined?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Geo_reference with only a zone set defaults to hemisphere='undefined', so get_epsg() returns None. Use the isSouthHemisphere parameter to infer the UTM EPSG (32700+zone or 32600+zone) when no explicit EPSG is stored. Co-Authored-By: Claude Sonnet 4.6 --- anuga/geospatial_data/geospatial_data.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/anuga/geospatial_data/geospatial_data.py b/anuga/geospatial_data/geospatial_data.py index b1b9b0dba..c63bd22d5 100644 --- a/anuga/geospatial_data/geospatial_data.py +++ b/anuga/geospatial_data/geospatial_data.py @@ -355,6 +355,9 @@ def get_data_points(self, 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())) From 05fb821b6814f7042a79cc1b85a2a8c7497c71f8 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 13 May 2026 19:05:19 +1000 Subject: [PATCH 25/28] fix: apply NumPy >= 2.0 scalar fix to main's epsg_to_ll/ll_to_epsg main's redfearn.py was missing the e.item()/n.item() workaround that prevent pyproj's scalar path from calling float(ndim>0 array). Co-Authored-By: Claude Sonnet 4.6 --- anuga/coordinate_transforms/redfearn.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/anuga/coordinate_transforms/redfearn.py b/anuga/coordinate_transforms/redfearn.py index 9b261df06..76cf66cf3 100644 --- a/anuga/coordinate_transforms/redfearn.py +++ b/anuga/coordinate_transforms/redfearn.py @@ -30,7 +30,16 @@ def epsg_to_ll(easting, northing, epsg): """ from pyproj import Transformer t = Transformer.from_crs(epsg, 4326, always_xy=True) - lon, lat = t.transform(easting, northing) + 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) @@ -53,7 +62,12 @@ def ll_to_epsg(lat, lon, epsg): """ from pyproj import Transformer t = Transformer.from_crs(4326, epsg, always_xy=True) - easting, northing = t.transform(lon, lat) + 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): From 2eae982c0efe3bbbad88e5a956f0c3cb5b38d68a Mon Sep 17 00:00:00 2001 From: gaoflow Date: Mon, 1 Jun 2026 21:22:50 +0200 Subject: [PATCH 26/28] Fix NameError in util.get_gauges_from_file get_gauges_from_file called gauge_get_from_file, which lives in abstract_2d_finite_volumes/gauge.py and was never imported into util, so every call raised NameError. Import it locally inside the function (gauge imports from util at module level, so a top-level import would be circular) and add a regression test that parses a small gauge file. --- .../tests/test_util.py | 26 +++++++++++++++++-- anuga/abstract_2d_finite_volumes/util.py | 2 ++ 2 files changed, 26 insertions(+), 2 deletions(-) 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) From 1fd1c71eea055d377d2888abb15f02f536c420e4 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Wed, 3 Jun 2026 18:55:09 +1000 Subject: [PATCH 27/28] fix: strip whitespace from gauge names in gauge_get_from_file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit loc.strip(r'\n') only stripped newlines, leaving leading spaces from comma-separated fields (e.g. "308500, 6193000, gauge_a" → " gauge_a"). Co-Authored-By: Claude Sonnet 4.6 --- anuga/abstract_2d_finite_volumes/gauge.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From e4700a7b8e2aba666488d572ab58ba8f5bbfa416 Mon Sep 17 00:00:00 2001 From: Stephen Roberts Date: Sat, 6 Jun 2026 20:17:10 +1000 Subject: [PATCH 28/28] Repair wheels on all platforms so they are self-contained MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Previously only Linux wheels were repaired (auditwheel); macOS and Windows wheels were copied to dist/ unrepaired. The macOS build links conda's libomp via an @rpath into the build env, so `import anuga` from the published arm64 wheel fails on a plain pip install with "Library not loaded: @rpath/libomp.dylib". Run repairwheel for every platform — it dispatches to auditwheel (Linux), delocate (macOS, bundles libomp.dylib and rewrites the rpath) and delvewheel (Windows, bundles dependent DLLs) — so all published wheels are self-contained. Conda lib dirs that exist on the OS are passed via -l so dependent libs are found at repair time. Co-Authored-By: Claude Opus 4.8 --- .github/workflows/python-publish-pypi.yml | 24 +++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/.github/workflows/python-publish-pypi.yml b/.github/workflows/python-publish-pypi.yml index 293ce2d9a..045b64013 100644 --- a/.github/workflows/python-publish-pypi.yml +++ b/.github/workflows/python-publish-pypi.yml @@ -70,20 +70,24 @@ jobs: export PKG_CONFIG_PATH="$CONDA_PREFIX/lib/pkgconfig:$CONDA_PREFIX/share/pkgconfig" pip wheel --no-build-isolation --no-deps -w dist-wheel . - - name: Repair wheel (Linux — make manylinux-compatible) - if: runner.os == 'Linux' + - name: Repair wheel (bundle external libs; self-contained per platform) shell: bash -el {0} run: | conda activate anuga_env pip install repairwheel - repairwheel -o dist dist-wheel/*.whl - - - name: Stage wheel (non-Linux) - if: runner.os != 'Linux' - shell: bash -el {0} - run: | - mkdir -p dist - cp dist-wheel/*.whl dist/ + # 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