From 74b736d0b630637805217a2b678b6cc508647099 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 7 Oct 2025 09:56:07 +0100 Subject: [PATCH 01/39] tests: Ruff fixes for FCI runtest --- tests/MMS/spatial/fci/runtest | 238 +++++++++++++++------------------- 1 file changed, 105 insertions(+), 133 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 7a9d6e655e..1f0b297ebc 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -6,19 +6,28 @@ # Cores: 2 # requires: zoidberg -from boututils.run_wrapper import build_and_log, launch_safe -from boutdata.collect import collect +import pathlib +import pickle +import sys + import boutconfig as conf +import zoidberg as zb +from boutdata.collect import collect +from boututils.run_wrapper import build_and_log, launch_safe +from numpy import arange, array, linspace, log, polyfit +from scipy.interpolate import RectBivariateSpline as RBS -from numpy import array, log, polyfit, linspace, arange -import pickle +def myRBS(a, b, c): + mx, _ = c.shape + kx = max(mx - 1, 1) + kx = min(kx, 3) + return RBS(a, b, c, kx=kx) -from sys import stdout -import zoidberg as zb +zb.poloidal_grid.RectBivariateSpline = myRBS -nx = 4 # Not changed for these tests +nx = 3 # Not changed for these tests # Resolution in y and z nlist = [8, 16, 32, 64, 128] @@ -31,7 +40,6 @@ directory = "data" nproc = 2 mthread = 2 - success = True error_2 = {} @@ -46,115 +54,88 @@ failures = [] build_and_log("FCI MMS test") for nslice in nslices: - for method in [ - "hermitespline", - "lagrange4pt", - "bilinear", - # "monotonichermitespline", - ]: - error_2[nslice] = [] - error_inf[nslice] = [] - - # Which central difference scheme to use and its expected order - order = nslice * 2 - method_orders[nslice] = {"name": "C{}".format(order), "order": order} - - for n in nlist: - # Define the magnetic field using new poloidal gridding method - # Note that the Bz and Bzprime parameters here must be the same as in mms.py - field = zb.field.Slab(Bz=0.05, Bzprime=0.1) - # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid( - nx, n, 0.1, 1.0, MXG=1 - ) - # Set the ylength and y locations - ylength = 10.0 - - if yperiodic: - ycoords = linspace(0.0, ylength, n, endpoint=False) - else: - # Doesn't include the end points - ycoords = (arange(n) + 0.5) * ylength / float(n) - - # Create the grid - grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) - # Make and write maps - maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) - zb.write_maps( - grid, - field, - maps, - new_names=False, - metric2d=conf.isMetric2D(), - quiet=True, - ) - - args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={} NXPE={}".format( - n, - nslice, - yperiodic, - method_orders[nslice]["name"], - 2 if conf.has["petsc"] and method == "hermitespline" else 1, - ) - args += f" mesh:paralleltransform:xzinterpolation:type={method}" - - # Command to run - cmd = "./fci_mms " + args - - print("Running command: " + cmd) - - # Launch using MPI - s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) - - # Save output to log file - with open("run.log." + str(n), "w") as f: - f.write(out) - - if s: - print("Run failed!\nOutput was:\n") - print(out) - exit(s) - - # Collect data - l_2 = collect( - "l_2", - tind=[1, 1], - info=False, - path=directory, - xguards=False, - yguards=False, - ) - l_inf = collect( - "l_inf", - tind=[1, 1], - info=False, - path=directory, - xguards=False, - yguards=False, - ) - - error_2[nslice].append(l_2) - error_inf[nslice].append(l_inf) - - print("Errors : l-2 {:f} l-inf {:f}".format(l_2, l_inf)) - - dx = 1.0 / array(nlist) - - # Calculate convergence order - fit = polyfit(log(dx), log(error_2[nslice]), 1) - order = fit[0] - stdout.write("Convergence order = {:f} (fit)".format(order)) - - order = log(error_2[nslice][-2] / error_2[nslice][-1]) / log(dx[-2] / dx[-1]) - stdout.write(", {:f} (small spacing)".format(order)) - - # Should be close to the expected order - if order > method_orders[nslice]["order"] * 0.95: - print("............ PASS\n") + error_2[nslice] = [] + error_inf[nslice] = [] + + # Which central difference scheme to use and its expected order + order = nslice * 2 + name = f"C{order}" + method_orders[nslice] = {"name": name, "order": order} + + for n in nlist: + # Define the magnetic field using new poloidal gridding method + # Note that the Bz and Bzprime parameters here must be the same as in mms.py + field = zb.field.Slab(Bz=0.05, Bzprime=0.1) + # Create rectangular poloidal grids + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.0, MXG=1) + # Set the ylength and y locations + ylength = 10.0 + + if yperiodic: + ycoords = linspace(0.0, ylength, n, endpoint=False) else: - print("............ FAIL\n") - success = False - failures.append(method_orders[nslice]["name"]) + # Doesn't include the end points + ycoords = (arange(n) + 0.5) * ylength / float(n) + + # Create the grid + grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) + # Make and write maps + maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) + zb.write_maps( + grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True + ) + + # Command to run + args = f" MZ={n} MYG={nslice} mesh:paralleltransform:y_periodic={yperiodic} mesh:ddy:first={name}" + cmd = f"./fci_mms {args}" + + print(f"Running command: {cmd}") + + # Launch using MPI + s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) + + # Save output to log file + pathlib.Path(f"run.log.{n}").write_text(out) + + if s: + print(f"Run failed!\nOutput was:\n{out}") + sys.exit(s) + + # Collect data + l_2 = collect( + "l_2", tind=[1, 1], info=False, path=directory, xguards=False, yguards=False + ) + l_inf = collect( + "l_inf", + tind=[1, 1], + info=False, + path=directory, + xguards=False, + yguards=False, + ) + + error_2[nslice].append(l_2) + error_inf[nslice].append(l_inf) + + print(f"Errors : l-2 {l_2:f} l-inf {l_inf:f}") + + dx = 1.0 / array(nlist) + + # Calculate convergence order + fit = polyfit(log(dx), log(error_2[nslice]), 1) + order = fit[0] + print(f"Convergence order = {order:f} (fit)", end="") + + order = log(error_2[nslice][-2] / error_2[nslice][-1]) / log(dx[-2] / dx[-1]) + print(f", {order:f} (small spacing)") + + # Should be close to the expected order + if order > order * 0.95: + print("............ PASS\n") + else: + print("............ FAIL\n") + success = False + failures.append(name) with open("fci_mms.pkl", "wb") as output: @@ -164,7 +145,7 @@ with open("fci_mms.pkl", "wb") as output: pickle.dump(error_inf[nslice], output) # Do we want to show the plot as well as save it to file. -showPlot = True +show_plot = True if False: try: @@ -174,18 +155,9 @@ if False: fig, ax = plt.subplots(1, 1) for nslice in nslices: - ax.plot( - dx, - error_2[nslice], - "-", - label="{} $l_2$".format(method_orders[nslice]["name"]), - ) - ax.plot( - dx, - error_inf[nslice], - "--", - label="{} $l_\\inf$".format(method_orders[nslice]["name"]), - ) + name = method_orders[nslice]["name"] + ax.plot(dx, error_2[nslice], "-", label=f"{name} $l_2$") + ax.plot(dx, error_inf[nslice], "--", label=f"{name} $l_\\inf$") ax.legend(loc="upper left") ax.grid() ax.set_yscale("log") @@ -198,7 +170,7 @@ if False: print("Plot saved to fci_mms.pdf") - if showPlot: + if show_plot: plt.show() plt.close() except ImportError: @@ -206,9 +178,9 @@ if False: if success: print("All tests passed") - exit(0) + sys.exit(0) else: print("Some tests failed:") for failure in failures: - print("\t" + failure) - exit(1) + print(f"\t{failure}") + sys.exit(1) From eb589e16948c795bc4137e85e72d77cd452c3142 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 7 Oct 2025 12:05:01 +0100 Subject: [PATCH 02/39] tests: Expand FCI MMS test to `Grad2_par2` --- CMakeLists.txt | 1 + tests/MMS/spatial/fci/data/BOUT.inp | 5 +- tests/MMS/spatial/fci/fci_mms.cxx | 48 ++++++++++----- tests/MMS/spatial/fci/mms.py | 5 +- tests/MMS/spatial/fci/runtest | 96 +++++++++++++++-------------- 5 files changed, 90 insertions(+), 65 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 44d100e444..eabfc055ee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -376,6 +376,7 @@ if (zoidberg_FOUND EQUAL 0) else() set(zoidberg_FOUND OFF) endif() +message(STATUS "Found Zoidberg for FCI tests: ${zoidberg_FOUND}") option(BOUT_GENERATE_FIELDOPS "Automatically re-generate the Field arithmetic operators from the Python templates. \ Requires Python3, clang-format, and Jinja2. Turn this OFF to skip generating them if, for example, \ diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 5f2001a906..f065989524 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -1,7 +1,6 @@ - input_field = sin(y - 2*z) + sin(y - z) - -solution = (6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) +grad_par_solution = (6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) +grad2_par2_solution = (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))/sqrt((0.01*x + 0.045)^2 + 1.0) MXG = 1 MYG = 1 diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 18405a7f88..48b18f04ef 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -1,5 +1,4 @@ #include "bout/bout.hxx" -#include "bout/derivs.hxx" #include "bout/field_factory.hxx" int main(int argc, char** argv) { @@ -8,29 +7,50 @@ int main(int argc, char** argv) { using bout::globals::mesh; Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; - Field3D solution{FieldFactory::get()->create3D("solution", Options::getRoot(), mesh)}; // Communicate to calculate parallel transform mesh->communicate(input); - Field3D result{Grad_par(input)}; - Field3D error{result - solution}; - Options dump; // Add mesh geometry variables mesh->outputVars(dump); - dump["l_2"] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); - dump["l_inf"] = max(abs(error), true, "RGN_NOBNDRY"); + auto* factory = FieldFactory::get(); + { + Field3D solution{factory->create3D("grad_par_solution", Options::getRoot(), mesh)}; + Field3D result{Grad_par(input)}; + Field3D error{result - solution}; + + dump["grad_par_l_2"] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump["grad_par_l_inf"] = max(abs(error), true, "RGN_NOBNDRY"); - dump["result"] = result; - dump["error"] = error; - dump["input"] = input; - dump["solution"] = solution; + dump["grad_par_result"] = result; + dump["grad_par_error"] = error; + dump["grad_par_input"] = input; + dump["grad_par_solution"] = solution; - for (int slice = 1; slice < mesh->ystart; ++slice) { - dump[fmt::format("input.ynext(-{})", slice)] = input.ynext(-slice); - dump[fmt::format("input.ynext({})", slice)] = input.ynext(slice); + for (int slice = 1; slice < mesh->ystart; ++slice) { + dump[fmt::format("grad_par_input.ynext(-{})", slice)] = input.ynext(-slice); + dump[fmt::format("grad_par_input.ynext({})", slice)] = input.ynext(slice); + } + } + { + Field3D solution{factory->create3D("grad2_par2_solution", Options::getRoot(), mesh)}; + Field3D result{Grad2_par2(input)}; + Field3D error{result - solution}; + + dump["grad2_par2_l_2"] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump["grad2_par2_l_inf"] = max(abs(error), true, "RGN_NOBNDRY"); + + dump["grad2_par2_result"] = result; + dump["grad2_par2_error"] = error; + dump["grad2_par2_input"] = input; + dump["grad2_par2_solution"] = solution; + + for (int slice = 1; slice < mesh->ystart; ++slice) { + dump[fmt::format("grad2_par2_input.ynext(-{})", slice)] = input.ynext(-slice); + dump[fmt::format("grad2_par2_input.ynext({})", slice)] = input.ynext(slice); + } } bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 1e71135c90..ae48ef8f06 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -30,5 +30,6 @@ def FCI_ddy(f): ############################################ # Equations solved -print("input = " + exprToStr(f)) -print("solution = " + exprToStr(FCI_ddy(f))) +print(f"input_field = {exprToStr(f)}") +print(f"grad_par_solution = {exprToStr(FCI_ddy(f))}") +print(f"grad2_par2_solution = {exprToStr(FCI_ddy(FCI_ddy(f)))}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 1f0b297ebc..27b18baafd 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -9,6 +9,7 @@ import pathlib import pickle import sys +from collections import defaultdict import boutconfig as conf import zoidberg as zb @@ -19,6 +20,7 @@ from scipy.interpolate import RectBivariateSpline as RBS def myRBS(a, b, c): + """RectBivariateSpline, but automatically tune spline degree for small arrays""" mx, _ = c.shape kx = max(mx - 1, 1) kx = min(kx, 3) @@ -27,6 +29,13 @@ def myRBS(a, b, c): zb.poloidal_grid.RectBivariateSpline = myRBS + +def quiet_collect(name: str): + return collect( + name, tind=[1, 1], info=False, path=directory, xguards=False, yguards=False, + ) + + nx = 3 # Not changed for these tests # Resolution in y and z @@ -50,12 +59,29 @@ method_orders = {} yperiodic = True failures = [] +operators = ("grad_par", "grad2_par2") build_and_log("FCI MMS test") + +def assert_convergence(error, dx, name, order) -> bool: + fit = polyfit(log(dx), log(error), 1) + order = fit[0] + print(f"{name} convergence order = {order:f} (fit)", end="") + + order = log(error[-2] / error[-1]) / log(dx[-2] / dx[-1]) + print(f", {order:f} (small spacing)", end="") + + # Should be close to the expected order + success = order > order * 0.95 + print(f" ............ {'PASS' if success else 'FAIL'}") + + return success + + for nslice in nslices: - error_2[nslice] = [] - error_inf[nslice] = [] + error_2[nslice] = defaultdict(list) + error_inf[nslice] = defaultdict(list) # Which central difference scheme to use and its expected order order = nslice * 2 @@ -79,70 +105,48 @@ for nslice in nslices: # Create the grid grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) - # Make and write maps maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) zb.write_maps( - grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True + grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True, ) # Command to run args = f" MZ={n} MYG={nslice} mesh:paralleltransform:y_periodic={yperiodic} mesh:ddy:first={name}" cmd = f"./fci_mms {args}" - print(f"Running command: {cmd}") # Launch using MPI - s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) + status, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) # Save output to log file pathlib.Path(f"run.log.{n}").write_text(out) - if s: + if status: print(f"Run failed!\nOutput was:\n{out}") - sys.exit(s) + sys.exit(status) # Collect data - l_2 = collect( - "l_2", tind=[1, 1], info=False, path=directory, xguards=False, yguards=False - ) - l_inf = collect( - "l_inf", - tind=[1, 1], - info=False, - path=directory, - xguards=False, - yguards=False, - ) + for operator in operators: + l_2 = quiet_collect(f"{operator}_l_2") + l_inf = quiet_collect(f"{operator}_l_inf") - error_2[nslice].append(l_2) - error_inf[nslice].append(l_inf) + error_2[nslice][operator].append(l_2) + error_inf[nslice][operator].append(l_inf) - print(f"Errors : l-2 {l_2:f} l-inf {l_inf:f}") + print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") dx = 1.0 / array(nlist) + for operator in operators: + test_name = f"{operator} {name}" + success &= assert_convergence(error_2[nslice][operator], dx, test_name, order) + if not success: + failures.append(test_name) - # Calculate convergence order - fit = polyfit(log(dx), log(error_2[nslice]), 1) - order = fit[0] - print(f"Convergence order = {order:f} (fit)", end="") - - order = log(error_2[nslice][-2] / error_2[nslice][-1]) / log(dx[-2] / dx[-1]) - print(f", {order:f} (small spacing)") - - # Should be close to the expected order - if order > order * 0.95: - print("............ PASS\n") - else: - print("............ FAIL\n") - success = False - failures.append(name) - -with open("fci_mms.pkl", "wb") as output: +with pathlib.Path("fci_mms.pkl").open("wb") as output: pickle.dump(nlist, output) - for nslice in nslices: - pickle.dump(error_2[nslice], output) - pickle.dump(error_inf[nslice], output) + pickle.dump(error_2, output) + pickle.dump(error_inf, output) # Do we want to show the plot as well as save it to file. show_plot = True @@ -177,10 +181,10 @@ if False: print("No matplotlib") if success: - print("All tests passed") - sys.exit(0) + print("\nAll tests passed") else: - print("Some tests failed:") + print("\nSome tests failed:") for failure in failures: print(f"\t{failure}") - sys.exit(1) + +sys.exit(0 if success else 1) From c0cc84b78ffe127fea61414db92d92c23743906b Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 7 Oct 2025 15:10:32 +0100 Subject: [PATCH 03/39] tests: Generalise FCI MMS test to allow for more cases --- tests/MMS/spatial/fci/runtest | 302 +++++++++++++++++++++------------- 1 file changed, 188 insertions(+), 114 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 27b18baafd..c575f1afbc 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -6,10 +6,12 @@ # Cores: 2 # requires: zoidberg +import argparse +import json import pathlib -import pickle import sys -from collections import defaultdict +from time import time +from typing import Any import boutconfig as conf import zoidberg as zb @@ -18,6 +20,16 @@ from boututils.run_wrapper import build_and_log, launch_safe from numpy import arange, array, linspace, log, polyfit from scipy.interpolate import RectBivariateSpline as RBS +# Global parameters +DIRECTORY = "data" +NPROC = 2 +MTHREAD = 2 +OPERATORS = ("grad_par", "grad2_par2") +NX = 3 +# Resolution in y and z +NLIST = [8, 16, 32, 64, 128] +dx = 1.0 / array(NLIST) + def myRBS(a, b, c): """RectBivariateSpline, but automatically tune spline degree for small arrays""" @@ -30,38 +42,16 @@ def myRBS(a, b, c): zb.poloidal_grid.RectBivariateSpline = myRBS -def quiet_collect(name: str): +def quiet_collect(name: str) -> float: + # Index to return a plain (numpy) float rather than `BoutArray` return collect( - name, tind=[1, 1], info=False, path=directory, xguards=False, yguards=False, - ) - - -nx = 3 # Not changed for these tests - -# Resolution in y and z -nlist = [8, 16, 32, 64, 128] - -# Number of parallel slices (in each direction) -nslices = [1] - -directory = "data" - -nproc = 2 -mthread = 2 - -success = True - -error_2 = {} -error_inf = {} -method_orders = {} - -# Run with periodic Y? -yperiodic = True - -failures = [] -operators = ("grad_par", "grad2_par2") - -build_and_log("FCI MMS test") + name, + tind=[1, 1], + info=False, + path=DIRECTORY, + xguards=False, + yguards=False, + )[()] def assert_convergence(error, dx, name, order) -> bool: @@ -74,117 +64,201 @@ def assert_convergence(error, dx, name, order) -> bool: # Should be close to the expected order success = order > order * 0.95 - print(f" ............ {'PASS' if success else 'FAIL'}") + print(f"\t............ {'PASS' if success else 'FAIL'}") return success -for nslice in nslices: - error_2[nslice] = defaultdict(list) - error_inf[nslice] = defaultdict(list) +def run_fci_operators( + nslice: int, nz: int, yperiodic: bool, name: str +) -> dict[str, float]: + # Define the magnetic field using new poloidal gridding method + # Note that the Bz and Bzprime parameters here must be the same as in mms.py + field = zb.field.Slab(Bz=0.05, Bzprime=0.1) + # Create rectangular poloidal grids + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(NX, nz, 0.1, 1.0, MXG=1) + # Set the ylength and y locations + ylength = 10.0 + + if yperiodic: + ycoords = linspace(0.0, ylength, nz, endpoint=False) + else: + # Doesn't include the end points + ycoords = (arange(nz) + 0.5) * ylength / float(nz) + + # Create the grid + grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) + maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) + zb.write_maps( + grid, + field, + maps, + new_names=False, + metric2d=conf.isMetric2D(), + quiet=True, + ) - # Which central difference scheme to use and its expected order - order = nslice * 2 - name = f"C{order}" - method_orders[nslice] = {"name": name, "order": order} - - for n in nlist: - # Define the magnetic field using new poloidal gridding method - # Note that the Bz and Bzprime parameters here must be the same as in mms.py - field = zb.field.Slab(Bz=0.05, Bzprime=0.1) - # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.0, MXG=1) - # Set the ylength and y locations - ylength = 10.0 - - if yperiodic: - ycoords = linspace(0.0, ylength, n, endpoint=False) - else: - # Doesn't include the end points - ycoords = (arange(n) + 0.5) * ylength / float(n) - - # Create the grid - grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) - maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) - zb.write_maps( - grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True, - ) + # Command to run + args = f"MZ={nz} MYG={nslice} mesh:paralleltransform:y_periodic={yperiodic} mesh:ddy:first={name}" + cmd = f"./fci_mms {args}" + print(f"Running command: {cmd}", end="") + + # Launch using MPI + start = time() + status, out = launch_safe(cmd, nproc=NPROC, mthread=MTHREAD, pipe=True) + print(f" ... done in {time() - start:.3}s") + + # Save output to log file + pathlib.Path(f"run.log.{nz}").write_text(out) + + if status: + print(f"Run failed!\nOutput was:\n{out}") + sys.exit(status) + + return { + operator: { + "l_2": quiet_collect(f"{operator}_l_2"), + "l_inf": quiet_collect(f"{operator}_l_inf"), + } + for operator in OPERATORS + } - # Command to run - args = f" MZ={n} MYG={nslice} mesh:paralleltransform:y_periodic={yperiodic} mesh:ddy:first={name}" - cmd = f"./fci_mms {args}" - print(f"Running command: {cmd}") - # Launch using MPI - status, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) +def transpose( + errors: list[dict[str, dict[str, float]]], +) -> dict[str, dict[str, list[float]]]: + """Turn a list of {operator: error} into a dict of {operator: [errors]}""" - # Save output to log file - pathlib.Path(f"run.log.{n}").write_text(out) + kinds = ("l_2", "l_inf") + result = {operator: {kind: [] for kind in kinds} for operator in OPERATORS} + for error in errors: + for k, v in error.items(): + for kind in kinds: + result[k][kind].append(v[kind]) + return result - if status: - print(f"Run failed!\nOutput was:\n{out}") - sys.exit(status) - # Collect data - for operator in operators: - l_2 = quiet_collect(f"{operator}_l_2") - l_inf = quiet_collect(f"{operator}_l_inf") +def check_fci_operators(case: dict) -> bool: + failures = [] - error_2[nslice][operator].append(l_2) - error_inf[nslice][operator].append(l_inf) + nslice = case["nslice"] + yperiodic = case["yperiodic"] + name = case["name"] + order = case["order"] + + all_errors = [] + + for n in NLIST: + errors = run_fci_operators(nslice, n, yperiodic, name) + all_errors.append(errors) + + for operator in OPERATORS: + l_2 = errors[operator]["l_2"] + l_inf = errors[operator]["l_inf"] print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") - dx = 1.0 / array(nlist) - for operator in operators: + final_errors = transpose(all_errors) + for operator in OPERATORS: test_name = f"{operator} {name}" - success &= assert_convergence(error_2[nslice][operator], dx, test_name, order) + success = assert_convergence( + final_errors[operator]["l_2"], dx, test_name, order + ) if not success: failures.append(test_name) + return final_errors, failures -with pathlib.Path("fci_mms.pkl").open("wb") as output: - pickle.dump(nlist, output) - pickle.dump(error_2, output) - pickle.dump(error_inf, output) - -# Do we want to show the plot as well as save it to file. -show_plot = True -if False: +def make_plots(cases): try: - # Plot using matplotlib if available import matplotlib.pyplot as plt + except ImportError: + print("No matplotlib") + return - fig, ax = plt.subplots(1, 1) + fig, axes = plt.subplots(1, len(OPERATORS), figsize=(9, 4)) - for nslice in nslices: - name = method_orders[nslice]["name"] - ax.plot(dx, error_2[nslice], "-", label=f"{name} $l_2$") - ax.plot(dx, error_inf[nslice], "--", label=f"{name} $l_\\inf$") + for ax, operator in zip(axes, OPERATORS): + for name, case in cases.items(): + ax.loglog(dx, case[operator]["l_2"], "-", label=f"{name} $l_2$") + ax.loglog(dx, case[operator]["l_inf"], "--", label=f"{name} $l_\\inf$") ax.legend(loc="upper left") ax.grid() - ax.set_yscale("log") - ax.set_xscale("log") - ax.set_title("error scaling") + ax.set_title(f"Error scaling for {operator}") ax.set_xlabel(r"Mesh spacing $\delta x$") ax.set_ylabel("Error norm") - plt.savefig("fci_mms.pdf") + fig.tight_layout() + fig.savefig("fci_mms.pdf") + print("Plot saved to fci_mms.pdf") - print("Plot saved to fci_mms.pdf") + if args.show_plots: + plt.show() + plt.close() - if show_plot: - plt.show() - plt.close() - except ImportError: - print("No matplotlib") -if success: - print("\nAll tests passed") -else: - print("\nSome tests failed:") - for failure in failures: - print(f"\t{failure}") +def make_case(nslice: int, yperiodic: bool) -> dict[str, Any]: + """ + nslice: + Number of parallel slices (in each direction) + yperiodic: + Run with periodic Y + """ + order = nslice * 2 + return { + "nslice": nslice, + # Which central difference scheme to use and its expected order + "order": order, + "name": f"C{order}", + "yperiodic": yperiodic, + } + + +if __name__ == "__main__": + build_and_log("FCI MMS test") + + parser = argparse.ArgumentParser("Error scaling test for FCI operators") + parser.add_argument( + "--make-plots", action="store_true", help="Create plots of error scaling" + ) + parser.add_argument( + "--show-plots", + action="store_true", + help="Stop and show plots, implies --make-plots", + ) + parser.add_argument( + "--dump-errors", + type=str, + help="Output file to dump errors as JSON", + default="fci_operator_errors.json", + ) + + args = parser.parse_args() + + success = True + failures = [] + cases = { + "nslice=1": make_case(nslice=1, yperiodic=True), + } + + for case in cases.values(): + error2, failures_ = check_fci_operators(case) + case.update(error2) + failures.extend(failures_) + success &= len(failures) == 0 + + if args.dump_errors: + pathlib.Path(args.dump_errors).write_text(json.dumps(cases)) + + if args.make_plots or args.show_plots: + make_plots(cases) + + if success: + print("\nAll tests passed") + else: + print("\nSome tests failed:") + for failure in failures: + print(f"\t{failure}") -sys.exit(0 if success else 1) + sys.exit(0 if success else 1) From 2eda48cf04bb41abaa8fce6d61336fc5b5dddf2d Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 7 Oct 2025 15:30:58 +0100 Subject: [PATCH 04/39] tests: Add test for FCI `Div_par` --- tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 67 +++++++++++++---------------- tests/MMS/spatial/fci/mms.py | 4 ++ tests/MMS/spatial/fci/runtest | 2 +- 4 files changed, 36 insertions(+), 38 deletions(-) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index f065989524..3103b37e2d 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -1,6 +1,7 @@ input_field = sin(y - 2*z) + sin(y - z) grad_par_solution = (6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) grad2_par2_solution = (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))/sqrt((0.01*x + 0.045)^2 + 1.0) +div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530717959*cos(y - z) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) MXG = 1 MYG = 1 diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 48b18f04ef..bc995424c6 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -1,6 +1,33 @@ #include "bout/bout.hxx" +#include "bout/difops.hxx" +#include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include + +namespace { +auto fci_op_test(const std::string& name, Options& dump, const Field3D& input, + const Field3D& result) { + auto* mesh = input.getMesh(); + Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), + Options::getRoot(), mesh)}; + Field3D error{result - solution}; + + dump[fmt::format("{}_l_2", name)] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump[fmt::format("{}_l_inf", name)] = max(abs(error), true, "RGN_NOBNDRY"); + + dump[fmt::format("{}_result", name)] = result; + dump[fmt::format("{}_error", name)] = error; + dump[fmt::format("{}_input", name)] = input; + dump[fmt::format("{}_solution", name)] = solution; + + for (int slice = 1; slice < mesh->ystart; ++slice) { + dump[fmt::format("{}_input.ynext(-{})", name, slice)] = input.ynext(-slice); + dump[fmt::format("{}_input.ynext({})", name, slice)] = input.ynext(slice); + } +} +} // namespace + int main(int argc, char** argv) { BoutInitialise(argc, argv); @@ -15,43 +42,9 @@ int main(int argc, char** argv) { // Add mesh geometry variables mesh->outputVars(dump); - auto* factory = FieldFactory::get(); - { - Field3D solution{factory->create3D("grad_par_solution", Options::getRoot(), mesh)}; - Field3D result{Grad_par(input)}; - Field3D error{result - solution}; - - dump["grad_par_l_2"] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); - dump["grad_par_l_inf"] = max(abs(error), true, "RGN_NOBNDRY"); - - dump["grad_par_result"] = result; - dump["grad_par_error"] = error; - dump["grad_par_input"] = input; - dump["grad_par_solution"] = solution; - - for (int slice = 1; slice < mesh->ystart; ++slice) { - dump[fmt::format("grad_par_input.ynext(-{})", slice)] = input.ynext(-slice); - dump[fmt::format("grad_par_input.ynext({})", slice)] = input.ynext(slice); - } - } - { - Field3D solution{factory->create3D("grad2_par2_solution", Options::getRoot(), mesh)}; - Field3D result{Grad2_par2(input)}; - Field3D error{result - solution}; - - dump["grad2_par2_l_2"] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); - dump["grad2_par2_l_inf"] = max(abs(error), true, "RGN_NOBNDRY"); - - dump["grad2_par2_result"] = result; - dump["grad2_par2_error"] = error; - dump["grad2_par2_input"] = input; - dump["grad2_par2_solution"] = solution; - - for (int slice = 1; slice < mesh->ystart; ++slice) { - dump[fmt::format("grad2_par2_input.ynext(-{})", slice)] = input.ynext(-slice); - dump[fmt::format("grad2_par2_input.ynext({})", slice)] = input.ynext(slice); - } - } + fci_op_test("grad_par", dump, input, Grad_par(input)); + fci_op_test("grad2_par2", dump, input, Grad2_par2(input)); + fci_op_test("div_par", dump, input, Div_par(input)); bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index ae48ef8f06..5f5d48bee2 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -26,6 +26,9 @@ def FCI_ddy(f): return (Bt * diff(f, y) * 2.0 * pi / Ly + Bpx * diff(f, z) * 2.0 * pi / Lz) / B +def FCI_div_par(f): + return Bpx * FCI_ddy(f / Bpx) + ############################################ # Equations solved @@ -33,3 +36,4 @@ def FCI_ddy(f): print(f"input_field = {exprToStr(f)}") print(f"grad_par_solution = {exprToStr(FCI_ddy(f))}") print(f"grad2_par2_solution = {exprToStr(FCI_ddy(FCI_ddy(f)))}") +print(f"div_par_solution = {exprToStr(FCI_div_par(f))}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index c575f1afbc..c98d924e7b 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -24,7 +24,7 @@ from scipy.interpolate import RectBivariateSpline as RBS DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ("grad_par", "grad2_par2") +OPERATORS = ("grad_par", "grad2_par2", "div_par") NX = 3 # Resolution in y and z NLIST = [8, 16, 32, 64, 128] From 636bf288f616f67e2d84552e569a24d64d8f5485 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 7 Oct 2025 15:47:01 +0100 Subject: [PATCH 05/39] tests: Make MMS script update input in-place --- tests/MMS/spatial/fci/data/BOUT.inp | 7 +++--- tests/MMS/spatial/fci/mms.py | 36 +++++++++++++++++++++-------- 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 3103b37e2d..f44773667c 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -3,13 +3,12 @@ grad_par_solution = (6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y grad2_par2_solution = (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))/sqrt((0.01*x + 0.045)^2 + 1.0) div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530717959*cos(y - z) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) -MXG = 1 -MYG = 1 -NXPE = 1 - [mesh] symmetricglobalx = true file = fci.grid.nc +MXG = 1 +MYG = 1 +NXPE = 1 [mesh:ddy] first = C2 diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 5f5d48bee2..e434ea3ecf 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -3,11 +3,15 @@ # Generate manufactured solution and sources for FCI test # -from boutdata.mms import * +from math import pi +import warnings -from sympy import sin, cos, sqrt +from boututils.boutwarnings import AlwaysWarning +from boutdata.data import BoutOptionsFile +from boutdata.mms import diff, exprToStr, x, y, z +from sympy import sin, sqrt, Expr -from math import pi +warnings.simplefilter("ignore", AlwaysWarning) f = sin(y - z) + sin(y - 2 * z) @@ -23,17 +27,31 @@ B = sqrt(Bpx**2 + Bt**2) -def FCI_ddy(f): +def FCI_ddy(f: Expr) -> Expr: return (Bt * diff(f, y) * 2.0 * pi / Ly + Bpx * diff(f, z) * 2.0 * pi / Lz) / B -def FCI_div_par(f): + +def FCI_div_par(f: Expr) -> Expr: return Bpx * FCI_ddy(f / Bpx) ############################################ # Equations solved -print(f"input_field = {exprToStr(f)}") -print(f"grad_par_solution = {exprToStr(FCI_ddy(f))}") -print(f"grad2_par2_solution = {exprToStr(FCI_ddy(FCI_ddy(f)))}") -print(f"div_par_solution = {exprToStr(FCI_div_par(f))}") +input_field = exprToStr(f) +grad_par_solution = exprToStr(FCI_ddy(f)) +grad2_par2_solution = exprToStr(FCI_ddy(FCI_ddy(f))) +div_par_solution = exprToStr(FCI_div_par(f)) + +print(f"input_field = {input_field}") +print(f"grad_par_solution = {grad_par_solution}") +print(f"grad2_par2_solution = {grad2_par2_solution}") +print(f"div_par_solution = {div_par_solution}") +print(f"div_par_K_grad_par_solution = {div_par_K_grad_par_solution}") + +options = BoutOptionsFile("data/BOUT.inp") +options["input_field"] = input_field +options["grad_par_solution"] = grad_par_solution +options["grad2_par2_solution"] = grad2_par2_solution +options["div_par_solution"] = div_par_solution +options.write("data/BOUT.inp", overwrite=True) From 16318eccf35685ca4caf7fe9263a06c4d884abe4 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 7 Oct 2025 16:34:44 +0100 Subject: [PATCH 06/39] tests: Add test for FCI `Div_par_K_Grad_par` --- tests/MMS/spatial/fci/data/BOUT.inp | 2 ++ tests/MMS/spatial/fci/fci_mms.cxx | 5 ++++- tests/MMS/spatial/fci/mms.py | 23 ++++++++++++++++++----- tests/MMS/spatial/fci/runtest | 5 +++-- 4 files changed, 27 insertions(+), 8 deletions(-) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index f44773667c..e31d84cb35 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -2,6 +2,8 @@ input_field = sin(y - 2*z) + sin(y - z) grad_par_solution = (6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) grad2_par2_solution = (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))/sqrt((0.01*x + 0.045)^2 + 1.0) div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530717959*cos(y - z) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) +K = cos(y - z) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index bc995424c6..a0b70e9da6 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -2,6 +2,7 @@ #include "bout/difops.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/options.hxx" #include @@ -34,9 +35,10 @@ int main(int argc, char** argv) { using bout::globals::mesh; Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; + Field3D K{FieldFactory::get()->create3D("K", Options::getRoot(), mesh)}; // Communicate to calculate parallel transform - mesh->communicate(input); + mesh->communicate(input, K); Options dump; // Add mesh geometry variables @@ -45,6 +47,7 @@ int main(int argc, char** argv) { fci_op_test("grad_par", dump, input, Grad_par(input)); fci_op_test("grad2_par2", dump, input, Grad2_par2(input)); fci_op_test("div_par", dump, input, Div_par(input)); + fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index e434ea3ecf..79ec661507 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -9,11 +9,12 @@ from boututils.boutwarnings import AlwaysWarning from boutdata.data import BoutOptionsFile from boutdata.mms import diff, exprToStr, x, y, z -from sympy import sin, sqrt, Expr +from sympy import sin, cos, sqrt, Expr warnings.simplefilter("ignore", AlwaysWarning) f = sin(y - z) + sin(y - 2 * z) +K = cos(z - y) Lx = 0.1 Ly = 10.0 @@ -27,23 +28,33 @@ B = sqrt(Bpx**2 + Bt**2) -def FCI_ddy(f: Expr) -> Expr: +def FCI_grad_par(f: Expr) -> Expr: return (Bt * diff(f, y) * 2.0 * pi / Ly + Bpx * diff(f, z) * 2.0 * pi / Lz) / B +def FCI_grad2_par2(f: Expr) -> Expr: + return FCI_grad_par(FCI_grad_par(f)) + + def FCI_div_par(f: Expr) -> Expr: - return Bpx * FCI_ddy(f / Bpx) + return Bpx * FCI_grad_par(f / Bpx) + + +def FCI_div_par_K_grad_par(f: Expr, K: Expr) -> Expr: + return (K * FCI_grad2_par2(f)) + (FCI_div_par(K) * FCI_grad_par(f)) ############################################ # Equations solved input_field = exprToStr(f) -grad_par_solution = exprToStr(FCI_ddy(f)) -grad2_par2_solution = exprToStr(FCI_ddy(FCI_ddy(f))) +grad_par_solution = exprToStr(FCI_grad_par(f)) +grad2_par2_solution = exprToStr(FCI_grad2_par2(f)) div_par_solution = exprToStr(FCI_div_par(f)) +div_par_K_grad_par_solution = exprToStr(FCI_div_par_K_grad_par(f, K)) print(f"input_field = {input_field}") +print(f"K = {K}") print(f"grad_par_solution = {grad_par_solution}") print(f"grad2_par2_solution = {grad2_par2_solution}") print(f"div_par_solution = {div_par_solution}") @@ -51,7 +62,9 @@ def FCI_div_par(f: Expr) -> Expr: options = BoutOptionsFile("data/BOUT.inp") options["input_field"] = input_field +options["K"] = K options["grad_par_solution"] = grad_par_solution options["grad2_par2_solution"] = grad2_par2_solution options["div_par_solution"] = div_par_solution +options["div_par_K_grad_par_solution"] = div_par_K_grad_par_solution options.write("data/BOUT.inp", overwrite=True) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index c98d924e7b..5a19877a45 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -24,7 +24,7 @@ from scipy.interpolate import RectBivariateSpline as RBS DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ("grad_par", "grad2_par2", "div_par") +OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par") NX = 3 # Resolution in y and z NLIST = [8, 16, 32, 64, 128] @@ -177,7 +177,8 @@ def make_plots(cases): print("No matplotlib") return - fig, axes = plt.subplots(1, len(OPERATORS), figsize=(9, 4)) + num_operators = len(OPERATORS) + fig, axes = plt.subplots(1, num_operators, figsize=(num_operators * 4, 4)) for ax, operator in zip(axes, OPERATORS): for name, case in cases.items(): From 14b4b11a097320b6c936a27fe81b24161cc924ff Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 7 Oct 2025 17:30:48 +0100 Subject: [PATCH 07/39] Many small fixes for FCI - Include all relevant headers, remove unused ones, add some forward declarations - Make data `private` instead of `protected` - Add getter instead of `const` member - Change member reference to pointer - Move ctor to .cxx file - Use `std::array` instead of C array - Bunch of other small clang-tidy fixes --- src/mesh/parallel/fci.cxx | 178 ++++++++++++++++++++++++++------------ src/mesh/parallel/fci.hxx | 83 ++++++------------ 2 files changed, 150 insertions(+), 111 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index e8d3af1cdb..e44a74966e 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -37,79 +37,99 @@ **************************************************************************/ #include "fci.hxx" + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" +#include "bout/field_data.hxx" +#include "bout/mesh.hxx" +#include "bout/msg_stack.hxx" +#include "bout/options.hxx" #include "bout/parallel_boundary_op.hxx" #include "bout/parallel_boundary_region.hxx" -#include -#include -#include -#include -#include +#include "bout/paralleltransform.hxx" +#include "bout/region.hxx" + +#include +#include +#include +#include +#include +#include #include +#include + +using namespace std::string_view_literals; -FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& options, - int offset_, const std::shared_ptr& inner_boundary, +FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, + Options& options, int offset, + const std::shared_ptr& inner_boundary, const std::shared_ptr& outer_boundary, bool zperiodic) - : map_mesh(mesh), offset(offset_), - region_no_boundary(map_mesh.getRegion("RGN_NOBNDRY")), + : map_mesh(&mesh), offset_(offset), + region_no_boundary(map_mesh->getRegion("RGN_NOBNDRY")), corner_boundary_mask(map_mesh) { - TRACE("Creating FCIMAP for direction {:d}", offset); + TRACE("Creating FCIMAP for direction {:d}", offset_); - if (offset == 0) { + if (offset_ == 0) { throw BoutException( "FCIMap called with offset = 0; You probably didn't mean to do that"); } auto& interpolation_options = options["xzinterpolation"]; - interp = - XZInterpolationFactory::getInstance().create(&interpolation_options, &map_mesh); - interp->setYOffset(offset); + interp = XZInterpolationFactory::getInstance().create(&interpolation_options, map_mesh); + interp->setYOffset(offset_); interp_corner = - XZInterpolationFactory::getInstance().create(&interpolation_options, &map_mesh); - interp_corner->setYOffset(offset); + XZInterpolationFactory::getInstance().create(&interpolation_options, map_mesh); + interp_corner->setYOffset(offset_); // Index-space coordinates of forward/backward points - Field3D xt_prime{&map_mesh}, zt_prime{&map_mesh}; + Field3D xt_prime{map_mesh}; + Field3D zt_prime{map_mesh}; // Real-space coordinates of grid points - Field3D R{&map_mesh}, Z{&map_mesh}; + Field3D R{map_mesh}; + Field3D Z{map_mesh}; // Real-space coordinates of forward/backward points - Field3D R_prime{&map_mesh}, Z_prime{&map_mesh}; + Field3D R_prime{map_mesh}; + Field3D Z_prime{map_mesh}; - map_mesh.get(R, "R", 0.0, false); - map_mesh.get(Z, "Z", 0.0, false); + map_mesh->get(R, "R", 0.0, false); + map_mesh->get(Z, "Z", 0.0, false); // Get a unique name for a field based on the sign/magnitude of the offset - const auto parallel_slice_field_name = [&](std::string field) -> std::string { - const std::string direction = (offset > 0) ? "forward" : "backward"; + const auto parallel_slice_field_name = [&](std::string_view field) -> std::string { + const auto direction = (offset_ > 0) ? "forward"sv : "backward"sv; // We only have a suffix for parallel slices beyond the first // This is for backwards compatibility - const std::string slice_suffix = - (std::abs(offset) > 1) ? "_" + std::to_string(std::abs(offset)) : ""; - return direction + "_" + field + slice_suffix; + const auto slice_suffix = + (std::abs(offset_) > 1) ? fmt::format("_{}", std::abs(offset_)) : ""; + return fmt::format("{}_{}{}", direction, field, slice_suffix); }; // If we can't read in any of these fields, things will silently not // work, so best throw - if (map_mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) { + if (map_mesh->get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("xt_prime")); } - if (map_mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) { + if (map_mesh->get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("zt_prime")); } - if (map_mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) { + if (map_mesh->get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("R")); } - if (map_mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) { + if (map_mesh->get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", parallel_slice_field_name("Z")); @@ -157,7 +177,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& interp->calcWeights(xt_prime, zt_prime); } - const int ncz = map_mesh.LocalNz; + const int ncz = map_mesh->LocalNz; BoutMask to_remove(map_mesh); const int xend = @@ -167,15 +187,16 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { // z is periodic, so make sure the z-index wraps around if (zperiodic) { - zt_prime[i] = zt_prime[i] - - ncz * (static_cast(zt_prime[i] / static_cast(ncz))); + zt_prime[i] = + zt_prime[i] + - (ncz * (static_cast(zt_prime[i] / static_cast(ncz)))); if (zt_prime[i] < 0.0) { zt_prime[i] += ncz; } } - if ((xt_prime[i] >= map_mesh.xstart) and (xt_prime[i] <= xend)) { + if ((xt_prime[i] >= map_mesh->xstart) and (xt_prime[i] <= xend)) { // Not a boundary continue; } @@ -215,7 +236,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& const BoutReal dR_dz = 0.5 * (R[i_zp] - R[i_zm]); const BoutReal dZ_dz = 0.5 * (Z[i_zp] - Z[i_zm]); - const BoutReal det = dR_dx * dZ_dz - dR_dz * dZ_dx; // Determinant of 2x2 matrix + const BoutReal det = (dR_dx * dZ_dz) - (dR_dz * dZ_dx); // Determinant of 2x2 matrix const BoutReal dR = R_prime[i] - R[i]; const BoutReal dZ = Z_prime[i] - Z[i]; @@ -228,9 +249,9 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& // outer boundary. However, if any of the surrounding points are negative, // that also means inner. So to differentiate between inner and outer we // need at least 2 points in the domain. - ASSERT2(map_mesh.xend - map_mesh.xstart >= 2); - auto boundary = (xt_prime[i] < map_mesh.xstart) ? inner_boundary : outer_boundary; - boundary->add_point(x, y, z, x + dx, y + 0.5 * offset, + ASSERT2(map_mesh->xend - map_mesh->xstart >= 2); + auto boundary = (xt_prime[i] < map_mesh->xstart) ? inner_boundary : outer_boundary; + boundary->add_point(x, y, z, x + dx, y + (0.5 * offset_), z + dz, // Intersection point in local index space 0.5, // Distance to intersection 1 // Default to that there is a point in the other direction @@ -240,13 +261,14 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& interp->setRegion(region_no_boundary); - const auto region = fmt::format("RGN_YPAR_{:+d}", offset); - if (not map_mesh.hasRegion3D(region)) { + const auto region = fmt::format("RGN_YPAR_{:+d}", offset_); + if (not map_mesh->hasRegion3D(region)) { // The valid region for this slice - map_mesh.addRegion3D( - region, Region(map_mesh.xstart, map_mesh.xend, map_mesh.ystart + offset, - map_mesh.yend + offset, 0, map_mesh.LocalNz - 1, - map_mesh.LocalNy, map_mesh.LocalNz)); + map_mesh->addRegion3D(region, Region(map_mesh->xstart, map_mesh->xend, + map_mesh->ystart + offset_, + map_mesh->yend + offset_, 0, + map_mesh->LocalNz - 1, map_mesh->LocalNy, + map_mesh->LocalNz)); } } @@ -254,7 +276,7 @@ Field3D FCIMap::integrate(Field3D& f) const { TRACE("FCIMap::integrate"); ASSERT1(f.getDirectionY() == YDirectionType::Standard); - ASSERT1(&map_mesh == f.getMesh()); + ASSERT1(map_mesh == f.getMesh()); // Cell centre values Field3D centre = interp->interpolate(f); @@ -269,7 +291,7 @@ Field3D FCIMap::integrate(Field3D& f) const { #endif BOUT_FOR(i, region_no_boundary) { - const auto inext = i.yp(offset); + const auto inext = i.yp(offset_); const BoutReal f_c = centre[inext]; const auto izm = i.zm(); const int x = i.x(); @@ -278,7 +300,7 @@ Field3D FCIMap::integrate(Field3D& f) const { const int zm = izm.z(); if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z) || corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm) - || (x == map_mesh.xstart)) { + || (x == map_mesh->xstart)) { // One of the corners leaves the domain. // Use the cell centre value, since boundary conditions are not // currently applied to corners. @@ -299,10 +321,60 @@ Field3D FCIMap::integrate(Field3D& f) const { return result; } +FCITransform::FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool zperiodic, + Options* opt) + : ParallelTransform(mesh, opt), R{&mesh}, Z{&mesh} { + + // check the coordinate system used for the grid data source + FCITransform::checkInputGrid(); + + // Real-space coordinates of grid cells + mesh.get(R, "R", 0.0, false); + mesh.get(Z, "Z", 0.0, false); + + auto forward_boundary_xin = + std::make_shared("FCI_forward", BNDRY_PAR_FWD_XIN, +1, &mesh); + auto backward_boundary_xin = + std::make_shared("FCI_backward", BNDRY_PAR_BKWD_XIN, -1, &mesh); + auto forward_boundary_xout = + std::make_shared("FCI_forward", BNDRY_PAR_FWD_XOUT, +1, &mesh); + auto backward_boundary_xout = + std::make_shared("FCI_backward", BNDRY_PAR_BKWD_XOUT, -1, &mesh); + + // Add the boundary region to the mesh's vector of parallel boundaries + mesh.addBoundaryPar(forward_boundary_xin, BoundaryParType::xin_fwd); + mesh.addBoundaryPar(backward_boundary_xin, BoundaryParType::xin_bwd); + mesh.addBoundaryPar(forward_boundary_xout, BoundaryParType::xout_fwd); + mesh.addBoundaryPar(backward_boundary_xout, BoundaryParType::xout_bwd); + + field_line_maps.reserve(static_cast(mesh.ystart) * 2); + for (int offset = 1; offset < mesh.ystart + 1; ++offset) { + field_line_maps.emplace_back(mesh, dy, options, offset, forward_boundary_xin, + forward_boundary_xout, zperiodic); + field_line_maps.emplace_back(mesh, dy, options, -offset, backward_boundary_xin, + backward_boundary_xout, zperiodic); + } + ASSERT0(mesh.ystart == 1); + const std::array bndries = {forward_boundary_xin, forward_boundary_xout, + backward_boundary_xin, backward_boundary_xout}; + for (const auto& bndry : bndries) { + for (const auto& bndry2 : bndries) { + if (bndry->dir == bndry2->dir) { + continue; + } + for (bndry->first(); !bndry->isDone(); bndry->next()) { + if (bndry2->contains(*bndry)) { + bndry->setValid(0); + } + } + } + } +} + void FCITransform::checkInputGrid() { std::string parallel_transform; if (mesh.isDataSourceGridFile() - && !mesh.get(parallel_transform, "parallel_transform")) { + && (mesh.get(parallel_transform, "parallel_transform") == 0)) { if (parallel_transform != "fci") { throw BoutException( "Incorrect parallel transform type '" + parallel_transform @@ -310,8 +382,8 @@ void FCITransform::checkInputGrid() { "to generate metric components for FCITransform. Should be 'fci'."); } } // else: parallel_transform variable not found in grid input, indicates older input - // file or grid from options so must rely on the user having ensured the type is - // correct + // file or grid from options so must rely on the user having ensured the type is + // correct } void FCITransform::calcParallelSlices(Field3D& f) { @@ -327,8 +399,8 @@ void FCITransform::calcParallelSlices(Field3D& f) { // Interpolate f onto yup and ydown fields for (const auto& map : field_line_maps) { - f.ynext(map.offset) = map.interpolate(f); - f.ynext(map.offset).setRegion(fmt::format("RGN_YPAR_{:+d}", map.offset)); + f.ynext(map.offset()) = map.interpolate(f); + f.ynext(map.offset()).setRegion(fmt::format("RGN_YPAR_{:+d}", map.offset())); } } @@ -345,7 +417,7 @@ void FCITransform::integrateParallelSlices(Field3D& f) { // Integrate f onto yup and ydown fields for (const auto& map : field_line_maps) { - f.ynext(map.offset) = map.integrate(f); + f.ynext(map.offset()) = map.integrate(f); } } diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 1a02f558e1..65529a4c4e 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -26,6 +26,11 @@ #ifndef BOUT_FCITRANSFORM_H #define BOUT_FCITRANSFORM_H +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/coordinates.hxx" +#include "bout/region.hxx" #include #include #include @@ -33,25 +38,26 @@ #include #include +#include #include +class BoundaryRegionPar; +class FieldPerp; +class Field2D; +class Field3D; +class Options; + /// Field line map - contains the coefficients for interpolation class FCIMap { /// Interpolation objects std::unique_ptr interp; // Cell centre std::unique_ptr interp_corner; // Cell corner at (x+1, z+1) -public: - FCIMap() = delete; - FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, int offset, - const std::shared_ptr& inner_boundary, - const std::shared_ptr& outer_boundary, bool zperiodic); - // The mesh this map was created on - Mesh& map_mesh; + Mesh* map_mesh; /// Direction of map - const int offset; + int offset_; /// region containing all points where the field line has not left the /// domain @@ -59,8 +65,17 @@ public: /// If any of the integration area has left the domain BoutMask corner_boundary_mask; +public: + FCIMap() = delete; + FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, int offset, + const std::shared_ptr& inner_boundary, + const std::shared_ptr& outer_boundary, bool zperiodic); + + /// Direction of map + int offset() const { return offset_; } + Field3D interpolate(Field3D& f) const { - ASSERT1(&map_mesh == f.getMesh()); + ASSERT1(map_mesh == f.getMesh()); return interp->interpolate(f); } @@ -72,55 +87,7 @@ class FCITransform : public ParallelTransform { public: FCITransform() = delete; FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool zperiodic = true, - Options* opt = nullptr) - : ParallelTransform(mesh, opt), R{&mesh}, Z{&mesh} { - - // check the coordinate system used for the grid data source - FCITransform::checkInputGrid(); - - // Real-space coordinates of grid cells - mesh.get(R, "R", 0.0, false); - mesh.get(Z, "Z", 0.0, false); - - auto forward_boundary_xin = - std::make_shared("FCI_forward", BNDRY_PAR_FWD_XIN, +1, &mesh); - auto backward_boundary_xin = std::make_shared( - "FCI_backward", BNDRY_PAR_BKWD_XIN, -1, &mesh); - auto forward_boundary_xout = - std::make_shared("FCI_forward", BNDRY_PAR_FWD_XOUT, +1, &mesh); - auto backward_boundary_xout = std::make_shared( - "FCI_backward", BNDRY_PAR_BKWD_XOUT, -1, &mesh); - - // Add the boundary region to the mesh's vector of parallel boundaries - mesh.addBoundaryPar(forward_boundary_xin, BoundaryParType::xin_fwd); - mesh.addBoundaryPar(backward_boundary_xin, BoundaryParType::xin_bwd); - mesh.addBoundaryPar(forward_boundary_xout, BoundaryParType::xout_fwd); - mesh.addBoundaryPar(backward_boundary_xout, BoundaryParType::xout_bwd); - - field_line_maps.reserve(mesh.ystart * 2); - for (int offset = 1; offset < mesh.ystart + 1; ++offset) { - field_line_maps.emplace_back(mesh, dy, options, offset, forward_boundary_xin, - forward_boundary_xout, zperiodic); - field_line_maps.emplace_back(mesh, dy, options, -offset, backward_boundary_xin, - backward_boundary_xout, zperiodic); - } - ASSERT0(mesh.ystart == 1); - std::shared_ptr bndries[]{ - forward_boundary_xin, forward_boundary_xout, backward_boundary_xin, - backward_boundary_xout}; - for (auto& bndry : bndries) { - for (const auto& bndry2 : bndries) { - if (bndry->dir == bndry2->dir) { - continue; - } - for (bndry->first(); !bndry->isDone(); bndry->next()) { - if (bndry2->contains(*bndry)) { - bndry->setValid(0); - } - } - } - } - } + Options* opt = nullptr); void calcParallelSlices(Field3D& f) override; From da0d4c4a8df2d98f924c87d7f2ac4af290bfb284 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 8 Oct 2025 10:11:22 +0100 Subject: [PATCH 08/39] tests: Fix for 3D metric in FCI test The Div_par operators use parallel slices of B -- with 2D metrics, these are just the field itself, in 3D we need the actual slices. While `Coordinates::geometry` does communicate the fields, it puts off calculating the parallel slices because that needs the fully constructed `Coordinates`. Upcoming changes should fix this and remove the need to explicitly communicate `Coordinates` members. --- src/mesh/coordinates.cxx | 2 +- tests/MMS/spatial/fci/fci_mms.cxx | 9 ++++++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 3dfee6a553..12e465ffb6 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1577,7 +1577,7 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, // Need Bxy at location of f, which might be different from location of this // Coordinates object - auto Bxy_floc = f.getCoordinates()->Bxy; + const auto& Bxy_floc = f.getCoordinates()->Bxy; if (!f.hasParallelSlices()) { // No yup/ydown fields. The Grad_par operator will diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index a0b70e9da6..b30ec05e9a 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -1,4 +1,5 @@ #include "bout/bout.hxx" +#include "bout/build_config.hxx" #include "bout/difops.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" @@ -37,7 +38,13 @@ int main(int argc, char** argv) { Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; Field3D K{FieldFactory::get()->create3D("K", Options::getRoot(), mesh)}; - // Communicate to calculate parallel transform + // Communicate to calculate parallel transform. + if constexpr (bout::build::use_metric_3d) { + // Div_par operators require B parallel slices: + // Coordinates::geometry doesn't ensure this (yet) + auto& Bxy = mesh->getCoordinates()->Bxy; + mesh->communicate(Bxy); + } mesh->communicate(input, K); Options dump; From 2f83692721c64579910175148517b63b9cff7677 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 8 Oct 2025 11:18:59 +0100 Subject: [PATCH 09/39] tests: Add test for FCI `Laplace_par` --- tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 1 + tests/MMS/spatial/fci/mms.py | 7 +++++++ tests/MMS/spatial/fci/runtest | 2 +- 4 files changed, 10 insertions(+), 1 deletion(-) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index e31d84cb35..93e2101473 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -4,6 +4,7 @@ grad2_par2_solution = (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01 div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530717959*cos(y - z) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) K = cos(y - z) +laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index b30ec05e9a..b9d335a3c4 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -55,6 +55,7 @@ int main(int argc, char** argv) { fci_op_test("grad2_par2", dump, input, Grad2_par2(input)); fci_op_test("div_par", dump, input, Div_par(input)); fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); + fci_op_test("laplace_par", dump, input, Laplace_par(input)); bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 79ec661507..178c158572 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -44,6 +44,10 @@ def FCI_div_par_K_grad_par(f: Expr, K: Expr) -> Expr: return (K * FCI_grad2_par2(f)) + (FCI_div_par(K) * FCI_grad_par(f)) +def FCI_Laplace_par(f: Expr) -> Expr: + return FCI_div_par(FCI_grad_par(f)) + + ############################################ # Equations solved @@ -52,6 +56,7 @@ def FCI_div_par_K_grad_par(f: Expr, K: Expr) -> Expr: grad2_par2_solution = exprToStr(FCI_grad2_par2(f)) div_par_solution = exprToStr(FCI_div_par(f)) div_par_K_grad_par_solution = exprToStr(FCI_div_par_K_grad_par(f, K)) +Laplace_par_solution = exprToStr(FCI_Laplace_par(f)) print(f"input_field = {input_field}") print(f"K = {K}") @@ -59,6 +64,7 @@ def FCI_div_par_K_grad_par(f: Expr, K: Expr) -> Expr: print(f"grad2_par2_solution = {grad2_par2_solution}") print(f"div_par_solution = {div_par_solution}") print(f"div_par_K_grad_par_solution = {div_par_K_grad_par_solution}") +print(f"laplace_par_solution = {Laplace_par_solution}") options = BoutOptionsFile("data/BOUT.inp") options["input_field"] = input_field @@ -67,4 +73,5 @@ def FCI_div_par_K_grad_par(f: Expr, K: Expr) -> Expr: options["grad2_par2_solution"] = grad2_par2_solution options["div_par_solution"] = div_par_solution options["div_par_K_grad_par_solution"] = div_par_K_grad_par_solution +options["laplace_par_solution"] = Laplace_par_solution options.write("data/BOUT.inp", overwrite=True) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 5a19877a45..612486bdf4 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -24,7 +24,7 @@ from scipy.interpolate import RectBivariateSpline as RBS DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par") +OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par", "laplace_par") NX = 3 # Resolution in y and z NLIST = [8, 16, 32, 64, 128] From 0cde140be713e61b3ab0636aaa3fb5fccc1218e7 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 8 Oct 2025 11:22:57 +0100 Subject: [PATCH 10/39] Reduce duplication in FCI mms script --- tests/MMS/spatial/fci/mms.py | 36 ++++++++++++++---------------------- 1 file changed, 14 insertions(+), 22 deletions(-) diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 178c158572..b28e337ac0 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -51,27 +51,19 @@ def FCI_Laplace_par(f: Expr) -> Expr: ############################################ # Equations solved -input_field = exprToStr(f) -grad_par_solution = exprToStr(FCI_grad_par(f)) -grad2_par2_solution = exprToStr(FCI_grad2_par2(f)) -div_par_solution = exprToStr(FCI_div_par(f)) -div_par_K_grad_par_solution = exprToStr(FCI_div_par_K_grad_par(f, K)) -Laplace_par_solution = exprToStr(FCI_Laplace_par(f)) - -print(f"input_field = {input_field}") -print(f"K = {K}") -print(f"grad_par_solution = {grad_par_solution}") -print(f"grad2_par2_solution = {grad2_par2_solution}") -print(f"div_par_solution = {div_par_solution}") -print(f"div_par_K_grad_par_solution = {div_par_K_grad_par_solution}") -print(f"laplace_par_solution = {Laplace_par_solution}") - options = BoutOptionsFile("data/BOUT.inp") -options["input_field"] = input_field -options["K"] = K -options["grad_par_solution"] = grad_par_solution -options["grad2_par2_solution"] = grad2_par2_solution -options["div_par_solution"] = div_par_solution -options["div_par_K_grad_par_solution"] = div_par_K_grad_par_solution -options["laplace_par_solution"] = Laplace_par_solution + +for name, expr in ( + ("input_field", f), + ("K", K), + ("grad_par_solution", FCI_grad_par(f)), + ("grad2_par2_solution", FCI_grad2_par2(f)), + ("div_par_solution", FCI_div_par(f)), + ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), + ("laplace_par_solution", FCI_Laplace_par(f)), +): + expr_str = exprToStr(expr) + print(f"{name} = {expr_str}") + options[name] = expr_str + options.write("data/BOUT.inp", overwrite=True) From 37dc27341019c5bf962269e4f015a777014535c4 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 8 Oct 2025 11:25:56 +0100 Subject: [PATCH 11/39] Remove circular include --- include/bout/difops.hxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index 71053d454a..c2fdac195d 100644 --- a/include/bout/difops.hxx +++ b/include/bout/difops.hxx @@ -40,7 +40,9 @@ #include "bout/field3d.hxx" #include "bout/bout_types.hxx" -#include "bout/solver.hxx" +#include "bout/coordinates.hxx" + +class Solver; /*! * Parallel derivative (central differencing) in Y From 9ac46e845b65adb16da3f6e2240bc6cb177636c4 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 8 Oct 2025 16:48:20 +0100 Subject: [PATCH 12/39] tests: Add cases for FCI interpolation methods --- tests/MMS/spatial/fci/runtest | 44 +++++++++++++++-------------------- 1 file changed, 19 insertions(+), 25 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 612486bdf4..6ac0b0bca8 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -99,7 +99,7 @@ def run_fci_operators( ) # Command to run - args = f"MZ={nz} MYG={nslice} mesh:paralleltransform:y_periodic={yperiodic} mesh:ddy:first={name}" + args = f"MZ={nz} MYG={nslice} mesh:paralleltransform:y_periodic={yperiodic} {name}" cmd = f"./fci_mms {args}" print(f"Running command: {cmd}", end="") @@ -138,18 +138,18 @@ def transpose( return result -def check_fci_operators(case: dict) -> bool: +def check_fci_operators(name: str, case: dict) -> bool: failures = [] nslice = case["nslice"] yperiodic = case["yperiodic"] - name = case["name"] order = case["order"] + args = case["args"] all_errors = [] for n in NLIST: - errors = run_fci_operators(nslice, n, yperiodic, name) + errors = run_fci_operators(nslice, n, yperiodic, args) all_errors.append(errors) for operator in OPERATORS: @@ -170,7 +170,7 @@ def check_fci_operators(case: dict) -> bool: return final_errors, failures -def make_plots(cases): +def make_plots(cases: dict[str, dict]): try: import matplotlib.pyplot as plt except ImportError: @@ -199,23 +199,6 @@ def make_plots(cases): plt.close() -def make_case(nslice: int, yperiodic: bool) -> dict[str, Any]: - """ - nslice: - Number of parallel slices (in each direction) - yperiodic: - Run with periodic Y - """ - order = nslice * 2 - return { - "nslice": nslice, - # Which central difference scheme to use and its expected order - "order": order, - "name": f"C{order}", - "yperiodic": yperiodic, - } - - if __name__ == "__main__": build_and_log("FCI MMS test") @@ -240,11 +223,22 @@ if __name__ == "__main__": success = True failures = [] cases = { - "nslice=1": make_case(nslice=1, yperiodic=True), + "nslice=1 hermitespline": { + "nslice": 1, + "order": 2, + "yperiodic": True, + "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=hermitespline", + }, + "nslice=1 lagrange4pt": { + "nslice": 1, + "order": 2, + "yperiodic": True, + "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=lagrange4pt", + }, } - for case in cases.values(): - error2, failures_ = check_fci_operators(case) + for name, case in cases.items(): + error2, failures_ = check_fci_operators(name, case) case.update(error2) failures.extend(failures_) success &= len(failures) == 0 From 5974a6e89d457a8e6cbcdad177837aeda924ad50 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 9 Oct 2025 17:23:45 +0100 Subject: [PATCH 13/39] tests: Increase nx for hermitespline interpolation boundary problem --- src/mesh/interpolation/hermite_spline_xz.cxx | 14 +++++++++++++- tests/MMS/spatial/fci/runtest | 4 +++- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 5020a5b9a3..27a4f1d614 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -171,7 +171,19 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z BoutReal t_x = delta_x(x, y, z) - static_cast(i_corn); BoutReal t_z = delta_z(x, y, z) - static_cast(k_corner(x, y, z)); - // NOTE: A (small) hack to avoid one-sided differences + // NOTE: A (small) hack to avoid one-sided differences. We need at + // least 2 interior points due to an awkwardness with the + // boundaries. The splines need derivatives in x, but we don't + // know the value in the boundaries, so _any_ interpolation in the + // last interior cell can't be done. Instead, we fudge the + // interpolation in the last cell to be at the extreme right-hand + // edge of the previous cell (that is, exactly on the last + // interior point). However, this doesn't work with only one + // interior point, because we have to do something similar to the + // _first_ cell, and these two fudges cancel out and we end up + // indexing into the boundary anyway. + // TODO(peter): Can we remove this if we apply (dirichlet?) BCs to + // the X derivatives? Note that we need at least _2_ if (i_corn >= xend) { i_corn = xend - 1; t_x = 1.0; diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 6ac0b0bca8..5393437ff7 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -25,7 +25,9 @@ DIRECTORY = "data" NPROC = 2 MTHREAD = 2 OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par", "laplace_par") -NX = 3 +# Note that we need at least _2_ interior points for hermite spline +# interpolation due to an awkwardness with the boundaries +NX = 4 # Resolution in y and z NLIST = [8, 16, 32, 64, 128] dx = 1.0 / array(NLIST) From 6ab945e6783b754acf6e84d4fbefcfcd16d94c6f Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 9 Oct 2025 17:24:28 +0100 Subject: [PATCH 14/39] tests: Add monotonichermitespline, decrease resolution scan --- tests/MMS/spatial/fci/runtest | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 5393437ff7..e3d10d989b 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -11,7 +11,6 @@ import json import pathlib import sys from time import time -from typing import Any import boutconfig as conf import zoidberg as zb @@ -29,7 +28,7 @@ OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par", "laplace # interpolation due to an awkwardness with the boundaries NX = 4 # Resolution in y and z -NLIST = [8, 16, 32, 64, 128] +NLIST = [8, 16, 32, 64] dx = 1.0 / array(NLIST) @@ -237,6 +236,12 @@ if __name__ == "__main__": "yperiodic": True, "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=lagrange4pt", }, + "nslice=1 monotonichermitespline": { + "nslice": 1, + "order": 2, + "yperiodic": True, + "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=monotonichermitespline", + }, } for name, case in cases.items(): From 16c3718187378e28c37c5c8304e2d9a5e326d794 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 3 Nov 2025 15:52:25 +0000 Subject: [PATCH 15/39] tests: Small refactor of FCI tests --- .../test-fci-boundary/get_par_bndry.cxx | 22 ++-- tests/integrated/test-fci-boundary/runtest | 36 ++---- tests/integrated/test-fci-mpi/fci_mpi.cxx | 59 +++++----- tests/integrated/test-fci-mpi/runtest | 107 +++++++++++------- 4 files changed, 116 insertions(+), 108 deletions(-) diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index ac0f5de2a6..ba282d8988 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -1,5 +1,5 @@ #include "bout/bout.hxx" -#include "bout/derivs.hxx" +#include "bout/field3d.hxx" #include "bout/field_factory.hxx" #include "bout/parallel_boundary_region.hxx" @@ -8,24 +8,24 @@ int main(int argc, char** argv) { using bout::globals::mesh; - std::vector fields; - fields.resize(static_cast(BoundaryParType::SIZE)); + std::vector fields (static_cast(BoundaryParType::SIZE), Field3D{0.0}); + Options dump; for (int i = 0; i < fields.size(); i++) { - fields[i] = Field3D{0.0}; + fields[i].allocate(); + const auto boundary = static_cast(i); + const auto boundary_name = toString(boundary); mesh->communicate(fields[i]); - for (const auto& bndry_par : - mesh->getBoundariesPar(static_cast(i))) { - output.write("{:s} region\n", toString(static_cast(i))); + for (const auto& bndry_par : mesh->getBoundariesPar(boundary)) { + output.write("{:s} region\n", boundary_name); for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { fields[i][bndry_par->ind()] += 1; - output.write("{:s} increment\n", toString(static_cast(i))); + output.write("{:s} increment\n", boundary_name); } } - output.write("{:s} done\n", toString(static_cast(i))); + output.write("{:s} done\n", boundary_name); - dump[fmt::format("field_{:s}", toString(static_cast(i)))] = - fields[i]; + dump[fmt::format("field_{:s}", boundary_name)] = fields[i]; } bout::writeDefaultOutputFile(dump); diff --git a/tests/integrated/test-fci-boundary/runtest b/tests/integrated/test-fci-boundary/runtest index 1b1460da53..e749055185 100755 --- a/tests/integrated/test-fci-boundary/runtest +++ b/tests/integrated/test-fci-boundary/runtest @@ -1,29 +1,15 @@ #!/usr/bin/env python3 # # Python script to run and analyse MMS test -# -# Cores: 2 -# only working with cmake -# requires: False from boututils.run_wrapper import launch_safe from boututils.datafile import DataFile -from boutdata.collect import collect as _collect +from boutdata.collect import collect import numpy as np -def collect(var): - return _collect( - var, - info=False, - path=directory, - xguards=False, - yguards=False, - ) - - -nprocs = [1] # , 2, 4] +nprocs = [1] mthread = 2 directory = "data" @@ -43,11 +29,6 @@ regions = { } regions = {k: v.astype(int) for k, v in regions.items()} -# for x in "xout", "xin": -# regions[x] = np.logical_or(regions[f"{x}_fwd"], regions[f"{x}_bwd"]) -# for x in "fwd", "bwd": -# regions[x] = np.logical_or(regions[f"xin_{x}"], regions[f"xout_{x}"]) -# regions["all"] = np.logical_or(regions["xin"], regions["xout"]) for x in "xout", "xin": regions[x] = regions[f"{x}_fwd"] + regions[f"{x}_bwd"] for x in "fwd", "bwd": @@ -56,15 +37,18 @@ regions["all"] = regions["xin"] + regions["xout"] for nproc in nprocs: cmd = "./get_par_bndry" - - # Launch using MPI _, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) for k, v in regions.items(): - # Collect data - data = collect(f"field_{k}") + data = collect( + f"field_{k}", + info=False, + path=directory, + xguards=False, + yguards=False, + ) assert np.allclose(data, v), ( - k + " does not match", + f"{k} does not match", np.sum(data), np.sum(v), np.max(data), diff --git a/tests/integrated/test-fci-mpi/fci_mpi.cxx b/tests/integrated/test-fci-mpi/fci_mpi.cxx index 94520dd4a6..cc4fba8ffe 100644 --- a/tests/integrated/test-fci-mpi/fci_mpi.cxx +++ b/tests/integrated/test-fci-mpi/fci_mpi.cxx @@ -1,38 +1,37 @@ +#include "fmt/format.h" #include "bout/bout.hxx" -#include "bout/derivs.hxx" #include "bout/field_factory.hxx" +namespace { +auto fci_mpi_test(int num, Options& dump) { + using bout::globals::mesh; + Field3D input{FieldFactory::get()->create3D(fmt::format("input_{:d}:function", num), + Options::getRoot(), mesh)}; + mesh->communicate(input); + + input.applyParallelBoundary("parallel_neumann_o2"); + + for (int slice = -mesh->ystart; slice <= mesh->ystart; ++slice) { + if (slice == 0) { + continue; + } + Field3D tmp{0.}; + BOUT_FOR(i, tmp.getRegion("RGN_NOBNDRY")) { + tmp[i] = input.ynext(slice)[i.yp(slice)]; + } + dump[fmt::format("output_{:d}_{:+d}", num, slice)] = tmp; + } +} +} // namespace + int main(int argc, char** argv) { BoutInitialise(argc, argv); - { - using bout::globals::mesh; - Options* options = Options::getRoot(); - int i = 0; - const std::string default_str{"not_set"}; - Options dump; - while (true) { - std::string temp_str; - options->get(fmt::format("input_{:d}:function", i), temp_str, default_str); - if (temp_str == default_str) { - break; - } - Field3D input{FieldFactory::get()->create3D(fmt::format("input_{:d}:function", i), - Options::getRoot(), mesh)}; - // options->get(fmt::format("input_{:d}:boundary_perp", i), temp_str, s"free_o3"); - mesh->communicate(input); - input.applyParallelBoundary("parallel_neumann_o2"); - for (int slice = -mesh->ystart; slice <= mesh->ystart; ++slice) { - if (slice != 0) { - Field3D tmp{0.}; - BOUT_FOR(i, tmp.getRegion("RGN_NOBNDRY")) { - tmp[i] = input.ynext(slice)[i.yp(slice)]; - } - dump[fmt::format("output_{:d}_{:+d}", i, slice)] = tmp; - } - } - ++i; - } - bout::writeDefaultOutputFile(dump); + Options dump; + + for (auto num : {0, 1, 2, 3}) { + fci_mpi_test(num, dump); } + + bout::writeDefaultOutputFile(dump); BoutFinalise(); } diff --git a/tests/integrated/test-fci-mpi/runtest b/tests/integrated/test-fci-mpi/runtest index 6676f8f7a5..828d8d4a50 100755 --- a/tests/integrated/test-fci-mpi/runtest +++ b/tests/integrated/test-fci-mpi/runtest @@ -1,57 +1,82 @@ #!/usr/bin/env python3 # # Python script to run and analyse MMS test -# - -# Cores: 8 -# requires: metric_3d -from boututils.run_wrapper import build_and_log, launch_safe, shell_safe +from boututils.run_wrapper import build_and_log, launch_safe from boutdata.collect import collect -import boutconfig as conf import itertools +import sys -import numpy as np +import numpy.testing as npt # Resolution in x and y -nlist = [1, 2, 4] +NLIST = [1, 2, 4] +MAXCORES = 8 +NSLICES = [1] -maxcores = 8 +build_and_log("FCI MMS test") -nslices = [1] +COLLECT_KW = dict(info=False, xguards=False, yguards=False, path="data") -success = True -build_and_log("FCI MMS test") +def run_case(nxpe: int, nype: int, mthread: int): + cmd = f"./fci_mpi NXPE={nxpe} NYPE={nype}" + print(f"Running command: {cmd}") + + _, out = launch_safe(cmd, nproc=nxpe * nype, mthread=mthread, pipe=True) + + # Save output to log file + with open(f"run.log.{nxpe}.{nype}.{nslice}.log", "w") as f: + f.write(out) + + +def test_case(nxpe: int, nype: int, mthread: int, ref: dict) -> bool: + run_case(nxpe, nype, mthread) + + failures = [] + + for name, val in ref.items(): + try: + npt.assert_allclose(val, collect(name, **COLLECT_KW)) + except AssertionError as e: + failures.append((nxpe, nype, name, e)) -for nslice in nslices: - for NXPE, NYPE in itertools.product(nlist, nlist): - if NXPE * NYPE > maxcores: + return failures + + +failures = [] + +for nslice in NSLICES: + # reference data! + run_case(1, 1, MAXCORES) + + ref = {} + for i in range(4): + for yp in range(1, nslice + 1): + for y in [-yp, yp]: + name = f"output_{i}_{y:+d}" + ref[name] = collect(name, **COLLECT_KW) + + for nxpe, nype in itertools.product(NLIST, NLIST): + if (nxpe, nype) == (1, 1): + # reference case, done above continue - args = f"NXPE={NXPE} NYPE={NYPE}" - # Command to run - cmd = f"./fci_mpi {args}" - - print(f"Running command: {cmd}") - - mthread = maxcores // (NXPE * NYPE) - # Launch using MPI - _, out = launch_safe(cmd, nproc=NXPE * NYPE, mthread=mthread, pipe=True) - - # Save output to log file - with open(f"run.log.{NXPE}.{NYPE}.{nslice}.log", "w") as f: - f.write(out) - - collect_kw = dict(info=False, xguards=False, yguards=False, path="data") - if NXPE == NYPE == 1: - # reference data! - ref = {} - for i in range(4): - for yp in range(1, nslice + 1): - for y in [-yp, yp]: - name = f"output_{i}_{y:+d}" - ref[name] = collect(name, **collect_kw) - else: - for name, val in ref.items(): - assert np.allclose(val, collect(name, **collect_kw)) + if nxpe * nype > MAXCORES: + continue + + mthread = MAXCORES // (nxpe * nype) + failures_ = test_case(nxpe, nype, mthread, ref) + failures.extend(failures_) + + +success = len(failures) == 0 +if success: + print("\nAll tests passed") +else: + print("\nSome tests failed:") + for (nxpe, nype, name, error) in failures: + print("----------") + print(f"case {nxpe=} {nype=} {name=}\n{error}") + +sys.exit(0 if success else 1) From 3d3ff654b60308e249ae8410a849bfec728f26f1 Mon Sep 17 00:00:00 2001 From: ZedThree <1486942+ZedThree@users.noreply.github.com> Date: Mon, 3 Nov 2025 17:05:43 +0000 Subject: [PATCH 16/39] Apply black changes --- tests/integrated/test-fci-mpi/runtest | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrated/test-fci-mpi/runtest b/tests/integrated/test-fci-mpi/runtest index 828d8d4a50..c18ab0391d 100755 --- a/tests/integrated/test-fci-mpi/runtest +++ b/tests/integrated/test-fci-mpi/runtest @@ -75,7 +75,7 @@ if success: print("\nAll tests passed") else: print("\nSome tests failed:") - for (nxpe, nype, name, error) in failures: + for nxpe, nype, name, error in failures: print("----------") print(f"case {nxpe=} {nype=} {name=}\n{error}") From 79f9d0fe70c00bc578404351e8379f1fa79a3dc0 Mon Sep 17 00:00:00 2001 From: ZedThree <1486942+ZedThree@users.noreply.github.com> Date: Tue, 4 Nov 2025 11:58:36 +0000 Subject: [PATCH 17/39] Apply clang-format changes --- src/mesh/parallel/fci.cxx | 4 ++-- tests/integrated/test-fci-boundary/get_par_bndry.cxx | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index e44a74966e..7c7ade1c8d 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -180,8 +180,8 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, const int ncz = map_mesh->LocalNz; BoutMask to_remove(map_mesh); - const int xend = - map_mesh.xstart + (map_mesh.xend - map_mesh.xstart + 1) * map_mesh.getNXPE() - 1; + const int xend = map_mesh->xstart + + ((map_mesh->xend - map_mesh->xstart + 1) * map_mesh->getNXPE()) - 1; // Serial loop because call to BoundaryRegionPar::addPoint // (probably?) can't be done in parallel BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index ba282d8988..8183d989d1 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -8,7 +8,7 @@ int main(int argc, char** argv) { using bout::globals::mesh; - std::vector fields (static_cast(BoundaryParType::SIZE), Field3D{0.0}); + std::vector fields(static_cast(BoundaryParType::SIZE), Field3D{0.0}); Options dump; for (int i = 0; i < fields.size(); i++) { From cee3e445cb7da3d7ce24cbae32cd5ccf53a9509c Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 4 Nov 2025 16:13:29 +0000 Subject: [PATCH 18/39] Apply clang-tidy fixes to FCI tests --- tests/MMS/spatial/fci/fci_mms.cxx | 13 ++++++++++--- .../integrated/test-fci-boundary/get_par_bndry.cxx | 8 ++++++++ tests/integrated/test-fci-mpi/fci_mpi.cxx | 6 +++++- 3 files changed, 23 insertions(+), 4 deletions(-) diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index b9d335a3c4..7967452f3d 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -1,19 +1,26 @@ #include "bout/bout.hxx" #include "bout/build_config.hxx" #include "bout/difops.hxx" +#include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/globals.hxx" #include "bout/options.hxx" +#include "bout/options_io.hxx" +#include "bout/utils.hxx" +#include + +#include #include namespace { auto fci_op_test(const std::string& name, Options& dump, const Field3D& input, const Field3D& result) { auto* mesh = input.getMesh(); - Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), - Options::getRoot(), mesh)}; - Field3D error{result - solution}; + const Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), + Options::getRoot(), mesh)}; + const Field3D error{result - solution}; dump[fmt::format("{}_l_2", name)] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); dump[fmt::format("{}_l_inf", name)] = max(abs(error), true, "RGN_NOBNDRY"); diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index 8183d989d1..8e3cfac2f7 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -1,8 +1,16 @@ #include "bout/bout.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" +#include "bout/output.hxx" #include "bout/parallel_boundary_region.hxx" +#include + +#include + int main(int argc, char** argv) { BoutInitialise(argc, argv); diff --git a/tests/integrated/test-fci-mpi/fci_mpi.cxx b/tests/integrated/test-fci-mpi/fci_mpi.cxx index cc4fba8ffe..94e8e878ef 100644 --- a/tests/integrated/test-fci-mpi/fci_mpi.cxx +++ b/tests/integrated/test-fci-mpi/fci_mpi.cxx @@ -1,6 +1,11 @@ #include "fmt/format.h" #include "bout/bout.hxx" +#include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" +#include "bout/region.hxx" namespace { auto fci_mpi_test(int num, Options& dump) { @@ -8,7 +13,6 @@ auto fci_mpi_test(int num, Options& dump) { Field3D input{FieldFactory::get()->create3D(fmt::format("input_{:d}:function", num), Options::getRoot(), mesh)}; mesh->communicate(input); - input.applyParallelBoundary("parallel_neumann_o2"); for (int slice = -mesh->ystart; slice <= mesh->ystart; ++slice) { From 88d8db852528f91dffe9d16e93cc3f107baadebb Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 4 Nov 2025 16:14:55 +0000 Subject: [PATCH 19/39] Fix comment formatting --- src/mesh/parallel/fci.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 7c7ade1c8d..055e239725 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -381,7 +381,8 @@ void FCITransform::checkInputGrid() { + "' used " "to generate metric components for FCITransform. Should be 'fci'."); } - } // else: parallel_transform variable not found in grid input, indicates older input + } + // else: parallel_transform variable not found in grid input, indicates older input // file or grid from options so must rely on the user having ensured the type is // correct } From 9d9922c522a1fbe377377e9498ea839b9f7ce82d Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 14:27:05 +0000 Subject: [PATCH 20/39] Apply suggestion from @dschwoerer Co-authored-by: David Bold --- src/mesh/parallel/fci.cxx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 055e239725..6243bbb67a 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -107,9 +107,10 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, const auto direction = (offset_ > 0) ? "forward"sv : "backward"sv; // We only have a suffix for parallel slices beyond the first // This is for backwards compatibility - const auto slice_suffix = - (std::abs(offset_) > 1) ? fmt::format("_{}", std::abs(offset_)) : ""; - return fmt::format("{}_{}{}", direction, field, slice_suffix); + if (std::abs(offset_) == 1) { + return fmt::format("{}_{}", direction, field); + } + return fmt::format("{}_{}_{}", direction, field, std::abs(offset_)); }; // If we can't read in any of these fields, things will silently not From e961e75696aa7fcaf194c6994732725df9bc900c Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 17:13:18 +0000 Subject: [PATCH 21/39] tests: Fix typo that made tests always pass --- tests/MMS/spatial/fci/runtest | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index e3d10d989b..b1321c0c5d 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -55,7 +55,7 @@ def quiet_collect(name: str) -> float: )[()] -def assert_convergence(error, dx, name, order) -> bool: +def assert_convergence(error, dx, name, expected) -> bool: fit = polyfit(log(dx), log(error), 1) order = fit[0] print(f"{name} convergence order = {order:f} (fit)", end="") @@ -64,7 +64,7 @@ def assert_convergence(error, dx, name, order) -> bool: print(f", {order:f} (small spacing)", end="") # Should be close to the expected order - success = order > order * 0.95 + success = order > expected * 0.95 print(f"\t............ {'PASS' if success else 'FAIL'}") return success From 2674a58d4e77fda0db5bc213164b245c2854dc8b Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 17:28:58 +0000 Subject: [PATCH 22/39] tests: Remove monotonichermitespline FCI case --- tests/MMS/spatial/fci/runtest | 6 ------ tests/integrated/test-interpolate/runtest | 1 + 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index b1321c0c5d..2bfde26f51 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -236,12 +236,6 @@ if __name__ == "__main__": "yperiodic": True, "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=lagrange4pt", }, - "nslice=1 monotonichermitespline": { - "nslice": 1, - "order": 2, - "yperiodic": True, - "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=monotonichermitespline", - }, } for name, case in cases.items(): diff --git a/tests/integrated/test-interpolate/runtest b/tests/integrated/test-interpolate/runtest index f5460aff2a..da10472610 100755 --- a/tests/integrated/test-interpolate/runtest +++ b/tests/integrated/test-interpolate/runtest @@ -25,6 +25,7 @@ methods = { "hermitespline": 3, "lagrange4pt": 3, "bilinear": 2, + "monotonichermitespline": 2, } From f00db393cf180eba7a44136ada78483e30f965eb Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 21 Nov 2025 15:23:30 +0000 Subject: [PATCH 23/39] Add `rtol`, `atol` to `MonotonicHermiteSpline` This allows the user to control the clipping, allowing some overshoot Also restore the FCI MMS test with this interpolator --- include/bout/interpolation_xz.hxx | 40 +++++++++++++++---- .../monotonic_hermite_spline_xz.cxx | 13 +++--- tests/MMS/spatial/fci/runtest | 11 +++++ tests/integrated/test-interpolate/runtest | 1 - 4 files changed, 48 insertions(+), 17 deletions(-) diff --git a/include/bout/interpolation_xz.hxx b/include/bout/interpolation_xz.hxx index 6c7419f7e4..d12eecc3d3 100644 --- a/include/bout/interpolation_xz.hxx +++ b/include/bout/interpolation_xz.hxx @@ -24,6 +24,7 @@ #ifndef BOUT_INTERP_XZ_H #define BOUT_INTERP_XZ_H +#include "bout/bout_types.hxx" #include "bout/mask.hxx" #define USE_NEW_WEIGHTS 1 @@ -166,7 +167,8 @@ protected: #endif public: - XZHermiteSpline(Mesh* mesh = nullptr) : XZHermiteSpline(0, mesh) {} + XZHermiteSpline(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr) + : XZHermiteSpline(0, mesh) {} XZHermiteSpline(int y_offset = 0, Mesh* mesh = nullptr); XZHermiteSpline(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) : XZHermiteSpline(y_offset, mesh) { @@ -210,9 +212,29 @@ public: /// but also degrades accuracy near maxima and minima. /// Perhaps should only impose near boundaries, since that is where /// problems most obviously occur. +/// +/// You can control how tight the clipping to the range of the neighbouring cell +/// values through ``rtol`` and ``atol``: +/// +/// diff = (max_of_neighours - min_of_neighours) * rtol + atol +/// +/// and the interpolated value is instead clipped to the range +/// ``[min_of_neighours - diff, max_of_neighours + diff]`` class XZMonotonicHermiteSpline : public XZHermiteSpline { + /// Absolute tolerance for clipping + BoutReal atol = 0.0; + /// Relative tolerance for clipping + BoutReal rtol = 1.0; + public: - XZMonotonicHermiteSpline(Mesh* mesh = nullptr) : XZHermiteSpline(0, mesh) { + XZMonotonicHermiteSpline(Mesh* mesh = nullptr, Options* options = nullptr) + : XZHermiteSpline(0, mesh), + atol{(*options)["atol"] + .doc("Absolute tolerance for clipping overshoot") + .withDefault(0.0)}, + rtol{(*options)["rtol"] + .doc("Relative tolerance for clipping overshoot") + .withDefault(1.0)} { if (localmesh->getNXPE() > 1) { throw BoutException("Do not support MPI splitting in X"); } @@ -248,7 +270,8 @@ class XZLagrange4pt : public XZInterpolation { Field3D t_x, t_z; public: - XZLagrange4pt(Mesh* mesh = nullptr) : XZLagrange4pt(0, mesh) {} + XZLagrange4pt(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr) + : XZLagrange4pt(0, mesh) {} XZLagrange4pt(int y_offset = 0, Mesh* mesh = nullptr); XZLagrange4pt(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) : XZLagrange4pt(y_offset, mesh) { @@ -284,7 +307,8 @@ class XZBilinear : public XZInterpolation { Field3D w0, w1, w2, w3; public: - XZBilinear(Mesh* mesh = nullptr) : XZBilinear(0, mesh) {} + XZBilinear(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr) + : XZBilinear(0, mesh) {} XZBilinear(int y_offset = 0, Mesh* mesh = nullptr); XZBilinear(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) : XZBilinear(y_offset, mesh) { @@ -308,7 +332,7 @@ public: }; class XZInterpolationFactory - : public Factory { + : public Factory { public: static constexpr auto type_name = "XZInterpolation"; static constexpr auto section_name = "xzinterpolation"; @@ -316,10 +340,10 @@ public: static constexpr auto default_type = "hermitespline"; ReturnType create(Options* options = nullptr, Mesh* mesh = nullptr) const { - return Factory::create(getType(options), mesh); + return Factory::create(getType(options), mesh, options); } - ReturnType create(const std::string& type, [[maybe_unused]] Options* options) const { - return Factory::create(type, nullptr); + ReturnType create(const std::string& type, Options* options) const { + return Factory::create(type, nullptr, options); } static void ensureRegistered(); diff --git a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx index f23bfd499e..f206ed1e0f 100644 --- a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx +++ b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx @@ -25,7 +25,7 @@ #include "bout/interpolation_xz.hxx" #include "bout/mesh.hxx" -#include +#include Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, const std::string& region) const { @@ -80,7 +80,6 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, // Perhaps should only impose near boundaries, since that is where // problems most obviously occur. const BoutReal localmax = BOUTMAX(f[ic], f[icxp], f[iczp], f[icxpzp]); - const BoutReal localmin = BOUTMIN(f[ic], f[icxp], f[iczp], f[icxpzp]); ASSERT2(std::isfinite(localmax) || i.x() < localmesh->xstart @@ -88,12 +87,10 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, ASSERT2(std::isfinite(localmin) || i.x() < localmesh->xstart || i.x() > localmesh->xend); - if (result > localmax) { - result = localmax; - } - if (result < localmin) { - result = localmin; - } + const auto diff = ((localmax - localmin) * rtol) + atol; + + result = std::min(result, localmax + diff); + result = std::max(result, localmin - diff); f_interp[iyp] = result; } diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 2bfde26f51..34340e53f4 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -236,6 +236,17 @@ if __name__ == "__main__": "yperiodic": True, "args": "mesh:ddy:first=C2 mesh:paralleltransform:xzinterpolation:type=lagrange4pt", }, + "nslice=1 monotonichermitespline": { + "nslice": 1, + "order": 2, + "yperiodic": True, + "args": ( + "mesh:ddy:first=C2 " + "mesh:paralleltransform:xzinterpolation:type=monotonichermitespline " + "mesh:paralleltransform:xzinterpolation:rtol=1e-3 " + "mesh:paralleltransform:xzinterpolation:atol=5e-3" + ), + }, } for name, case in cases.items(): diff --git a/tests/integrated/test-interpolate/runtest b/tests/integrated/test-interpolate/runtest index da10472610..f5460aff2a 100755 --- a/tests/integrated/test-interpolate/runtest +++ b/tests/integrated/test-interpolate/runtest @@ -25,7 +25,6 @@ methods = { "hermitespline": 3, "lagrange4pt": 3, "bilinear": 2, - "monotonichermitespline": 2, } From 0efe4e3d68486f8a7282ce3d6d695876867ca3c9 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 3 Dec 2025 14:34:36 +0000 Subject: [PATCH 24/39] Add missing header --- include/bout/interpolation_xz.hxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/bout/interpolation_xz.hxx b/include/bout/interpolation_xz.hxx index d12eecc3d3..4dd24259fd 100644 --- a/include/bout/interpolation_xz.hxx +++ b/include/bout/interpolation_xz.hxx @@ -24,8 +24,9 @@ #ifndef BOUT_INTERP_XZ_H #define BOUT_INTERP_XZ_H -#include "bout/bout_types.hxx" -#include "bout/mask.hxx" +#include +#include +#include #define USE_NEW_WEIGHTS 1 #if BOUT_HAS_PETSC From e44fbb5af3b47cde3fe08a744fe6a3223565359d Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 15:09:29 +0000 Subject: [PATCH 25/39] Port `FV::Div_par_mod` from Hermes-3 --- include/bout/fv_ops.hxx | 253 +++++++++++++++++++++++++++- tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 5 + tests/MMS/spatial/fci/mms.py | 2 + tests/MMS/spatial/fci/runtest | 9 +- 5 files changed, 266 insertions(+), 4 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 94007a57a2..1519313fc7 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -5,12 +5,16 @@ #ifndef BOUT_FV_OPS_H #define BOUT_FV_OPS_H +#include + +#include "bout/bout_types.hxx" #include "bout/field3d.hxx" #include "bout/globals.hxx" -#include "bout/vector2d.hxx" - +#include "bout/mesh.hxx" +#include "bout/output_bout_types.hxx" // NOLINT(unused-includes) +#include "bout/region.hxx" #include "bout/utils.hxx" -#include +#include "bout/vector2d.hxx" namespace FV { /*! @@ -524,5 +528,248 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { * X-Z Finite Volume diffusion operator */ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); + +/// Finite volume parallel divergence +/// +/// NOTE: Modified version, applies limiter to velocity and field +/// Performs better (smaller overshoots) than Div_par +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// @param[out] flow_ylow Flow at the lower Y cell boundary +/// Already includes area factor * flux +template +Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true) { + + Coordinates* coord = f_in.getCoordinates(); + + if (f_in.isFci()) { + // Use mid-point (cell boundary) averages + if (flow_ylow.isAllocated()) { + flow_ylow = emptyFrom(flow_ylow); + } + + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) + * (coord->J[i] + coord->J.yup()[iyp]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp])) + - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) + * (coord->J[i] + coord->J.ydown()[iym]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym]))) + / (coord->dy[i] * coord->J[i]); + } + return result; + } + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); + + const Mesh* mesh = f_in.getMesh(); + + CellEdges cellboundary; + + ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); + ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); + const bool are_unaligned = + ((f_in.getDirectionY() == YDirectionType::Standard) + and (v_in.getDirectionY() == YDirectionType::Standard) + and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); + + const Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; + const Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; + const Field3D wave_speed = + are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes +#if not(BOUT_USE_METRIC_3D) + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j) + coord->J(i, j + 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j) + coord->J(i, j - 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1)); + + const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j); +#endif + for (int k = 0; k < mesh->LocalNz; k++) { +#if BOUT_USE_METRIC_3D + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + const BoutReal area_lc = + common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k); +#endif + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); // Calculate sv.R and sv.L + + //////////////////////////////////////////// + // Right boundary + + BoutReal flux = BoutNaN; + + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { + // Last point in domain + + // Calculate velocity at right boundary (y+1/2) + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); + + const BoutReal bndryval = 0.5 * (s.c + s.p); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); + } + + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(v(i, j, k)), fabs(v(i, j + 1, k))); + + flux = s.R * 0.5 * (sv.R + amax); + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + flow_ylow(i, j + 1, k) += flux * area_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { + // First point in domain + const BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(v(i, j, k)), fabs(v(i, j - 1, k))); + + flux = s.L * 0.5 * (sv.L - amax); + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + + flow_ylow(i, j, k) += flux * area_lc; + } + } + } + if (are_unaligned) { + flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); + } + return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; +} } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 93e2101473..9171178d5d 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -5,6 +5,7 @@ div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530 div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 7967452f3d..ca8ee8cd9f 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -4,6 +4,7 @@ #include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" #include "bout/globals.hxx" #include "bout/options.hxx" #include "bout/options_io.hxx" @@ -64,6 +65,10 @@ int main(int argc, char** argv) { fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); fci_op_test("laplace_par", dump, input, Laplace_par(input)); + // Finite volume methods + Field3D flow_ylow; + fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + bout::writeDefaultOutputFile(dump); BoutFinalise(); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index b28e337ac0..2fb2bd6aa6 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -16,6 +16,7 @@ f = sin(y - z) + sin(y - 2 * z) K = cos(z - y) + Lx = 0.1 Ly = 10.0 Lz = 1.0 @@ -61,6 +62,7 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("div_par_solution", FCI_div_par(f)), ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), + ("FV_div_par_mod_solution", FCI_div_par(f * K)), ): expr_str = exprToStr(expr) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 34340e53f4..7e960cb024 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -23,7 +23,14 @@ from scipy.interpolate import RectBivariateSpline as RBS DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par", "laplace_par") +OPERATORS = ( + "grad_par", + "grad2_par2", + "div_par", + "div_par_K_grad_par", + "laplace_par", + "FV_div_par_mod", +) # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries NX = 4 From db5027890888fca325ac39dccdc9d5da1b447396 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 15:43:38 +0000 Subject: [PATCH 26/39] Port `FV::Div_par_fvv` from Hermes-3 --- include/bout/fv_ops.hxx | 208 ++++++++++++++++++++++++++++ tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 1 + tests/MMS/spatial/fci/mms.py | 1 + tests/MMS/spatial/fci/runtest | 1 + 5 files changed, 212 insertions(+) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 1519313fc7..0960476911 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -771,5 +771,213 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, } return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } + +/// This operator calculates Div_par(f v v) +/// It is used primarily (only?) in the parallel momentum equation. +/// +/// This operator is used rather than Div(f fv) so that the values of +/// f and v are consistent with other advection equations: The product +/// fv is not interpolated to cell boundaries. +template +Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true) { + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + const Mesh* mesh = f_in.getMesh(); + const Coordinates* coord = f_in.getCoordinates(); + CellEdges cellboundary; + + if (f_in.isFci()) { + // FCI version, using yup/down fields + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& B = coord->Bxy; + const auto& B_up = coord->Bxy.yup(); + const auto& B_down = coord->Bxy.ydown(); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + const auto& g_22 = coord->g_22; + const auto& dy = coord->dy; + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Maximum local wave speed + const BoutReal amax = + BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); + + const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) + - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]); + + // Penalty terms. This implementation is very dissipative. + BoutReal penalty = + (amax * (f_in[i] * v_in[i] - f_up[iyp] * v_up[iyp]) / (B[i] + B_up[iyp])) + + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym]) + / (B[i] + B_down[iym])); + + if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) { + if (term * penalty > 0) { + penalty = term; + } else { + penalty = -term; + } + } + + result[i] = B[i] * (term + penalty) / (2 * dy[i] * sqrt(g_22[i])); + +#if CHECK > 0 + if (!std::isfinite(result[i])) { + throw BoutException("Non-finite value in Div_par_fvv at {}\n" + "fup {} vup {} fdown {} vdown {} amax {}\n", + "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], + v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], + dy[i], sqrt(g_22[i])); + } +#endif + } + return result; + } + + ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); + + /// Ensure that f, v and wave_speed are field aligned + Field3D f = toFieldAligned(f_in, "RGN_NOX"); + Field3D v = toFieldAligned(v_in, "RGN_NOX"); + Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes + + for (int k = 0; k < mesh->LocalNz; k++) { + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); + + //////////////////////////////////////////// + // Right boundary + + // Calculate velocity at right boundary (y+1/2) + const BoutReal v_mid_r = 0.5 * (sv.c + sv.p); + // And mid-point density at right boundary + const BoutReal n_mid_r = 0.5 * (s.c + s.p); + BoutReal flux = NAN; + + if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + // Last point in domain + + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_r * v_mid_r * v_mid_r; + } else { + // Add flux due to difference in boundary values + flux = (s.R * sv.R * sv.R) // Use right cell edge values + + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid_r + * (sv.R - v_mid_r)); // Damp differences in velocity, not flux + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(sv.c), fabs(sv.p)); + + flux = s.R * 0.5 * (sv.R + amax) * sv.R; + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + const BoutReal v_mid_l = 0.5 * (sv.c + sv.m); + const BoutReal n_mid_l = 0.5 * (s.c + s.m); + + if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + // First point in domain + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_l * v_mid_l * v_mid_l; + } else { + // Add flux due to difference in boundary values + flux = (s.L * sv.L * sv.L) + - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid_l + * (sv.L - v_mid_l)); + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(sv.c), fabs(sv.m)); + + flux = s.L * 0.5 * (sv.L - amax) * sv.L; + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + } + } + } + return fromFieldAligned(result, "RGN_NOBNDRY"); +} } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 9171178d5d..99edee6ecc 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -6,6 +6,7 @@ div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0. K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_fvv_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(2*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)^2/(0.01*x + 0.045)) - 1.25663706143592*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)^2/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index ca8ee8cd9f..13744c4965 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -68,6 +68,7 @@ int main(int argc, char** argv) { // Finite volume methods Field3D flow_ylow; fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + fci_op_test("FV_div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 2fb2bd6aa6..62089c7e21 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -63,6 +63,7 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), ("FV_div_par_mod_solution", FCI_div_par(f * K)), + ("FV_div_par_fvv_solution", FCI_div_par(f * K * K)), ): expr_str = exprToStr(expr) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 7e960cb024..c0f0a45132 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -30,6 +30,7 @@ OPERATORS = ( "div_par_K_grad_par", "laplace_par", "FV_div_par_mod", + "FV_div_par_fvv", ) # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries From 9a603dab6c67be5754eae13dee421a6c08417999 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 16:12:41 +0000 Subject: [PATCH 27/39] Port `Div_par_K_Grad_par_mod` from Hermes-3 --- include/bout/difops.hxx | 4 ++ src/mesh/difops.cxx | 102 ++++++++++++++++++++++++++++ tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 6 +- tests/MMS/spatial/fci/mms.py | 1 + tests/MMS/spatial/fci/runtest | 1 + 6 files changed, 114 insertions(+), 1 deletion(-) diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index c2fdac195d..c415980d18 100644 --- a/include/bout/difops.hxx +++ b/include/bout/difops.hxx @@ -195,6 +195,10 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field2D& f, Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); +/// Version with energy flow diagnostic +Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow, + bool bndry_flux = true); + /*! * Perpendicular Laplacian operator * diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 42fa4d6ca5..5508c64ab5 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -367,6 +367,108 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc) + Div_par(kY, outloc) * Grad_par(f, outloc); } +Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, + Field3D& flow_ylow, bool bndry_flux) { + TRACE("FV::Div_par_K_Grad_par_mod"); + + ASSERT2(Kin.getLocation() == fin.getLocation()); + + Mesh* mesh = Kin.getMesh(); + Coordinates* coord = fin.getCoordinates(); + + if (Kin.hasParallelSlices() && fin.hasParallelSlices()) { + // Using parallel slices. + // Note: Y slices may use different coordinate systems + // -> Only B, dy and g_22 can be used in yup/ydown + // Others (e.g J) may not be averaged between y planes. + + const auto& K_up = Kin.yup(); + const auto& K_down = Kin.ydown(); + + const auto& f_up = fin.yup(); + const auto& f_down = fin.ydown(); + + Field3D result{zeroFrom(fin)}; + flow_ylow = zeroFrom(fin); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Upper cell edge + const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary + const BoutReal J_up = 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary + const BoutReal g_22_up = 0.5 * (coord->g_22[i] + coord->g_22.yup()[iyp]); + const BoutReal gradient_up = 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); + + const BoutReal flux_up = c_up * J_up * gradient_up / g_22_up; + + // Lower cell edge + const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary + const BoutReal J_down = 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary + const BoutReal g_22_down = 0.5 * (coord->g_22[i] + coord->g_22.ydown()[iym]); + const BoutReal gradient_down = 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); + + const BoutReal flux_down = c_down * J_down * gradient_down / g_22_down; + + result[i] = (flux_up - flux_down) / (coord->dy[i] * coord->J[i]); + } + + return result; + } + + // Calculate in field-aligned coordinates + const auto& K = toFieldAligned(Kin, "RGN_NOX"); + const auto& f = toFieldAligned(fin, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Calculate flux at upper surface + + const auto iyp = i.yp(); + const auto iym = i.ym(); + + if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) + || (i.y() != mesh->yend)) { + + BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary + BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + + BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + + BoutReal flux = c * J * gradient / g_22; + + result[i] += flux / (coord->dy[i] * coord->J[i]); + } + + // Calculate flux at lower surface + if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) + || (i.y() != mesh->ystart)) { + BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary + BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + + BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + + BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); + + BoutReal flux = c * J * gradient / g_22; + + result[i] -= flux / (coord->dy[i] * coord->J[i]); + flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; + } + } + + // Shifted to field aligned coordinates, so need to shift back + result = fromFieldAligned(result, "RGN_NOBNDRY"); + flow_ylow = fromFieldAligned(flow_ylow); + + return result; +} + + /******************************************************************************* * Delp2 * perpendicular Laplacian operator diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 99edee6ecc..76ac3035c9 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -7,6 +7,7 @@ K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) FV_div_par_fvv_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(2*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)^2/(0.01*x + 0.045)) - 1.25663706143592*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)^2/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +div_par_K_grad_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 13744c4965..17408aeeab 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -59,14 +59,18 @@ int main(int argc, char** argv) { // Add mesh geometry variables mesh->outputVars(dump); + // Dummy variable for *_mod overloads + Field3D flow_ylow; + fci_op_test("grad_par", dump, input, Grad_par(input)); fci_op_test("grad2_par2", dump, input, Grad2_par2(input)); fci_op_test("div_par", dump, input, Div_par(input)); fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); + fci_op_test("div_par_K_grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(K, input, flow_ylow)); fci_op_test("laplace_par", dump, input, Laplace_par(input)); // Finite volume methods - Field3D flow_ylow; fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); fci_op_test("FV_div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 62089c7e21..801a8d3f26 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -61,6 +61,7 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("grad2_par2_solution", FCI_grad2_par2(f)), ("div_par_solution", FCI_div_par(f)), ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), + ("div_par_K_grad_par_mod_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), ("FV_div_par_mod_solution", FCI_div_par(f * K)), ("FV_div_par_fvv_solution", FCI_div_par(f * K * K)), diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index c0f0a45132..73babc9691 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -28,6 +28,7 @@ OPERATORS = ( "grad2_par2", "div_par", "div_par_K_grad_par", + "div_par_K_grad_par_mod", "laplace_par", "FV_div_par_mod", "FV_div_par_fvv", From a5b8d5df215bf50ffa1671215a14ab4eeca13be1 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 16:59:34 +0000 Subject: [PATCH 28/39] Add MMS tests for finite volume operators --- tests/MMS/CMakeLists.txt | 1 + .../MMS/spatial/finite-volume/CMakeLists.txt | 6 + tests/MMS/spatial/finite-volume/data/BOUT.inp | 23 ++ tests/MMS/spatial/finite-volume/fv_mms.cxx | 61 +++++ tests/MMS/spatial/finite-volume/makefile | 6 + tests/MMS/spatial/finite-volume/mms.py | 62 ++++++ tests/MMS/spatial/finite-volume/runtest | 209 ++++++++++++++++++ 7 files changed, 368 insertions(+) create mode 100644 tests/MMS/spatial/finite-volume/CMakeLists.txt create mode 100644 tests/MMS/spatial/finite-volume/data/BOUT.inp create mode 100644 tests/MMS/spatial/finite-volume/fv_mms.cxx create mode 100644 tests/MMS/spatial/finite-volume/makefile create mode 100755 tests/MMS/spatial/finite-volume/mms.py create mode 100755 tests/MMS/spatial/finite-volume/runtest diff --git a/tests/MMS/CMakeLists.txt b/tests/MMS/CMakeLists.txt index 0c42da7074..cd639c9059 100644 --- a/tests/MMS/CMakeLists.txt +++ b/tests/MMS/CMakeLists.txt @@ -8,6 +8,7 @@ add_subdirectory(spatial/d2dx2) add_subdirectory(spatial/d2dz2) add_subdirectory(spatial/diffusion) add_subdirectory(spatial/fci) +add_subdirectory(spatial/finite-volume) add_subdirectory(time) add_subdirectory(time-petsc) add_subdirectory(wave-1d) diff --git a/tests/MMS/spatial/finite-volume/CMakeLists.txt b/tests/MMS/spatial/finite-volume/CMakeLists.txt new file mode 100644 index 0000000000..6d9c839a05 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/CMakeLists.txt @@ -0,0 +1,6 @@ +bout_add_mms_test(MMS-spatial-finite-volume + SOURCES fv_mms.cxx + USE_RUNTEST + USE_DATA_BOUT_INP + PROCESSORS 2 +) diff --git a/tests/MMS/spatial/finite-volume/data/BOUT.inp b/tests/MMS/spatial/finite-volume/data/BOUT.inp new file mode 100644 index 0000000000..afab4a34d5 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/data/BOUT.inp @@ -0,0 +1,23 @@ +input_field = 0.1*sin(2.0*y) + 1 +K = 0.1*cos(3.0*y) + 1 +FV_Div_par_mod_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) +FV_Div_par_fvv_solution = -0.376991118430775*(0.1*sin(2.0*y) + 1)*(0.1*cos(3.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y) +FV_Div_par_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) +FV_Div_par_K_Grad_par_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) +FV_Div_par_K_Grad_par_mod_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) + +[mesh] +MXG = 0 + +nx = 1 +ny = 128 +nz = 1 + +Ly = 10 + +dy = Ly / ny +J = 1 # Identity metric + +[mesh:ddy] +first = C2 +second = C2 diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx new file mode 100644 index 0000000000..6b45ef3259 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -0,0 +1,61 @@ +#include "bout/bout.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" +#include "bout/utils.hxx" + +#include + +#include +#include + +namespace { +auto fv_op_test(const std::string& name, Options& dump, const Field3D& input, + const Field3D& result) { + auto* mesh = input.getMesh(); + const Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), + Options::getRoot(), mesh)}; + const Field3D error{result - solution}; + + dump[fmt::format("{}_l_2", name)] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump[fmt::format("{}_l_inf", name)] = max(abs(error), true, "RGN_NOBNDRY"); + + dump[fmt::format("{}_result", name)] = result; + dump[fmt::format("{}_error", name)] = error; + dump[fmt::format("{}_input", name)] = input; + dump[fmt::format("{}_solution", name)] = solution; +} +} // namespace + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + using bout::globals::mesh; + + Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; + Field3D K{FieldFactory::get()->create3D("K", Options::getRoot(), mesh)}; + + // Communicate to calculate parallel transform. + mesh->communicate(input, K); + + Options dump; + // Add mesh geometry variables + mesh->outputVars(dump); + + // Dummy variable for *_mod overloads + Field3D flow_ylow; + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, K, K)); + fv_op_test("FV_Div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(K, input)); + fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, FV::Div_par_K_Grad_par(K, input)); + + bout::writeDefaultOutputFile(dump); + + BoutFinalise(); +} diff --git a/tests/MMS/spatial/finite-volume/makefile b/tests/MMS/spatial/finite-volume/makefile new file mode 100644 index 0000000000..88ba6c77e7 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../../.. + +SOURCEC = fci_mms.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py new file mode 100755 index 0000000000..c8b473138d --- /dev/null +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +# +# Generate manufactured solution and sources for FCI test +# + +from math import pi +import warnings + +from boututils.boutwarnings import AlwaysWarning +from boutdata.data import BoutOptionsFile +from boutdata.mms import exprToStr, y, Grad_par, Div_par, Metric +from sympy import sin, cos, Expr + +warnings.simplefilter("ignore", AlwaysWarning) + +# Length of the y domain +Ly = 10.0 + +# Identity +metric = Metric() + +# Define solution in terms of input x,y,z +f = 1 + 0.1 * sin(2 * y) +K = 1 + 0.1 * cos(3 * y) + +# Turn solution into real x and z coordinates +replace = [(y, metric.y * 2 * pi / Ly)] + +f = f.subs(replace) +K = K.subs(replace) + +# Substitute back to get input y coordinates +replace = [ (metric.y, y*Ly/(2*pi) ) ] + + +def Grad2_par2(f: Expr) -> Expr: + return Grad_par(Grad_par(f)) + + +def Div_par_K_Grad_par(f: Expr, K: Expr) -> Expr: + return (K * Grad2_par2(f)) + (Div_par(K) * Grad_par(f)) + + +############################################ +# Equations solved + +options = BoutOptionsFile("data/BOUT.inp") + +for name, expr in ( + ("input_field", f), + ("K", K), + ("FV_Div_par_solution", Div_par(f * K)), + ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, K)), + ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, K)), + ("FV_Div_par_mod_solution", Div_par(f * K)), + ("FV_Div_par_fvv_solution", Div_par(f * K * K)), +): + expr_str = exprToStr(expr.subs(replace)) + print(f"{name} = {expr_str}") + options[name] = expr_str + +options.write("data/BOUT.inp", overwrite=True) diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest new file mode 100755 index 0000000000..a836f5e735 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/runtest @@ -0,0 +1,209 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# + +import argparse +import json +import pathlib +import sys +from time import time + +from boutdata.collect import collect +from boututils.run_wrapper import build_and_log, launch_safe +from numpy import array, log, polyfit + +# Global parameters +DIRECTORY = "data" +NPROC = 2 +MTHREAD = 2 +OPERATORS = ( + "FV_Div_par", + "FV_Div_par_K_Grad_par", + "FV_Div_par_K_Grad_par_mod", + "FV_Div_par_mod", + "FV_Div_par_fvv", +) +# Resolution in y and z +NLIST = [8, 16, 32, 64] +dx = 1.0 / array(NLIST) + + +def quiet_collect(name: str) -> float: + # Index to return a plain (numpy) float rather than `BoutArray` + return collect( + name, + tind=[1, 1], + info=False, + path=DIRECTORY, + xguards=False, + yguards=False, + )[()] + + +def assert_convergence(error, dx, name, expected) -> bool: + fit = polyfit(log(dx), log(error), 1) + order = fit[0] + print(f"{name} convergence order = {order:f} (fit)", end="") + + order = log(error[-2] / error[-1]) / log(dx[-2] / dx[-1]) + print(f", {order:f} (small spacing)", end="") + + # Should be close to the expected order + success = order > expected * 0.95 + print(f"\t............ {'PASS' if success else 'FAIL'}") + + return success + + +def run_fv_operators(nz: int, name: str) -> dict[str, float]: + # Command to run + args = f"MZ={nz} mesh:ny={nz} {name}" + cmd = f"./fv_mms {args}" + print(f"Running command: {cmd}", end="") + + # Launch using MPI + start = time() + status, out = launch_safe(cmd, nproc=NPROC, mthread=MTHREAD, pipe=True) + print(f" ... done in {time() - start:.3}s") + + # Save output to log file + pathlib.Path(f"run.log.{nz}").write_text(out) + + if status: + print(f"Run failed!\nOutput was:\n{out}") + sys.exit(status) + + return { + operator: { + "l_2": quiet_collect(f"{operator}_l_2"), + "l_inf": quiet_collect(f"{operator}_l_inf"), + } + for operator in OPERATORS + } + + +def transpose( + errors: list[dict[str, dict[str, float]]], +) -> dict[str, dict[str, list[float]]]: + """Turn a list of {operator: error} into a dict of {operator: [errors]}""" + + kinds = ("l_2", "l_inf") + result = {operator: {kind: [] for kind in kinds} for operator in OPERATORS} + for error in errors: + for k, v in error.items(): + for kind in kinds: + result[k][kind].append(v[kind]) + return result + + +def check_fv_operators(name: str, case: dict) -> bool: + failures = [] + + order = case["order"] + args = case["args"] + + all_errors = [] + + for n in NLIST: + errors = run_fv_operators(n, args) + all_errors.append(errors) + + for operator in OPERATORS: + l_2 = errors[operator]["l_2"] + l_inf = errors[operator]["l_inf"] + + print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") + + final_errors = transpose(all_errors) + for operator in OPERATORS: + test_name = f"{operator} {name}" + success = assert_convergence( + final_errors[operator]["l_2"], dx, test_name, order + ) + if not success: + failures.append(test_name) + + return final_errors, failures + + +def make_plots(cases: dict[str, dict]): + try: + import matplotlib.pyplot as plt + except ImportError: + print("No matplotlib") + return + + num_operators = len(OPERATORS) + fig, axes = plt.subplots(1, num_operators, figsize=(num_operators * 4, 4)) + + for ax, operator in zip(axes, OPERATORS): + for name, case in cases.items(): + ax.loglog(dx, case[operator]["l_2"], "-", label=f"{name} $l_2$") + ax.loglog(dx, case[operator]["l_inf"], "--", label=f"{name} $l_\\inf$") + ax.legend(loc="upper left") + ax.grid() + ax.set_title(f"Error scaling for {operator}") + ax.set_xlabel(r"Mesh spacing $\delta x$") + ax.set_ylabel("Error norm") + + fig.tight_layout() + fig.savefig("fv_mms.pdf") + print("Plot saved to fv_mms.pdf") + + if args.show_plots: + plt.show() + plt.close() + + +if __name__ == "__main__": + build_and_log("Finite volume MMS test") + + parser = argparse.ArgumentParser("Error scaling test for finite volume operators") + parser.add_argument( + "--make-plots", action="store_true", help="Create plots of error scaling" + ) + parser.add_argument( + "--show-plots", + action="store_true", + help="Stop and show plots, implies --make-plots", + ) + parser.add_argument( + "--dump-errors", + type=str, + help="Output file to dump errors as JSON", + default="fv_operator_errors.json", + ) + + args = parser.parse_args() + + success = True + failures = [] + + cases = { + "default": { + "order": 2, + "args": "", + }, + } + + for name, case in cases.items(): + error2, failures_ = check_fv_operators(name, case) + case.update(error2) + failures.extend(failures_) + success &= len(failures) == 0 + + if args.dump_errors: + pathlib.Path(args.dump_errors).write_text(json.dumps(cases)) + + if args.make_plots or args.show_plots: + make_plots(cases) + + if success: + print("\nAll tests passed") + else: + print("\nSome tests failed:") + for failure in failures: + print(f"\t{failure}") + + sys.exit(0 if success else 1) From bae7618750ebad6cb635aab462cf6999bb0da210 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 14:20:14 +0000 Subject: [PATCH 29/39] Add tests for all FV slope limiters --- tests/MMS/spatial/finite-volume/data/BOUT.inp | 12 +-- tests/MMS/spatial/finite-volume/fv_mms.cxx | 51 +++++++++---- tests/MMS/spatial/finite-volume/mms.py | 18 +++-- tests/MMS/spatial/finite-volume/runtest | 73 +++++++++---------- 4 files changed, 87 insertions(+), 67 deletions(-) diff --git a/tests/MMS/spatial/finite-volume/data/BOUT.inp b/tests/MMS/spatial/finite-volume/data/BOUT.inp index afab4a34d5..029011e437 100644 --- a/tests/MMS/spatial/finite-volume/data/BOUT.inp +++ b/tests/MMS/spatial/finite-volume/data/BOUT.inp @@ -1,10 +1,10 @@ input_field = 0.1*sin(2.0*y) + 1 -K = 0.1*cos(3.0*y) + 1 -FV_Div_par_mod_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) -FV_Div_par_fvv_solution = -0.376991118430775*(0.1*sin(2.0*y) + 1)*(0.1*cos(3.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y) -FV_Div_par_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) -FV_Div_par_K_Grad_par_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) -FV_Div_par_K_Grad_par_mod_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) +FV_Div_par_mod_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_fvv_solution = -0.376991118430775*(0.1*cos(3.0*y) + 1)*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2 +FV_Div_par_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_K_Grad_par_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +FV_Div_par_K_Grad_par_mod_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +v = (0.1*cos(3.0*y) + 1)/(0.1*sin(2.0*y) + 1) [mesh] MXG = 0 diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index 6b45ef3259..09f9986e72 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -15,19 +15,20 @@ namespace { auto fv_op_test(const std::string& name, Options& dump, const Field3D& input, - const Field3D& result) { + const Field3D& result, std::string suffix = "") { auto* mesh = input.getMesh(); const Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), Options::getRoot(), mesh)}; const Field3D error{result - solution}; - dump[fmt::format("{}_l_2", name)] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); - dump[fmt::format("{}_l_inf", name)] = max(abs(error), true, "RGN_NOBNDRY"); + dump[fmt::format("{}{}_l_2", name, suffix)] = + sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump[fmt::format("{}{}_l_inf", name, suffix)] = max(abs(error), true, "RGN_NOBNDRY"); - dump[fmt::format("{}_result", name)] = result; - dump[fmt::format("{}_error", name)] = error; - dump[fmt::format("{}_input", name)] = input; - dump[fmt::format("{}_solution", name)] = solution; + dump[fmt::format("{}{}_result", name, suffix)] = result; + dump[fmt::format("{}{}_error", name, suffix)] = error; + dump[fmt::format("{}{}_input", name, suffix)] = input; + dump[fmt::format("{}{}_solution", name, suffix)] = solution; } } // namespace @@ -37,23 +38,45 @@ int main(int argc, char** argv) { using bout::globals::mesh; Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; - Field3D K{FieldFactory::get()->create3D("K", Options::getRoot(), mesh)}; + Field3D v{FieldFactory::get()->create3D("v", Options::getRoot(), mesh)}; // Communicate to calculate parallel transform. - mesh->communicate(input, K); + mesh->communicate(input, v); Options dump; // Add mesh geometry variables mesh->outputVars(dump); + dump["v"] = v; // Dummy variable for *_mod overloads Field3D flow_ylow; - fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, K, K)); - fv_op_test("FV_Div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); - fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); - fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(K, input)); - fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, FV::Div_par_K_Grad_par(K, input)); + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MC"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MC"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), "_MC"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Upwind"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Upwind"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Upwind"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Fromm"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Fromm"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Fromm"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MinMod"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MinMod"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_MinMod"); + + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(v, input)); + fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(v, input, flow_ylow)); bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py index c8b473138d..dfcfce9a09 100755 --- a/tests/MMS/spatial/finite-volume/mms.py +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -21,13 +21,15 @@ # Define solution in terms of input x,y,z f = 1 + 0.1 * sin(2 * y) -K = 1 + 0.1 * cos(3 * y) +fv = 1 + 0.1 * cos(3 * y) + # Turn solution into real x and z coordinates replace = [(y, metric.y * 2 * pi / Ly)] f = f.subs(replace) -K = K.subs(replace) +fv = fv.subs(replace) +v = fv / f # Substitute back to get input y coordinates replace = [ (metric.y, y*Ly/(2*pi) ) ] @@ -48,12 +50,12 @@ def Div_par_K_Grad_par(f: Expr, K: Expr) -> Expr: for name, expr in ( ("input_field", f), - ("K", K), - ("FV_Div_par_solution", Div_par(f * K)), - ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, K)), - ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, K)), - ("FV_Div_par_mod_solution", Div_par(f * K)), - ("FV_Div_par_fvv_solution", Div_par(f * K * K)), + ("v", v), + ("FV_Div_par_solution", Div_par(f * v)), + ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_mod_solution", Div_par(f * v)), + ("FV_Div_par_fvv_solution", Div_par(f * v * v)), ): expr_str = exprToStr(expr.subs(replace)) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest index a836f5e735..5a02be1e50 100755 --- a/tests/MMS/spatial/finite-volume/runtest +++ b/tests/MMS/spatial/finite-volume/runtest @@ -17,13 +17,27 @@ from numpy import array, log, polyfit DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ( - "FV_Div_par", - "FV_Div_par_K_Grad_par", - "FV_Div_par_K_Grad_par_mod", - "FV_Div_par_mod", - "FV_Div_par_fvv", -) +OPERATORS = { + # Slope-limiters necessarily reduce the accuracy in places + "FV_Div_par_MC": 1.5, + "FV_Div_par_mod_MC": 1.5, + "FV_Div_par_fvv_MC": 1.5, + + "FV_Div_par_Upwind": 1, + "FV_Div_par_mod_Upwind": 1, + "FV_Div_par_fvv_Upwind": 1, + + "FV_Div_par_Fromm": 1.5, + "FV_Div_par_mod_Fromm": 1.5, + "FV_Div_par_fvv_Fromm": 1.5, + + "FV_Div_par_MinMod": 1.5, + "FV_Div_par_mod_MinMod": 1.5, + "FV_Div_par_fvv_MinMod": 1.5, + + "FV_Div_par_K_Grad_par": 2, + "FV_Div_par_K_Grad_par_mod": 2, +} # Resolution in y and z NLIST = [8, 16, 32, 64] dx = 1.0 / array(NLIST) @@ -56,10 +70,9 @@ def assert_convergence(error, dx, name, expected) -> bool: return success -def run_fv_operators(nz: int, name: str) -> dict[str, float]: +def run_fv_operators(nz: int) -> dict[str, float]: # Command to run - args = f"MZ={nz} mesh:ny={nz} {name}" - cmd = f"./fv_mms {args}" + cmd = f"./fv_mms MZ={nz} mesh:ny={nz}" print(f"Running command: {cmd}", end="") # Launch using MPI @@ -97,16 +110,13 @@ def transpose( return result -def check_fv_operators(name: str, case: dict) -> bool: +def test_fv_operators() -> bool: failures = [] - order = case["order"] - args = case["args"] - all_errors = [] for n in NLIST: - errors = run_fv_operators(n, args) + errors = run_fv_operators(n) all_errors.append(errors) for operator in OPERATORS: @@ -116,13 +126,12 @@ def check_fv_operators(name: str, case: dict) -> bool: print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") final_errors = transpose(all_errors) - for operator in OPERATORS: - test_name = f"{operator} {name}" + for operator, order in OPERATORS.items(): success = assert_convergence( - final_errors[operator]["l_2"], dx, test_name, order + final_errors[operator]["l_2"], dx, operator, order ) if not success: - failures.append(test_name) + failures.append(operator) return final_errors, failures @@ -138,9 +147,8 @@ def make_plots(cases: dict[str, dict]): fig, axes = plt.subplots(1, num_operators, figsize=(num_operators * 4, 4)) for ax, operator in zip(axes, OPERATORS): - for name, case in cases.items(): - ax.loglog(dx, case[operator]["l_2"], "-", label=f"{name} $l_2$") - ax.loglog(dx, case[operator]["l_inf"], "--", label=f"{name} $l_\\inf$") + ax.loglog(dx, cases[operator]["l_2"], "-", label="$l_2$") + ax.loglog(dx, cases[operator]["l_inf"], "--", label="$l_\\inf$") ax.legend(loc="upper left") ax.grid() ax.set_title(f"Error scaling for {operator}") @@ -177,27 +185,14 @@ if __name__ == "__main__": args = parser.parse_args() - success = True - failures = [] - - cases = { - "default": { - "order": 2, - "args": "", - }, - } - - for name, case in cases.items(): - error2, failures_ = check_fv_operators(name, case) - case.update(error2) - failures.extend(failures_) - success &= len(failures) == 0 + error2, failures = test_fv_operators() + success = len(failures) == 0 if args.dump_errors: - pathlib.Path(args.dump_errors).write_text(json.dumps(cases)) + pathlib.Path(args.dump_errors).write_text(json.dumps(error2)) if args.make_plots or args.show_plots: - make_plots(cases) + make_plots(error2) if success: print("\nAll tests passed") From 39ae2bc211b569a6f03283117ae0bf6640800120 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 14:25:53 +0000 Subject: [PATCH 30/39] Port `Superbee` finite volume limiter from Hermes-3 Includes bug fix: ```diff - BoutReal gL = n.c - n.L; - BoutReal gR = n.R - n.c; + BoutReal gL = n.c - n.m; + BoutReal gR = n.p - n.c; ``` --- include/bout/fv_ops.hxx | 45 ++++++++++++++++++++++ tests/MMS/spatial/finite-volume/fv_mms.cxx | 7 ++++ tests/MMS/spatial/finite-volume/runtest | 4 ++ 3 files changed, 56 insertions(+) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 0960476911..feb6027c72 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -170,6 +170,51 @@ private: } }; +/// Superbee limiter +/// +/// This corresponds to the limiter function +/// φ(r) = max(0, min(2r, 1), min(r,2) +/// +/// The value at cell right (i.e. i + 1/2) is: +/// +/// n.R = n.c - φ(r) (n.c - (n.p + n.c)/2) +/// = n.c + φ(r) (n.p - n.c)/2 +/// +/// Four regimes: +/// a) r < 1/2 -> φ(r) = 2r +/// n.R = n.c + gL +/// b) 1/2 < r < 1 -> φ(r) = 1 +/// n.R = n.c + gR/2 +/// c) 1 < r < 2 -> φ(r) = r +/// n.R = n.c + gL/2 +/// d) 2 < r -> φ(r) = 2 +/// n.R = n.c + gR +/// +/// where the left and right gradients are: +/// gL = n.c - n.m +/// gR = n.p - n.c +/// +struct Superbee { + void operator()(Stencil1D& n) { + BoutReal gL = n.c - n.m; + BoutReal gR = n.p - n.c; + + // r = gL / gR + // Limiter is φ(r) + if (gL * gR < 0) { + // Different signs => Zero gradient + n.L = n.R = n.c; + } else { + BoutReal sign = SIGN(gL); + gL = fabs(gL); + gR = fabs(gR); + BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL)); + n.L = n.c - half_slope; + n.R = n.c + half_slope; + } + } +}; + /*! * Communicate fluxes between processors * Takes values in guard cells, and adds them to cells diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index 09f9986e72..edf4bbc16a 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -74,6 +74,13 @@ int main(int argc, char** argv) { fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), "_MinMod"); + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), + "_Superbee"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Superbee"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Superbee"); + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(v, input)); fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, Div_par_K_Grad_par_mod(v, input, flow_ylow)); diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest index 5a02be1e50..b38a6359ac 100755 --- a/tests/MMS/spatial/finite-volume/runtest +++ b/tests/MMS/spatial/finite-volume/runtest @@ -35,6 +35,10 @@ OPERATORS = { "FV_Div_par_mod_MinMod": 1.5, "FV_Div_par_fvv_MinMod": 1.5, + "FV_Div_par_Superbee": 1.5, + "FV_Div_par_mod_Superbee": 1.5, + "FV_Div_par_fvv_Superbee": 1.5, + "FV_Div_par_K_Grad_par": 2, "FV_Div_par_K_Grad_par_mod": 2, } From 4cd7909c868cd542186f40b5be489e3e0158bc0f Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 14:55:39 +0000 Subject: [PATCH 31/39] Fix clang-tidy warnings for fv_ops --- include/bout/fv_ops.hxx | 172 +++++++++---------- src/mesh/difops.cxx | 86 +++++----- src/mesh/fv_ops.cxx | 186 ++++++++++----------- tests/MMS/spatial/finite-volume/fv_mms.cxx | 1 + 4 files changed, 223 insertions(+), 222 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index feb6027c72..16f7548c20 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -5,33 +5,38 @@ #ifndef BOUT_FV_OPS_H #define BOUT_FV_OPS_H -#include - +#include "bout/assert.hxx" #include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_defines.hxx" +#include "bout/coordinates.hxx" +#include "bout/field.hxx" +#include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/globals.hxx" #include "bout/mesh.hxx" -#include "bout/output_bout_types.hxx" // NOLINT(unused-includes) +#include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner) #include "bout/region.hxx" #include "bout/utils.hxx" #include "bout/vector2d.hxx" +#include + namespace FV { /*! * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity */ -Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& x); +Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f); [[deprecated("Please use Div_a_Grad_perp instead")]] inline Field3D -Div_a_Laplace_perp(const Field3D& a, const Field3D& x) { - return Div_a_Grad_perp(a, x); +Div_a_Laplace_perp(const Field3D& a, const Field3D& f) { + return Div_a_Grad_perp(a, f); } /*! * Divergence of a parallel diffusion Div( k * Grad_par(f) ) */ -const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, - bool bndry_flux = true); +Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, bool bndry_flux = true); /*! * 4th-order derivative in Y, using derivatives @@ -53,7 +58,7 @@ const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, * * No fluxes through domain boundaries */ -const Field3D D4DY4(const Field3D& d, const Field3D& f); +Field3D D4DY4(const Field3D& d, const Field3D& f); /*! * 4th-order dissipation term @@ -71,18 +76,24 @@ const Field3D D4DY4(const Field3D& d, const Field3D& f); * f_2 | f_1 | f_0 | * f_b */ -const Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); +Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); /*! * Stencil used for Finite Volume calculations * which includes cell face values L and R */ struct Stencil1D { - // Cell centre values - BoutReal c, m, p, mm, pp; - - // Left and right cell face values - BoutReal L, R; + /// Cell centre values + BoutReal c; + BoutReal m; + BoutReal p; + BoutReal mm = BoutNaN; + BoutReal pp = BoutNaN; + + /// Left cell face value + BoutReal L = BoutNaN; + /// Right cell face value + BoutReal R = BoutNaN; }; /*! @@ -97,8 +108,8 @@ struct Upwind { */ struct Fromm { void operator()(Stencil1D& n) { - n.L = n.c - 0.25 * (n.p - n.m); - n.R = n.c + 0.25 * (n.p - n.m); + n.L = n.c - (0.25 * (n.p - n.m)); + n.R = n.c + (0.25 * (n.p - n.m)); } }; @@ -114,9 +125,9 @@ struct MinMod { void operator()(Stencil1D& n) { // Choose the gradient within the cell // as the minimum (smoothest) solution - BoutReal slope = _minmod(n.p - n.c, n.c - n.m); - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; + const BoutReal slope = _minmod(n.p - n.c, n.c - n.m); + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); } private: @@ -127,7 +138,7 @@ private: * returns zero, otherwise chooses the value * with the minimum magnitude. */ - BoutReal _minmod(BoutReal a, BoutReal b) { + static BoutReal _minmod(BoutReal a, BoutReal b) { if (a * b <= 0.0) { return 0.0; } @@ -149,17 +160,17 @@ private: */ struct MC { void operator()(Stencil1D& n) { - BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference - 0.5 * (n.p - n.m), // Central difference - 2. * (n.c - n.m)); // 2*left difference - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; + const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference + 0.5 * (n.p - n.m), // Central difference + 2. * (n.c - n.m)); // 2*left difference + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); } private: // Return zero if any signs are different // otherwise return the value with the minimum magnitude - BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { + static BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { // if any of the signs are different, return zero gradient if ((a * b <= 0.0) || (a * c <= 0.0)) { return 0.0; @@ -196,8 +207,8 @@ private: /// struct Superbee { void operator()(Stencil1D& n) { - BoutReal gL = n.c - n.m; - BoutReal gR = n.p - n.c; + const BoutReal gL = n.c - n.m; + const BoutReal gR = n.p - n.c; // r = gL / gR // Limiter is φ(r) @@ -205,10 +216,11 @@ struct Superbee { // Different signs => Zero gradient n.L = n.R = n.c; } else { - BoutReal sign = SIGN(gL); - gL = fabs(gL); - gR = fabs(gR); - BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL)); + const BoutReal sign = SIGN(gL); + const BoutReal abs_gL = fabs(gL); + const BoutReal abs_gR = fabs(gR); + const BoutReal half_slope = + sign * BOUTMAX(BOUTMIN(abs_gL, 0.5 * abs_gR), BOUTMIN(abs_gR, 0.5 * abs_gL)); n.L = n.c - half_slope; n.R = n.c + half_slope; } @@ -238,13 +250,13 @@ void communicateFluxes(Field3D& f); /// /// NB: Uses to/from FieldAligned coordinates template -const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, bool fixflux = true) { +Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, + bool fixflux = true) { ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - Mesh* mesh = f_in.getMesh(); + Mesh const* mesh = f_in.getMesh(); CellEdges cellboundary; @@ -264,29 +276,17 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, Field3D result{zeroFrom(f)}; - // Only need one guard cell, so no need to communicate fluxes - // Instead calculate in guard cells to preserve fluxes - int ys = mesh->ystart - 1; - int ye = mesh->yend + 1; - for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); - if (!mesh->firstY(i) || mesh->periodicY(i)) { - // Calculate in guard cell to get fluxes consistent between processors - ys = mesh->ystart - 1; - } else { - // Don't include the boundary cell. Note that this implies special - // handling of boundaries later - ys = mesh->ystart; - } - - if (!mesh->lastY(i) || mesh->periodicY(i)) { - // Calculate in guard cells - ye = mesh->yend + 1; - } else { - // Not in boundary cells - ye = mesh->yend; - } + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; for (int j = ys; j <= ye; j++) { // Pre-calculate factors which multiply fluxes @@ -295,16 +295,16 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_rp = + const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); // For left cell boundaries common_factor = (coord->J(i, j) + coord->J(i, j - 1)) / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_lm = + const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); #endif for (int k = 0; k < mesh->LocalNz; k++) { @@ -347,23 +347,23 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, // Calculate velocity at right boundary (y+1/2) BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - BoutReal flux; + BoutReal flux = NAN; - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { // Last point in domain - BoutReal bndryval = 0.5 * (s.c + s.p); + const BoutReal bndryval = 0.5 * (s.c + s.p); if (fixflux) { // Use mid-point to be consistent with boundary conditions flux = bndryval * vpar; } else { // Add flux due to difference in boundary values - flux = s.R * vpar + wave_speed(i, j, k) * (s.R - bndryval); + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); } } else { // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); if (vpar > amax) { // Supersonic flow out of this cell @@ -385,20 +385,20 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { // First point in domain - BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal bndryval = 0.5 * (s.c + s.m); if (fixflux) { // Use mid-point to be consistent with boundary conditions flux = bndryval * vpar; } else { // Add flux due to difference in boundary values - flux = s.L * vpar - wave_speed(i, j, k) * (s.L - bndryval); + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); } } else { // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); if (vpar < -amax) { // Supersonic out of this cell @@ -432,11 +432,11 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, * */ template -const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { +Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { ASSERT1(n_in.getLocation() == v.getLocation()); ASSERT1_FIELDS_COMPATIBLE(n_in, v.x); - Mesh* mesh = n_in.getMesh(); + const Mesh* mesh = n_in.getMesh(); CellEdges cellboundary; @@ -455,10 +455,10 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Calculate velocities - BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); - BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); - BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); - BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); + const BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); + const BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); + const BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); + const BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); // X direction Stencil1D s; @@ -473,7 +473,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { if ((i.x() == mesh->xend) && (mesh->lastX())) { // At right boundary in X if (bndry_flux) { - BoutReal flux; + BoutReal flux = NAN; if (vR > 0.0) { // Flux to boundary flux = vR * s.R; @@ -488,7 +488,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { // Not at a boundary if (vR > 0.0) { // Flux out into next cell - BoutReal flux = vR * s.R; + const BoutReal flux = vR * s.R; result[i] += flux / (coord->dx[i] * coord->J[i]); result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); } @@ -500,7 +500,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { // At left boundary in X if (bndry_flux) { - BoutReal flux; + BoutReal flux = NAN; if (vL < 0.0) { // Flux to boundary flux = vL * s.L; @@ -514,7 +514,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { } else { // Not at a boundary if (vL < 0.0) { - BoutReal flux = vL * s.L; + const BoutReal flux = vL * s.L; result[i] -= flux / (coord->dx[i] * coord->J[i]); result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); } @@ -531,12 +531,12 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { cellboundary(s); if (vU > 0.0) { - BoutReal flux = vU * s.R; + const BoutReal flux = vU * s.R; result[i] += flux / (coord->J[i] * coord->dz[i]); result[i.zp()] -= flux / (coord->J[i.zp()] * coord->dz[i.zp()]); } if (vD < 0.0) { - BoutReal flux = vD * s.L; + const BoutReal flux = vD * s.L; result[i] -= flux / (coord->J[i] * coord->dz[i]); result[i.zm()] += flux / (coord->J[i.zm()] * coord->dz[i.zm()]); } @@ -556,13 +556,13 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Y velocities on y boundaries - BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); - BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); + const BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); + const BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); // n (advected quantity) on y boundaries // Note: Use unshifted n_in variable - BoutReal nU = 0.5 * (n[i] + n[i.yp()]); - BoutReal nD = 0.5 * (n[i] + n[i.ym()]); + const BoutReal nU = 0.5 * (n[i] + n[i.yp()]); + const BoutReal nD = 0.5 * (n[i] + n[i.ym()]); yresult[i] = (nU * vU - nD * vD) / (coord->J[i] * coord->dy[i]); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 5508c64ab5..1e0651abca 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -25,20 +25,20 @@ #include "bout/build_defines.hxx" -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include // Delp2 uses same coefficients as inversion code - -#include -#include +#include "bout/assert.hxx" +#include "bout/derivs.hxx" +#include "bout/difops.hxx" +#include "bout/fft.hxx" +#include "bout/field2d.hxx" +#include "bout/globals.hxx" +#include "bout/interpolation.hxx" +#include "bout/invert_laplace.hxx" // Delp2 uses same coefficients as inversion code +#include "bout/msg_stack.hxx" +#include "bout/region.hxx" +#include "bout/solver.hxx" +#include "bout/unused.hxx" +#include "bout/utils.hxx" +#include "bout/vecops.hxx" #include @@ -367,14 +367,14 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc) + Div_par(kY, outloc) * Grad_par(f, outloc); } -Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, - Field3D& flow_ylow, bool bndry_flux) { +Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, Field3D& flow_ylow, + bool bndry_flux) { TRACE("FV::Div_par_K_Grad_par_mod"); ASSERT2(Kin.getLocation() == fin.getLocation()); - Mesh* mesh = Kin.getMesh(); - Coordinates* coord = fin.getCoordinates(); + const Mesh* mesh = Kin.getMesh(); + const Coordinates* coord = fin.getCoordinates(); if (Kin.hasParallelSlices() && fin.hasParallelSlices()) { // Using parallel slices. @@ -396,18 +396,22 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, const auto iym = i.ym(); // Upper cell edge - const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary - const BoutReal J_up = 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary + const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary + const BoutReal J_up = + 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary const BoutReal g_22_up = 0.5 * (coord->g_22[i] + coord->g_22.yup()[iyp]); - const BoutReal gradient_up = 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); + const BoutReal gradient_up = + 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); const BoutReal flux_up = c_up * J_up * gradient_up / g_22_up; // Lower cell edge - const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary - const BoutReal J_down = 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary + const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary + const BoutReal J_down = + 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary const BoutReal g_22_down = 0.5 * (coord->g_22[i] + coord->g_22.ydown()[iym]); - const BoutReal gradient_down = 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); + const BoutReal gradient_down = + 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); const BoutReal flux_down = c_down * J_down * gradient_down / g_22_down; @@ -426,35 +430,32 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Calculate flux at upper surface - + const auto ix = i.x(); + const auto iy = i.y(); const auto iyp = i.yp(); const auto iym = i.ym(); - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) - || (i.y() != mesh->yend)) { - - BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const bool is_periodic_y = mesh->periodicY(ix); - BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + if (bndry_flux || is_periodic_y || !mesh->lastY(ix) || (iy != mesh->yend)) { + const BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] += flux / (coord->dy[i] * coord->J[i]); } // Calculate flux at lower surface - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) - || (i.y() != mesh->ystart)) { - BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + if (bndry_flux || is_periodic_y || !mesh->firstY(ix) || (iy != mesh->ystart)) { + const BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - - BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); - - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] -= flux / (coord->dy[i] * coord->J[i]); flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; @@ -468,7 +469,6 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, return result; } - /******************************************************************************* * Delp2 * perpendicular Laplacian operator diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index fe5422b4d1..71e51561b0 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -1,8 +1,16 @@ -#include -#include -#include -#include -#include +#include "bout/fv_ops.hxx" + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_config.hxx" +#include "bout/coordinates.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" +#include "bout/globals.hxx" +#include "bout/msg_stack.hxx" +#include "bout/region.hxx" +#include "bout/utils.hxx" namespace { template @@ -34,28 +42,19 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { // Flux in x - int xs = mesh->xstart - 1; - int xe = mesh->xend; - - /* - if(mesh->firstX()) - xs += 1; - */ - /* - if(mesh->lastX()) - xe -= 1; - */ + const int xs = mesh->xstart - 1; + const int xe = mesh->xend; for (int i = xs; i <= xe; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { // Calculate flux from i to i+1 - BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) - * (coord->J(i, j, k) * coord->g11(i, j, k) - + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + const BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); @@ -179,15 +178,14 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { return result; } -const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, - bool bndry_flux) { +Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, bool bndry_flux) { TRACE("FV::Div_par_K_Grad_par"); ASSERT2(Kin.getLocation() == fin.getLocation()); - Mesh* mesh = Kin.getMesh(); + const Mesh* mesh = Kin.getMesh(); - bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); + const bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); const auto& K = use_parallel_slices ? Kin : toFieldAligned(Kin, "RGN_NOX"); const auto& f = use_parallel_slices ? fin : toFieldAligned(fin, "RGN_NOX"); @@ -211,13 +209,13 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) || (i.y() != mesh->yend)) { - BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); - BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + const BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] += flux / (coord->dy[i] * coord->J[i]); } @@ -225,14 +223,15 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, // Calculate flux at lower surface if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) || (i.y() != mesh->ystart)) { - BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - BoutReal gradient = 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); + const BoutReal gradient = + 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] -= flux / (coord->dy[i] * coord->J[i]); } @@ -246,10 +245,10 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, return result; } -const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { +Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { ASSERT1_FIELDS_COMPATIBLE(d_in, f_in); - Mesh* mesh = d_in.getMesh(); + const Mesh* mesh = d_in.getMesh(); Coordinates* coord = f_in.getCoordinates(); @@ -265,9 +264,9 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int i = mesh->xstart; i <= mesh->xend; i++) { // Check for boundaries - bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool yperiodic = mesh->periodicY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); // Always calculate fluxes at upper Y cell boundary const int ystart = @@ -283,15 +282,15 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int j = ystart; j <= yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { - BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); + const BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); // 3rd derivative at upper boundary - BoutReal d3fdy3 = + const BoutReal d3fdy3 = (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)) / dy3; - BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) - * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; + const BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) + * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; result(i, j, k) += flux / (coord->J(i, j, k) * coord->dy(i, j, k)); result(i, j + 1, k) -= flux / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); @@ -303,8 +302,8 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } -const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { - Mesh* mesh = f_in.getMesh(); +Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { + const Mesh* mesh = f_in.getMesh(); // Convert to field aligned coordinates const bool is_unaligned = (f_in.getDirectionY() == YDirectionType::Standard); @@ -315,10 +314,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { Coordinates* coord = f_in.getCoordinates(); for (int i = mesh->xstart; i <= mesh->xend; i++) { - bool yperiodic = mesh->periodicY(i); + const bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); for (int j = mesh->ystart; j <= mesh->yend; j++) { @@ -343,8 +342,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { // Not on domain boundary // 3rd derivative at right cell boundary - const BoutReal d3fdx3 = - (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)); + const BoutReal d3fdx3 = (f(i, j + 2, k) - (3. * f(i, j + 1, k)) + + (3. * f(i, j, k)) - f(i, j - 1, k)); result(i, j, k) += d3fdx3 * factor_rc; result(i, j + 1, k) -= d3fdx3 * factor_rp; @@ -365,10 +364,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); const BoutReal d3fdx3 = - -((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k)) // Boundary value f_b - - 6. * f(i, j, k) // f_0 - + 4. * f(i, j - 1, k) // f_1 - - (6. / 5) * f(i, j - 2, k) // f_2 + -(((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k))) // Boundary value f_b + - (6. * f(i, j, k)) // f_0 + + (4. * f(i, j - 1, k)) // f_1 + - ((6. / 5) * f(i, j - 2, k)) // f_2 ); result(i, j, k) += d3fdx3 * factor_rc; @@ -394,8 +393,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); // Not on a domain boundary - const BoutReal d3fdx3 = - (f(i, j + 1, k) - 3. * f(i, j, k) + 3. * f(i, j - 1, k) - f(i, j - 2, k)); + const BoutReal d3fdx3 = (f(i, j + 1, k) - (3. * f(i, j, k)) + + (3. * f(i, j - 1, k)) - f(i, j - 2, k)); result(i, j, k) -= d3fdx3 * factor_lc; result(i, j - 1, k) += d3fdx3 * factor_lm; @@ -412,10 +411,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { const BoutReal factor_lm = common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); const BoutReal d3fdx3 = - -(-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k)) // Boundary value f_b - + 6. * f(i, j, k) // f_0 - - 4. * f(i, j + 1, k) // f_1 - + (6. / 5) * f(i, j + 2, k) // f_2 + -((-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k))) // Boundary value f_b + + (6. * f(i, j, k)) // f_0 + - (4. * f(i, j + 1, k)) // f_1 + + ((6. / 5) * f(i, j + 2, k)) // f_2 ); result(i, j, k) -= d3fdx3 * factor_lc; @@ -438,8 +437,9 @@ void communicateFluxes(Field3D& f) { throw BoutException("communicateFluxes: Sorry!"); } - int size = mesh->LocalNy * mesh->LocalNz; - comm_handle xin, xout; + const int size = mesh->LocalNy * mesh->LocalNz; + comm_handle xin = nullptr; + comm_handle xout = nullptr; // Cache results to silence spurious compiler warning about xin, // xout possibly being uninitialised when used const bool not_first = mesh->periodicX || !mesh->firstX(); @@ -498,45 +498,45 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { // o --- gD --- o // Coordinates* coords = a.getCoordinates(outloc); - Mesh* mesh = f.getMesh(); + const Mesh* mesh = f.getMesh(); for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { // wrap k-index around as Z is (currently) periodic. - int kp = (k + 1) % (mesh->LocalNz); - int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); + const int kp = (k + 1) % (mesh->LocalNz); + const int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); // Calculate gradients on cell faces -- assumes constant grid spacing - BoutReal gR = - (coords->g11(i, j, k) + coords->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coords->dx(i + 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) - * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4. * coords->dz(i, j, k)); - - BoutReal gL = - (coords->g11(i - 1, j, k) + coords->g11(i, j, k)) - * (f(i, j, k) - f(i - 1, j, k)) - / (coords->dx(i - 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) - * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4 * coords->dz(i, j, k)); - - BoutReal gD = - coords->g13(i, j, k) - * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k); - - BoutReal gU = - coords->g13(i, j, k) - * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k); + const BoutReal gR = + ((coords->g11(i, j, k) + coords->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coords->dx(i + 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) + * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4. * coords->dz(i, j, k))); + + const BoutReal gL = + ((coords->g11(i - 1, j, k) + coords->g11(i, j, k)) + * (f(i, j, k) - f(i - 1, j, k)) + / (coords->dx(i - 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) + * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4 * coords->dz(i, j, k))); + + const BoutReal gD = + (coords->g13(i, j, k) + * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k)); + + const BoutReal gU = + (coords->g13(i, j, k) + * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k)); // Flow right BoutReal flux = gR * 0.25 * (coords->J(i + 1, j, k) + coords->J(i, j, k)) diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index edf4bbc16a..19f3f14610 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -1,4 +1,5 @@ #include "bout/bout.hxx" +#include "bout/difops.hxx" #include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" From 97499913736b01d1609ee3fc199ab677a13daaad Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 15:34:52 +0000 Subject: [PATCH 32/39] Apply ruff format --- tests/MMS/spatial/finite-volume/mms.py | 2 +- tests/MMS/spatial/finite-volume/runtest | 9 +-------- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py index dfcfce9a09..a95ecc1328 100755 --- a/tests/MMS/spatial/finite-volume/mms.py +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -32,7 +32,7 @@ v = fv / f # Substitute back to get input y coordinates -replace = [ (metric.y, y*Ly/(2*pi) ) ] +replace = [(metric.y, y * Ly / (2 * pi))] def Grad2_par2(f: Expr) -> Expr: diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest index b38a6359ac..bcd4672545 100755 --- a/tests/MMS/spatial/finite-volume/runtest +++ b/tests/MMS/spatial/finite-volume/runtest @@ -22,23 +22,18 @@ OPERATORS = { "FV_Div_par_MC": 1.5, "FV_Div_par_mod_MC": 1.5, "FV_Div_par_fvv_MC": 1.5, - "FV_Div_par_Upwind": 1, "FV_Div_par_mod_Upwind": 1, "FV_Div_par_fvv_Upwind": 1, - "FV_Div_par_Fromm": 1.5, "FV_Div_par_mod_Fromm": 1.5, "FV_Div_par_fvv_Fromm": 1.5, - "FV_Div_par_MinMod": 1.5, "FV_Div_par_mod_MinMod": 1.5, "FV_Div_par_fvv_MinMod": 1.5, - "FV_Div_par_Superbee": 1.5, "FV_Div_par_mod_Superbee": 1.5, "FV_Div_par_fvv_Superbee": 1.5, - "FV_Div_par_K_Grad_par": 2, "FV_Div_par_K_Grad_par_mod": 2, } @@ -131,9 +126,7 @@ def test_fv_operators() -> bool: final_errors = transpose(all_errors) for operator, order in OPERATORS.items(): - success = assert_convergence( - final_errors[operator]["l_2"], dx, operator, order - ) + success = assert_convergence(final_errors[operator]["l_2"], dx, operator, order) if not success: failures.append(operator) From 07ba3b8d99aa03b14ed7d500f72cc63be2688d7c Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 16:37:31 +0000 Subject: [PATCH 33/39] tests: Communicate some coordinates fields for 3D metrics --- tests/MMS/spatial/finite-volume/fv_mms.cxx | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index 19f3f14610..a2a977aeef 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -42,6 +42,15 @@ int main(int argc, char** argv) { Field3D v{FieldFactory::get()->create3D("v", Options::getRoot(), mesh)}; // Communicate to calculate parallel transform. + if constexpr (bout::build::use_metric_3d) { + // Div_par operators require B parallel slices: + // Coordinates::geometry doesn't ensure this (yet) + auto& Bxy = mesh->getCoordinates()->Bxy; + auto& J = mesh->getCoordinates()->J; + auto& dy = mesh->getCoordinates()->dy; + auto& g_22 = mesh->getCoordinates()->g_22; + mesh->communicate(Bxy, J, dy, g_22); + } mesh->communicate(input, v); Options dump; From 7857a6bd46a654b338ba691423503e64df21715a Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 17 Nov 2025 17:21:49 +0000 Subject: [PATCH 34/39] Include all headers --- tests/MMS/spatial/finite-volume/fv_mms.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index a2a977aeef..bb999bcc66 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -1,4 +1,5 @@ #include "bout/bout.hxx" +#include "bout/build_config.hxx" #include "bout/difops.hxx" #include "bout/field.hxx" #include "bout/field3d.hxx" From 706f42892a4c23d7ffd487534d76129c6d4c240b Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 3 Dec 2025 14:19:48 +0000 Subject: [PATCH 35/39] Div_par_fvv: FCI change numerical dissipation Damp gradients of velocity rather than gradients of momentum --- include/bout/fv_ops.hxx | 28 ++++++++++------------------ 1 file changed, 10 insertions(+), 18 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 16f7548c20..678e7499c7 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -858,24 +858,16 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const BoutReal amax = BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); - const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) - - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]); - - // Penalty terms. This implementation is very dissipative. - BoutReal penalty = - (amax * (f_in[i] * v_in[i] - f_up[iyp] * v_up[iyp]) / (B[i] + B_up[iyp])) - + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym]) - / (B[i] + B_down[iym])); - - if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) { - if (term * penalty > 0) { - penalty = term; - } else { - penalty = -term; - } - } - - result[i] = B[i] * (term + penalty) / (2 * dy[i] * sqrt(g_22[i])); + result[i] = + B[i] + * ((f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) + - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]) + // Penalty terms. This implementation is very dissipative. + // Note: This version adds a viscosity that damps gradients of velocity + + amax * (f_in[i] + f_up[iyp]) * (v_in[i] - v_up[iyp]) / (B[i] + B_up[iyp]) + + amax * (f_in[i] + f_down[iym]) * (v_in[i] - v_down[iym]) + / (B[i] + B_down[iym])) + / (2 * dy[i] * sqrt(g_22[i])); #if CHECK > 0 if (!std::isfinite(result[i])) { From 21cb28a93eb8cac2e203df9464c2784445ab7e4f Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 3 Dec 2025 14:20:12 +0000 Subject: [PATCH 36/39] Set expected order of `FV::Div_par_fvv` to 1 --- tests/MMS/spatial/fci/runtest | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 73babc9691..1e6d570c96 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -23,16 +23,16 @@ from scipy.interpolate import RectBivariateSpline as RBS DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ( - "grad_par", - "grad2_par2", - "div_par", - "div_par_K_grad_par", - "div_par_K_grad_par_mod", - "laplace_par", - "FV_div_par_mod", - "FV_div_par_fvv", -) +OPERATORS = { + "grad_par": 2, + "grad2_par2": 2, + "div_par": 2, + "div_par_K_grad_par": 2, + "div_par_K_grad_par_mod": 2, + "laplace_par": 2, + "FV_div_par_mod": 2, + "FV_div_par_fvv": 1, +} # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries NX = 4 @@ -169,10 +169,11 @@ def check_fci_operators(name: str, case: dict) -> bool: print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") final_errors = transpose(all_errors) - for operator in OPERATORS: + for operator, operator_order in OPERATORS.items(): test_name = f"{operator} {name}" + expected_order = min(order, operator_order) success = assert_convergence( - final_errors[operator]["l_2"], dx, test_name, order + final_errors[operator]["l_2"], dx, test_name, expected_order ) if not success: failures.append(test_name) From 972a19db7dbf622e99f608136e77da4bf35ba32a Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 4 Dec 2025 09:43:50 +0000 Subject: [PATCH 37/39] Communicate some more coordinate quantities for FCI MMS test --- tests/MMS/spatial/fci/fci_mms.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 17408aeeab..3a9b72070e 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -51,7 +51,10 @@ int main(int argc, char** argv) { // Div_par operators require B parallel slices: // Coordinates::geometry doesn't ensure this (yet) auto& Bxy = mesh->getCoordinates()->Bxy; - mesh->communicate(Bxy); + auto& J = mesh->getCoordinates()->J; + auto& g_22 = mesh->getCoordinates()->g_22; + auto& dy = mesh->getCoordinates()->dy; + mesh->communicate(Bxy, J, g_22, dy); } mesh->communicate(input, K); From b6c17057a4422c954965ddc4dcc271d5eee9dee2 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 4 Dec 2025 09:44:10 +0000 Subject: [PATCH 38/39] Guard against non-thread safe PETSc call --- src/mesh/interpolation/hermite_spline_xz.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 27a4f1d614..c58c50ddb5 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -25,6 +25,8 @@ #include "bout/globals.hxx" #include "bout/index_derivs_interface.hxx" #include "bout/interpolation_xz.hxx" +#include "bout/openmpwrap.hxx" +#include "bout/region.hxx" #include @@ -296,6 +298,7 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z k_corner(x, y, z) - 1 + k); vals[k] = newWeights[j * 4 + k][i]; } + BOUT_OMP(critical) MatSetValues(petscWeights, 1, idxn, 4, idxm, vals, INSERT_VALUES); } #endif From e0957087afec7921bf24859f7fdaef65f52a137b Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 4 Dec 2025 09:45:39 +0000 Subject: [PATCH 39/39] Remove some unused headers --- src/mesh/difops.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 1e0651abca..8cc1a7c28e 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -28,7 +28,6 @@ #include "bout/assert.hxx" #include "bout/derivs.hxx" #include "bout/difops.hxx" -#include "bout/fft.hxx" #include "bout/field2d.hxx" #include "bout/globals.hxx" #include "bout/interpolation.hxx" @@ -38,7 +37,6 @@ #include "bout/solver.hxx" #include "bout/unused.hxx" #include "bout/utils.hxx" -#include "bout/vecops.hxx" #include