diff --git a/.github/workflows/convergence.yml b/.github/workflows/convergence.yml new file mode 100644 index 0000000000..c652125ec4 --- /dev/null +++ b/.github/workflows/convergence.yml @@ -0,0 +1,36 @@ +name: Convergence + +on: + push: + branches: [master] + pull_request: + workflow_dispatch: + +env: + OMPI_MCA_rmaps_base_oversubscribe: 1 + +jobs: + convergence: + name: "Convergence" + runs-on: ubuntu-latest + timeout-minutes: 240 + + steps: + - uses: actions/checkout@v4 + + - name: Setup Ubuntu + run: | + sudo apt update -y + sudo apt install -y cmake gcc g++ python3 python3-dev \ + openmpi-bin libopenmpi-dev + + - name: Setup Python + uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Build MFC + run: ./mfc.sh build -j 4 + + - name: Run convergence tests + run: ./mfc.sh test --only Convergence --no-build -j 1 diff --git a/examples/1D_advection_convergence/case.py b/examples/1D_advection_convergence/case.py new file mode 100644 index 0000000000..f8fbef41f1 --- /dev/null +++ b/examples/1D_advection_convergence/case.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python3 +""" +1D periodic advection convergence case. + +Two identical fluids (same gamma, same density = 1) with a sine-wave volume +fraction. Since both EOS are identical, alpha_1 advects passively at u = 1 +with no acoustic coupling. After exactly one period (T = L/u = 1), the +exact solution equals the IC, so L2(q_cons_vf1(T) - q_cons_vf1(0)) is +purely the scheme's accumulated spatial truncation error. +""" + +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="1D advection convergence case") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=64, help="Grid points (default: 64)") +parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, or 5") +parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO") +parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4)") +parser.add_argument("--mp-weno", action="store_true", help="Enable MP-WENO limiter") +parser.add_argument("--muscl-lim", type=int, default=1, help="MUSCL limiter: 1=minmod 2=MC 3=VanAlbada 4=VanLeer 5=Superbee") +args = parser.parse_args() + +gamma = 1.4 +N = args.N +m = N - 1 +L = 1.0 +dx = L / N + +# Max wave speed: acoustic speed + convective speed +# c_sound = sqrt(gamma * p / rho) = sqrt(gamma) ≈ 1.183 (for p=1, rho=1) +c_max = math.sqrt(gamma) + 1.0 +dt = args.cfl * dx / c_max +T_end = 1.0 # exactly one period: u=1, L=1 +Nt = max(4, math.ceil(T_end / dt)) +dt = T_end / Nt # snap to land exactly on T_end + +if args.muscl: + scheme_params = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": args.muscl_lim, + } +else: + scheme_params = { + "recon_type": 1, + "weno_order": args.order, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F" if args.order == 1 else "T", + "null_weights": "F", + "mp_weno": "T" if args.mp_weno else "F", + } + +print( + json.dumps( + { + "run_time_info": "F", + "x_domain%beg": 0.0, + "x_domain%end": L, + "m": m, + "n": 0, + "p": 0, + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 2, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%length_x": L, + "patch_icpp(1)%vel(1)": 1.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": "0.5 + 0.2 * sin(2.0 * pi * x / lx)", + "patch_icpp(1)%alpha_rho(2)": "0.5 - 0.2 * sin(2.0 * pi * x / lx)", + "patch_icpp(1)%alpha(1)": "0.5 + 0.2 * sin(2.0 * pi * x / lx)", + "patch_icpp(1)%alpha(2)": "0.5 - 0.2 * sin(2.0 * pi * x / lx)", + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%pi_inf": 0.0, + "fluid_pp(2)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(2)%pi_inf": 0.0, + **scheme_params, + } + ) +) diff --git a/examples/1D_euler_convergence/case.py b/examples/1D_euler_convergence/case.py new file mode 100644 index 0000000000..d67091328e --- /dev/null +++ b/examples/1D_euler_convergence/case.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python3 +""" +1D single-fluid Euler convergence case. + +Single fluid with a density sine wave: rho = 1 + 0.2*sin(2*pi*x). +Constant velocity u=1 and pressure p=1. For this IC, the Euler equations +reduce to pure advection of all variables at speed u=1. After exactly one +period (T = L/u = 1), the exact solution equals the IC, so +L2(rho(T) - rho(0)) measures the accumulated scheme spatial truncation error. + +No non-conservative alpha equation — clean benchmark for WENO/MUSCL rates. +""" + +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="1D Euler convergence case") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=64, help="Grid points (default: 64)") +parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, 5, or 7") +parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO") +parser.add_argument("--teno", action="store_true", help="Use TENO instead of WENO") +parser.add_argument("--teno-ct", type=float, default=1e-6, help="TENO CT threshold (default: 1e-6)") +parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4)") +parser.add_argument("--no-mapped", action="store_true", help="Disable mapped WENO") +parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)") +parser.add_argument("--time-stepper", type=int, default=3, help="Time stepper: 1=Euler 2=RK2 3=RK3 (default: 3)") +parser.add_argument("--t-end", type=float, default=None, help="Override total simulation time (default: 1.0 = one period)") +args = parser.parse_args() + +gamma = 1.4 +N = args.N +m = N - 1 +L = 1.0 +dx = L / N + +# c_sound = sqrt(gamma * p / rho) = sqrt(gamma) for p=1, rho=1 +c_max = math.sqrt(gamma) + 1.0 # acoustic + convective +dt = args.cfl * dx / c_max +T_end = args.t_end if args.t_end is not None else 1.0 +Nt = max(1, math.ceil(T_end / dt)) +dt = T_end / Nt + +if args.muscl: + scheme_params = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": args.muscl_lim, + } +else: + scheme_params = { + "recon_type": 1, + "weno_order": args.order, + "weno_eps": 1.0e-40, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F" if (args.order == 1 or args.no_mapped or args.teno) else "T", + "null_weights": "F", + "mp_weno": "F", + "teno": "T" if args.teno else "F", + **({"teno_CT": args.teno_ct} if args.teno else {}), + } + +print( + json.dumps( + { + "run_time_info": "F", + "x_domain%beg": 0.0, + "x_domain%end": L, + "m": m, + "n": 0, + "p": 0, + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": args.time_stepper, + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%length_x": L, + "patch_icpp(1)%vel(1)": 1.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": "1.0 + 0.2 * sin(2.0 * pi * x / lx)", + "patch_icpp(1)%alpha(1)": 1.0, + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%pi_inf": 0.0, + **scheme_params, + } + ) +) diff --git a/examples/1D_sod_convergence/case.py b/examples/1D_sod_convergence/case.py new file mode 100644 index 0000000000..cf726519e9 --- /dev/null +++ b/examples/1D_sod_convergence/case.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python3 +""" +1D Sod shock tube convergence case. + +Standard Sod problem: rho_L=1, u_L=0, p_L=1; rho_R=0.125, u_R=0, p_R=0.1. +Discontinuity at x=0.5, gamma=1.4, T_end=0.2 (shock, contact, rarefaction). +Used for L1 self-convergence study; outflow BCs (-3) at both ends. +""" + +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="1D Sod shock tube convergence case") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=128, help="Grid points (default: 128)") +parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, 5, or 7 (default: 5)") +parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO") +parser.add_argument("--muscl-lim", type=int, default=1, help="MUSCL limiter: 1=minmod 2=MC 4=VanLeer 5=SUPERBEE (default: 1)") +parser.add_argument("--teno", action="store_true", help="Use TENO instead of WENO") +parser.add_argument("--teno-ct", type=float, default=1e-6, help="TENO CT threshold (default: 1e-6)") +parser.add_argument("--cfl", type=float, default=0.3, help="CFL number (default: 0.3)") +args = parser.parse_args() + +gamma = 1.4 +N = args.N +m = N - 1 +L = 1.0 +dx = L / N + +c_max = math.sqrt(gamma) + 1.0 # conservative: sound speed + max velocity +dt = args.cfl * dx / c_max +T_end = 0.2 +Nt = max(4, math.ceil(T_end / dt)) +dt = T_end / Nt + +if args.muscl: + scheme_params = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": args.muscl_lim, + } +else: + scheme_params = { + "recon_type": 1, + "weno_order": args.order, + "weno_eps": 1.0e-40, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F" if (args.order == 1 or args.teno) else "T", + "null_weights": "F", + "mp_weno": "F", + "teno": "T" if args.teno else "F", + **({"teno_CT": args.teno_ct} if args.teno else {}), + } + +print( + json.dumps( + { + "run_time_info": "F", + "x_domain%beg": 0.0, + "x_domain%end": L, + "m": m, + "n": 0, + "p": 0, + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.25, + "patch_icpp(1)%length_x": 0.5, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": 1.0, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.75, + "patch_icpp(2)%length_x": 0.5, + "patch_icpp(2)%vel(1)": 0.0, + "patch_icpp(2)%pres": 0.1, + "patch_icpp(2)%alpha_rho(1)": 0.125, + "patch_icpp(2)%alpha(1)": 1.0, + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%pi_inf": 0.0, + **scheme_params, + } + ) +) diff --git a/examples/2D_advection_convergence/case.py b/examples/2D_advection_convergence/case.py new file mode 100644 index 0000000000..39ce6f1c1e --- /dev/null +++ b/examples/2D_advection_convergence/case.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +""" +2D single-fluid Euler convergence case (smooth diagonal advection). + +Density wave rho = 1 + 0.2*sin(2*pi*(x+y)) with uniform velocity (u, v) = (1, 1) +on the unit square with periodic BCs. Pressure p = 1 throughout. For this IC +the Euler equations reduce to pure advection of all variables along the +diagonal at speed sqrt(2); the wave phase x+y - 2t returns to x+y after +T = 0.5. L2(rho(T) - rho(0)) measures accumulated scheme spatial truncation +error. + +The point of this case is to exercise WENO7 / TENO7 in 2D without the +primitive-to-conserved covariance floor that dominates the 2D isentropic +vortex (run_convergence.py) at testable resolutions. Here all primitives are +linear in the conserved variables, so there is no floor. +""" + +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="2D Euler diagonal-advection convergence case") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=64, help="Grid points per dim (default: 64)") +parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, 5, or 7") +parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO") +parser.add_argument("--teno", action="store_true", help="Use TENO instead of WENO") +parser.add_argument("--teno-ct", type=float, default=1e-6, help="TENO CT threshold (default: 1e-6)") +parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4; use 0.005 for WENO7/TENO7)") +parser.add_argument("--no-mapped", action="store_true", help="Disable mapped WENO") +parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)") +parser.add_argument("--time-stepper", type=int, default=3, help="Time stepper: 1=Euler 2=RK2 3=RK3 (default: 3)") +parser.add_argument("--t-end", type=float, default=None, help="Override total simulation time (default: 0.5 = one diagonal period)") +args = parser.parse_args() + +gamma = 1.4 +N = args.N +m = N - 1 +L = 1.0 +dx = L / N + +# Max wave speed: |u| + c with c = sqrt(gamma * p / rho_min). rho_min = 1 - 0.2. +c_max = math.sqrt(gamma / 0.8) + 1.0 +dt = args.cfl * dx / c_max +T_end = args.t_end if args.t_end is not None else 0.5 +Nt = max(1, math.ceil(T_end / dt)) +dt = T_end / Nt + +if args.muscl: + scheme_params = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": args.muscl_lim, + } +else: + scheme_params = { + "recon_type": 1, + "weno_order": args.order, + "weno_eps": 1.0e-40, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F" if (args.order == 1 or args.no_mapped or args.teno) else "T", + "null_weights": "F", + "mp_weno": "F", + "teno": "T" if args.teno else "F", + **({"teno_CT": args.teno_ct} if args.teno else {}), + } + +print( + json.dumps( + { + "run_time_info": "F", + "x_domain%beg": 0.0, + "x_domain%end": L, + "y_domain%beg": 0.0, + "y_domain%end": L, + "m": m, + "n": m, + "p": 0, + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": args.time_stepper, + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": L, + "patch_icpp(1)%length_y": L, + "patch_icpp(1)%vel(1)": 1.0, + "patch_icpp(1)%vel(2)": 1.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": "1.0 + 0.2 * sin(2.0 * pi * (x / lx + y / ly))", + "patch_icpp(1)%alpha(1)": 1.0, + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%pi_inf": 0.0, + **scheme_params, + } + ) +) diff --git a/examples/3D_advection_convergence/case.py b/examples/3D_advection_convergence/case.py new file mode 100644 index 0000000000..eeb33b071b --- /dev/null +++ b/examples/3D_advection_convergence/case.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python3 +""" +3D single-fluid Euler convergence case (smooth diagonal advection). + +Density wave rho = 1 + 0.2*sin(2*pi*(x+y+z)) with uniform velocity (1, 1, 1) +on the unit cube with periodic BCs. Pressure p = 1 throughout. The wave +phase x+y+z - 3t returns to x+y+z after T = 1/3. + +Usually invoked from the convergence test framework with --t-end set so +T = K / N for some integer K (cell-shift mode), allowing analytical +comparison via integer-cell np.roll of the IC. +""" + +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="3D Euler diagonal-advection convergence case") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=32, help="Grid points per dim (default: 32)") +parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, 5, or 7") +parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO") +parser.add_argument("--teno", action="store_true", help="Use TENO instead of WENO") +parser.add_argument("--teno-ct", type=float, default=1e-6, help="TENO CT threshold (default: 1e-6)") +parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4; use 0.005 for WENO7/TENO7)") +parser.add_argument("--no-mapped", action="store_true", help="Disable mapped WENO") +parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)") +parser.add_argument("--time-stepper", type=int, default=3, help="Time stepper: 1=Euler 2=RK2 3=RK3 (default: 3)") +parser.add_argument("--t-end", type=float, default=None, help="Override total simulation time (default: 1/3 = one diagonal period)") +args = parser.parse_args() + +gamma = 1.4 +N = args.N +m = N - 1 +L = 1.0 +dx = L / N + +c_max = math.sqrt(gamma / 0.8) + 1.0 +dt = args.cfl * dx / c_max +T_end = args.t_end if args.t_end is not None else 1.0 / 3.0 +Nt = max(1, math.ceil(T_end / dt)) +dt = T_end / Nt + +if args.muscl: + scheme_params = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": args.muscl_lim, + } +else: + scheme_params = { + "recon_type": 1, + "weno_order": args.order, + "weno_eps": 1.0e-40, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F" if (args.order == 1 or args.no_mapped or args.teno) else "T", + "null_weights": "F", + "mp_weno": "F", + "teno": "T" if args.teno else "F", + **({"teno_CT": args.teno_ct} if args.teno else {}), + } + +print( + json.dumps( + { + "run_time_info": "F", + "x_domain%beg": 0.0, + "x_domain%end": L, + "y_domain%beg": 0.0, + "y_domain%end": L, + "z_domain%beg": 0.0, + "z_domain%end": L, + "m": m, + "n": m, + "p": m, + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": args.time_stepper, + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + "bc_z%beg": -1, + "bc_z%end": -1, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%z_centroid": 0.5, + "patch_icpp(1)%length_x": L, + "patch_icpp(1)%length_y": L, + "patch_icpp(1)%length_z": L, + "patch_icpp(1)%vel(1)": 1.0, + "patch_icpp(1)%vel(2)": 1.0, + "patch_icpp(1)%vel(3)": 1.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": "1.0 + 0.2 * sin(2.0 * pi * (x / lx + y / ly + z / lz))", + "patch_icpp(1)%alpha(1)": 1.0, + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%pi_inf": 0.0, + **scheme_params, + } + ) +) diff --git a/src/common/include/2dHardcodedIC.fpp b/src/common/include/2dHardcodedIC.fpp index 3f4dc4bcbb..d6c3fed97f 100644 --- a/src/common/include/2dHardcodedIC.fpp +++ b/src/common/include/2dHardcodedIC.fpp @@ -8,6 +8,12 @@ real(wp) :: sinA, cosA real(wp) :: r_sq + ! # 283 - Gauss-averaged isentropic vortex (conserved-variable cell averages) + real(wp) :: gauss_xi(3), gauss_w(3), xq, yq, r2q, T_facq, wq + real(wp) :: rho_avg, rhou_avg, rhov_avg, E_avg + real(wp) :: rhoq, pq, uq, vq, Eq, vortex_eps + integer :: igq, jgq + ! # 291 - Shear/Thermal Layer Case real(wp) :: delta_shear, u_max, u_mean real(wp) :: T_wall, T_inf, P_atm, T_loc @@ -285,11 +291,11 @@ & 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*pi))*(5.0/(8.0*1.0*(1.4 + 1.0)*pi))*exp(2.0*1.0*(1.0 - (x_cc(i) & & - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**1.4 q_prim_vf(eqn_idx%mom%beg + 0)%sf(i, j, & - & 0) = 0.0 + (y_cc(j) - patch_icpp(1)%y_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1) & - & %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) + & 0) = patch_icpp(1)%vel(1) + (y_cc(j) - patch_icpp(1)%y_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) & + & - patch_icpp(1) %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, & - & 0) = 0.0 - (x_cc(i) - patch_icpp(1)%x_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1) & - & %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) + & 0) = patch_icpp(1)%vel(2) - (x_cc(i) - patch_icpp(1)%x_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) & + & - patch_icpp(1) %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) end if case (281) ! Acoustic pulse ! This is patch is hard-coded for test suite optimization used in the 2D_acoustic_pulse case: This analytic patch uses @@ -313,6 +319,46 @@ q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, & & 0) = 112.99092883944267*((0.1/0.3))*x_cc(i)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))) end if + case (283) ! Isentropic vortex: conserved-variable GL cell averages (3-pt tensor product) + ! GL averages of conserved variables (rho, rho*u, rho*v, E) eliminate the O(h^2) error that primitive-variable averaging + ! introduces through the nonlinear prim->cons conversion: cell_avg(rho*u) != cell_avg(rho)*cell_avg(u) by O(h^2). We back + ! out primitive values that reproduce the conserved averages exactly. Vortex strength eps is read from + ! patch_icpp(patch_id)%epsilon; defaults to 5. + if (patch_id == 1) then + vortex_eps = merge(patch_icpp(patch_id)%epsilon, 5._wp, patch_icpp(patch_id)%epsilon > 0._wp) + gauss_xi = [-sqrt(3._wp/5._wp), 0._wp, sqrt(3._wp/5._wp)] + gauss_w = [5._wp/9._wp, 8._wp/9._wp, 5._wp/9._wp] + rho_avg = 0._wp; rhou_avg = 0._wp; rhov_avg = 0._wp; E_avg = 0._wp + do igq = 1, 3 + do jgq = 1, 3 + xq = x_cc(i) + gauss_xi(igq)*(x_cb(i) - x_cb(i - 1))*0.5_wp + yq = y_cc(j) + gauss_xi(jgq)*(y_cb(j) - y_cb(j - 1))*0.5_wp + r2q = (xq - patch_icpp(patch_id)%x_centroid)**2._wp + (yq - patch_icpp(patch_id)%y_centroid)**2._wp + T_facq = 1._wp - (vortex_eps/(2._wp*pi))*(vortex_eps/(8._wp*(1.4_wp + 1._wp)*pi))*exp(2._wp*(1._wp - r2q)) + wq = gauss_w(igq)*gauss_w(jgq) + rhoq = T_facq**1.4_wp + pq = T_facq**2.4_wp + uq = patch_icpp(patch_id)%vel(1) + (yq - patch_icpp(patch_id)%y_centroid)*(vortex_eps/(2._wp*pi))*exp(1._wp & + & - r2q) + vq = patch_icpp(patch_id)%vel(2) - (xq - patch_icpp(patch_id)%x_centroid)*(vortex_eps/(2._wp*pi))*exp(1._wp & + & - r2q) + Eq = pq/0.4_wp + 0.5_wp*rhoq*(uq**2 + vq**2) + rho_avg = rho_avg + wq*rhoq + rhou_avg = rhou_avg + wq*(rhoq*uq) + rhov_avg = rhov_avg + wq*(rhoq*vq) + E_avg = E_avg + wq*Eq + end do + end do + rho_avg = rho_avg*0.25_wp + rhou_avg = rhou_avg*0.25_wp + rhov_avg = rhov_avg*0.25_wp + E_avg = E_avg*0.25_wp + ! Back out primitive vars so prim->cons conversion recovers the conserved averages + q_prim_vf(eqn_idx%cont%beg)%sf(i, j, 0) = rho_avg + q_prim_vf(eqn_idx%mom%beg + 0)%sf(i, j, 0) = rhou_avg/rho_avg + q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, 0) = rhov_avg/rho_avg + q_prim_vf(eqn_idx%E)%sf(i, j, 0) = (E_avg - 0.5_wp*(rhou_avg**2 + rhov_avg**2)/rho_avg)*0.4_wp + end if case (291) ! Isothermal Flat Plate T_inf = 1125.0_wp T_wall = 600.0_wp diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp index aca419474b..c8eff18847 100644 --- a/src/simulation/m_muscl.fpp +++ b/src/simulation/m_muscl.fpp @@ -171,7 +171,9 @@ contains slopeR = v_rs_ws_${XYZ}$_muscl(j, k, l, i) - v_rs_ws_${XYZ}$_muscl(j - 1, k, l, i) slope = 0._wp - if (muscl_lim == 1) then ! minmod + if (muscl_lim == 0) then ! unlimited (central difference) + slope = 5e-1_wp*(slopeL + slopeR) + else if (muscl_lim == 1) then ! minmod if (slopeL*slopeR > muscl_eps) then slope = min(abs(slopeL), abs(slopeR)) end if diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index d96d834755..05a3fee92c 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -133,7 +133,7 @@ "check_muscl": { "title": "MUSCL Reconstruction", "category": "Numerical Schemes", - "explanation": "muscl_order must be 1 or 2. Second order requires muscl_lim in {1,2,3,4,5}.", + "explanation": "muscl_order must be 1 or 2. Second order requires muscl_lim in {0,1,2,3,4,5}.", }, "check_time_stepping": { "title": "Time Stepping", @@ -808,7 +808,7 @@ def check_muscl_simulation(self): return self.prohibit(muscl_order == 2 and muscl_lim is None, "muscl_lim must be defined if using muscl_order = 2") - self.prohibit(muscl_lim is not None and (muscl_lim < 1 or muscl_lim > 5), "muscl_lim must be 1, 2, 3, 4, or 5") + self.prohibit(muscl_lim is not None and (muscl_lim < 0 or muscl_lim > 5), "muscl_lim must be 0 (unlimited), 1, 2, 3, 4, or 5") if muscl_eps is not None: self.prohibit(muscl_eps < 0, "muscl_eps must be >= 0 (use 0 for textbook limiter behavior)") diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index ccc4b10a73..18f107df74 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -609,8 +609,8 @@ def get_value_label(param_name: str, value: int) -> str: "value_labels": {1: "1st order", 2: "2nd order"}, }, "muscl_lim": { - "choices": [1, 2, 3, 4, 5], - "value_labels": {1: "minmod", 2: "MC", 3: "Van Albada", 4: "Van Leer", 5: "SUPERBEE"}, + "choices": [0, 1, 2, 3, 4, 5], + "value_labels": {0: "unlimited", 1: "minmod", 2: "MC", 3: "Van Albada", 4: "Van Leer", 5: "SUPERBEE"}, }, "int_comp": { "choices": [0, 1, 2], diff --git a/toolchain/mfc/test/case.py b/toolchain/mfc/test/case.py index c8cf0c2485..05c4bfa0bc 100644 --- a/toolchain/mfc/test/case.py +++ b/toolchain/mfc/test/case.py @@ -155,12 +155,16 @@ class TestCase(case.Case): trace: str override_tol: Optional[float] = None restart_check: bool = False + kind: str = "golden" + convergence_spec: Optional[dict] = None - def __init__(self, trace: str, mods: dict, ppn: int = None, override_tol: float = None, restart_check: bool = False) -> None: + def __init__(self, trace: str, mods: dict, ppn: int = None, override_tol: float = None, restart_check: bool = False, kind: str = "golden", convergence_spec: Optional[dict] = None) -> None: self.trace = trace self.ppn = ppn or 1 self.override_tol = override_tol self.restart_check = restart_check + self.kind = kind + self.convergence_spec = convergence_spec merge = {**BASE_CFG.copy(), **mods} merge = {key: val for key, val in merge.items() if val is not None} super().__init__(merge) @@ -361,11 +365,18 @@ class TestCaseBuilder: functor: Optional[Callable] override_tol: Optional[float] = None restart_check: bool = False + kind: str = "golden" + convergence_spec: Optional[dict] = None def get_uuid(self) -> str: return trace_to_uuid(self.trace) def to_case(self) -> TestCase: + if self.kind == "convergence": + # Convergence cases drive their own runs — the BASE_CFG mods/path + # machinery is unused. Trace + spec are the only inputs. + return TestCase(self.trace, {}, self.ppn, self.override_tol, self.restart_check, kind=self.kind, convergence_spec=self.convergence_spec) + dictionary = {} if self.path: dictionary.update(input.load(self.path, self.args, do_print=False).params) @@ -411,6 +422,15 @@ def define_case_f(trace: str, path: str, args: List[str] = None, ppn: int = None return TestCaseBuilder(trace, mods or {}, path, args or [], ppn or 1, functor, override_tol) +def define_convergence_case(trace: str, spec: dict, ppn: int = None) -> TestCaseBuilder: + """Register a convergence-rate test (kind='convergence'). + + `spec` must include a `runner` key naming a runner in + `toolchain/mfc/test/convergence.py` plus runner-specific arguments. + """ + return TestCaseBuilder(trace, {}, None, None, ppn or 1, None, None, False, kind="convergence", convergence_spec=spec) + + def define_case_d(stack: CaseGeneratorStack, newTrace: str, newMods: dict, ppn: int = None, functor: Callable = None, override_tol: float = None, restart_check: bool = False) -> TestCaseBuilder: mods: dict = {} diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index 5b79ffbb32..b72584a032 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -5,7 +5,161 @@ from mfc import common from ..state import ARG -from .case import CaseGeneratorStack, Nt, TestCaseBuilder, define_case_d, define_case_f +from .case import CaseGeneratorStack, Nt, TestCaseBuilder, define_case_d, define_case_f, define_convergence_case +from .convergence import ConvergenceSpec, run_dt_sweep, run_h_sweep, run_sod_l1 + +# Convergence test specs. +# One TestCase per (problem, scheme) pair. Trace prefix "Convergence ->" is +# the filter handle (`./mfc.sh test --only Convergence`); convergence cases +# are skipped by default. + +# Advection convergence cases. Cell-shift mode: T = K*h per resolution, compare +# q(T) to np.roll(q(0), +K) per dim. Cost is O(1) in N (Nt = K*c/CFL independent +# of resolution) — wins ~10-100x vs period mode (full advection period). +# WENO7/TENO7 stay in period mode: at typical N their cell-shift signal sinks +# below machine precision (h^8 < 1e-15 at N=64) before any rate develops. +# +# expected_order is always the scheme's spatial order p. The runner subtracts +# 1 from the displayed rate in cell-shift mode (where raw rate = p+1) so the +# reported "spatial order" matches expected_order in both modes. +_CONS_VARS_1D = [("density", 1), ("x-momentum", 2), ("energy", 3)] +_CONS_VARS_2D = [("density", 1), ("energy", 4)] +_CONS_VARS_3D = [("density", 1), ("energy", 5)] + +# (label, extra_args, expected_order, tol, resolutions) +# WENO3-JS at smooth extrema empirically achieves ~1.5 in MFC (Henrick mapping +# enabled). MUSCL2 unlimited central → effective spatial order 2. +_CONVERGENCE_1D_SCHEMES = [ + ("WENO5", ["--order", "5", "--cfl", "0.02"], 5, 0.3, [32, 64, 128]), + ("WENO3", ["--order", "3", "--cfl", "0.02"], 1.5, 0.3, [64, 128, 256]), + ("WENO1", ["--order", "1", "--cfl", "0.02"], 1, 0.2, [64, 128, 256]), + ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 2, 0.3, [64, 128, 256]), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 5, 0.3, [32, 64, 128]), +] +# WENO7/TENO7 in 1D: period mode (full period T=1.0, see 1D case.py). +_CONVERGENCE_1D_PERIOD_SCHEMES = [ + ("WENO7", ["--order", "7", "--cfl", "0.005"], 7, 0.5, [64, 128]), + ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9", "--cfl", "0.005"], 7, 0.5, [64, 128]), +] + +_CONVERGENCE_2D_SCHEMES = [ + ("WENO5", ["--order", "5", "--cfl", "0.02"], 5, 0.3, [32, 64, 96]), + ("WENO3", ["--order", "3", "--cfl", "0.02"], 1.5, 0.3, [32, 64, 128]), + ("WENO1", ["--order", "1", "--cfl", "0.02"], 1, 0.2, [32, 64, 128]), + ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 2, 0.3, [32, 64, 128]), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 5, 0.3, [32, 64, 96]), +] +_CONVERGENCE_2D_PERIOD_SCHEMES = [ + ("WENO7", ["--order", "7", "--cfl", "0.005"], 7, 0.5, [80, 96]), + ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9", "--cfl", "0.005"], 7, 0.5, [80, 96]), +] + +# 3D diagonal advection: only cell-shift mode (period T=1/3 with N^3 cells +# would dominate CI even at N=64). WENO7/TENO7 skipped — at N=64 with K=1 +# the spatial error signal is below machine precision. +_CONVERGENCE_3D_SCHEMES = [ + ("WENO5", ["--order", "5", "--cfl", "0.02"], 5, 0.3, [32, 64]), + ("WENO3", ["--order", "3", "--cfl", "0.02"], 1.5, 0.3, [32, 64]), + ("WENO1", ["--order", "1", "--cfl", "0.02"], 1, 0.2, [32, 64]), + ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 2, 0.3, [32, 64]), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 5, 0.3, [32, 64]), +] + +# Sod L1 self-convergence: any conservative monotone scheme converges at L1 +# rate ~1 (Godunov). SUPERBEE is over-compressive; min_N=128 skips its +# pre-asymptotic point. +_RES_SOD_DEFAULT = [128, 256, 512, 1024] +_CONVERGENCE_SOD_SCHEMES = [ + ("WENO1", ["--order", "1"], 1, 0.5, None), + ("WENO3", ["--order", "3"], 1, 0.3, None), + ("WENO5", ["--order", "5"], 1, 0.3, None), + ("WENO7", ["--order", "7"], 1, 0.3, None), + ("MUSCL-minmod", ["--muscl", "--muscl-lim", "1"], 1, 0.3, None), + ("MUSCL-MC", ["--muscl", "--muscl-lim", "2"], 1, 0.3, None), + ("MUSCL-VanLeer", ["--muscl", "--muscl-lim", "4"], 1, 0.3, None), + ("MUSCL-SUPERBEE", ["--muscl", "--muscl-lim", "5"], 1, 0.5, 128), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6"], 1, 0.3, None), + ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9"], 1, 0.3, None), +] + +# Temporal order: fixed N=512 / WENO5; vary CFL. +_CONVERGENCE_TEMPORAL_SCHEMES = [ + ("RK1", ["--order", "5", "--time-stepper", "1"], 1, 0.1, [0.10, 0.05]), + ("RK2", ["--order", "5", "--time-stepper", "2"], 2, 0.2, [0.50, 0.25]), + ("RK3", ["--order", "5", "--time-stepper", "3"], 3, 0.3, [0.50, 0.25]), +] + + +def add_convergence_cases(cases): + num_ranks = 4 + + def _h_sweep(case_path, ndim, cons_vars, extra_args, expected, tol, resolutions, cell_shift): + return ConvergenceSpec( + runner=run_h_sweep, + case_path=case_path, + extra_args=extra_args, + expected_order=expected, + tol=tol, + cons_vars=cons_vars, + resolutions=resolutions, + ndim=ndim, + cell_shift=cell_shift, + num_ranks=num_ranks, # ignored by run_h_sweep when cell_shift > 0 + ) + + advection_groups = [ + (_CONVERGENCE_1D_SCHEMES, "1D", "examples/1D_euler_convergence/case.py", 1, _CONS_VARS_1D, 1, 1), + (_CONVERGENCE_1D_PERIOD_SCHEMES, "1D", "examples/1D_euler_convergence/case.py", 1, _CONS_VARS_1D, 0, num_ranks), + (_CONVERGENCE_2D_SCHEMES, "2D", "examples/2D_advection_convergence/case.py", 2, _CONS_VARS_2D, 1, 1), + (_CONVERGENCE_2D_PERIOD_SCHEMES, "2D", "examples/2D_advection_convergence/case.py", 2, _CONS_VARS_2D, 0, num_ranks), + (_CONVERGENCE_3D_SCHEMES, "3D", "examples/3D_advection_convergence/case.py", 3, _CONS_VARS_3D, 1, 1), + ] + for schemes, dim_label, case_path, ndim, cons_vars, cell_shift, ppn in advection_groups: + for label, extra_args, expected, tol, resolutions in schemes: + cases.append( + define_convergence_case( + f"Convergence -> {dim_label} -> {label}", + spec=_h_sweep(case_path, ndim, cons_vars, extra_args, expected, tol, resolutions, cell_shift), + ppn=ppn, + ) + ) + + for label, extra_args, expected, tol, min_N in _CONVERGENCE_SOD_SCHEMES: + resolutions = [N for N in _RES_SOD_DEFAULT if min_N is None or N >= min_N] + cases.append( + define_convergence_case( + f"Convergence -> Sod -> {label}", + spec=ConvergenceSpec( + runner=run_sod_l1, + case_path="examples/1D_sod_convergence/case.py", + extra_args=extra_args, + expected_order=expected, + tol=tol, + resolutions=resolutions, + num_ranks=num_ranks, + ), + ppn=num_ranks, + ) + ) + + for label, extra_args, expected, tol, cfls in _CONVERGENCE_TEMPORAL_SCHEMES: + cases.append( + define_convergence_case( + f"Convergence -> Temporal -> {label}", + spec=ConvergenceSpec( + runner=run_dt_sweep, + case_path="examples/1D_euler_convergence/case.py", + extra_args=extra_args, + expected_order=expected, + tol=tol, + cons_vars=_CONS_VARS_1D, + cfls=cfls, + fixed_N=512, + num_ranks=num_ranks, + ), + ppn=num_ranks, + ) + ) def get_bc_mods(bc: int, dimInfo): @@ -1581,6 +1735,11 @@ def foreach_example(): "1D_titarevtorro_analytical", "2D_acoustic_pulse_analytical", "2D_isentropicvortex_analytical", + "1D_euler_convergence", + "1D_advection_convergence", + "1D_sod_convergence", + "2D_advection_convergence", + "3D_advection_convergence", "2D_zero_circ_vortex_analytical", "3D_TaylorGreenVortex_analytical", "3D_IGR_TaylorGreenVortex_nvidia", @@ -2232,6 +2391,8 @@ def kernel_golden_tests(): kernel_golden_tests() + add_convergence_cases(cases) + # Sanity Check 1 if stack.size() != 0: raise common.MFCException("list_cases: stack isn't fully pop'ed") diff --git a/toolchain/mfc/test/convergence.py b/toolchain/mfc/test/convergence.py new file mode 100644 index 0000000000..d6512f75ea --- /dev/null +++ b/toolchain/mfc/test/convergence.py @@ -0,0 +1,324 @@ +"""Convergence-rate tests for ./mfc.sh test. + +A "convergence case" runs MFC at several resolutions (or several CFLs), reads +back the conserved variable, fits log(error) vs log(h) for the rate, and +checks that the rate matches the scheme's nominal order. Each case is one +TestCase with kind="convergence" and a ConvergenceSpec attached; test.py +routes it here via run_case(). + +Three runner flavours, picked by setting spec.runner in cases.py: + + run_h_sweep vary N (fixed CFL). cell_shift > 0 makes T = K*h and compares + to np.roll(q(0), +K) — cheap (Nt = O(1) in N) and reports the + scheme spatial order p directly. cell_shift == 0 runs to one + full period and compares to q(0). + run_dt_sweep fix N, vary CFL. Measures the time-stepper order. + run_sod_l1 1D L1 self-convergence: compare each N against 2N after + cell-averaging the fine grid. + +Helpers below: _read_field reads a Fortran-unformatted record per rank and +concatenates; _l2_norm and _fit_slope are obvious one-liners; _run_mfc shells +out to ./mfc.sh run and stashes p_all/ in a tempdir. +""" + +import dataclasses +import json +import math +import os +import shutil +import struct +import subprocess +import sys +import tempfile +import typing + +import numpy as np + +from .. import common + +CONS_TOL = 1e-10 +MFC = ".\\mfc.bat" if os.name == "nt" else "./mfc.sh" + + +@dataclasses.dataclass +class ConvergenceSpec: + """One convergence case. `runner` is the function that drives it.""" + + runner: typing.Callable + case_path: str + expected_order: float # scheme order p (or temporal q for dt sweeps) + tol: float # passes if measured >= expected_order - tol + extra_args: typing.List[str] = dataclasses.field(default_factory=list) + cons_vars: typing.List[typing.Tuple[str, int]] = dataclasses.field(default_factory=list) + primary_idx: int = 1 # which q_cons_vf{idx} to use for the L2/L1 norm + num_ranks: int = 1 + # Spatial-sweep cases: + resolutions: typing.List[int] = dataclasses.field(default_factory=list) + ndim: int = 1 + domain_len: float = 1.0 + cell_shift: int = 0 # 0 -> period mode; >0 -> compare to K-cell-shifted IC + # Temporal-sweep cases: + cfls: typing.List[float] = dataclasses.field(default_factory=list) + fixed_N: int = 0 + + +# Low-level helpers. + + +def _read_field(run_dir: str, step: int, var_idx: int, num_ranks: int, expected_size: int) -> np.ndarray: + """Read q_cons_vf.dat across all ranks, in rank order.""" + chunks = [] + for rank in range(num_ranks): + path = os.path.join(run_dir, "p_all", f"p{rank}", str(step), f"q_cons_vf{var_idx}.dat") + with open(path, "rb") as f: + rec_len = struct.unpack("i", f.read(4))[0] + chunks.append(np.frombuffer(f.read(rec_len), dtype=np.float64).copy()) + f.read(4) + arr = np.concatenate(chunks) + if arr.size != expected_size: + raise common.MFCException(f"Expected {expected_size} cells, got {arr.size}") + return arr + + +def _l2_norm(diff: np.ndarray, cell_vol: float) -> float: + return float(np.sqrt(np.sum(diff**2) * cell_vol)) + + +def _fit_slope(errors: typing.List[float], xs: typing.List[float]) -> float: + """Least-squares slope of log(errors) vs log(xs).""" + return float(np.polyfit(np.log(xs), np.log(errors), 1)[0]) + + +def _pairwise_slope(err1: float, err2: float, x1: float, x2: float) -> float: + return math.log(err2 / err1) / math.log(x2 / x1) + + +def _run_mfc(case_path: str, tmpdir: str, run_tag: str, args: typing.List[str], num_ranks: int) -> typing.Tuple[dict, str]: + """Run case.py through ./mfc.sh and stash p_all/ in tmpdir/run_tag for reading.""" + cfg_run = subprocess.run( + [sys.executable, case_path, "--mfc", "{}"] + args, + capture_output=True, + text=True, + check=False, + ) + if cfg_run.returncode != 0: + raise common.MFCException(f"case.py failed:\n{cfg_run.stderr}") + cfg = json.loads(cfg_run.stdout) + + sim = subprocess.run( + [MFC, "run", case_path, "-t", "pre_process", "simulation", "-n", str(num_ranks), "--"] + args, + capture_output=True, + text=True, + check=False, + ) + if sim.returncode != 0: + raise common.MFCException(f"./mfc.sh run failed for {run_tag}\n{sim.stdout[-3000:]}\n{sim.stderr}") + + case_dir = os.path.dirname(case_path) + src = os.path.join(case_dir, "p_all") + dst = os.path.join(tmpdir, run_tag, "p_all") + if os.path.exists(dst): + shutil.rmtree(dst) + shutil.copytree(src, dst) + shutil.rmtree(src, ignore_errors=True) + shutil.rmtree(os.path.join(case_dir, "D"), ignore_errors=True) + return cfg, os.path.join(tmpdir, run_tag) + + +def _conservation_at_step(run_dir: str, Nt: int, cell_vol: float, cons_vars, num_ranks: int, cell_count: int) -> dict: + """Per-variable |∫q(T) - ∫q(0)| / |∫q(0)|.""" + out = {} + for name, idx in cons_vars: + s0 = float(np.sum(_read_field(run_dir, 0, idx, num_ranks, cell_count))) * cell_vol + sT = float(np.sum(_read_field(run_dir, Nt, idx, num_ranks, cell_count))) * cell_vol + out[name] = abs(sT - s0) / (abs(s0) + 1e-300) + return out + + +def _conservation_lines(history: typing.List[dict], cons_vars) -> typing.Tuple[bool, typing.List[str]]: + lines = [f"\n Conservation (need rel. error < {CONS_TOL:.0e}):"] + passed = True + for name, _ in cons_vars: + max_err = max(c[name] for c in history) + ok = max_err < CONS_TOL + passed = passed and ok + lines.append(f" {name:<14}: max = {max_err:.2e} {'OK' if ok else 'FAIL'}") + return passed, lines + + +def _table_line(values: typing.List[str], widths: typing.List[int]) -> str: + return " " + " ".join(f"{v:>{w}}" for v, w in zip(values, widths)) + + +# Runners. + + +def run_h_sweep(spec: ConvergenceSpec) -> typing.Tuple[bool, str]: + """Vary N at fixed CFL. + + cell_shift > 0: T = K*h with v=1 — compare q(T) to np.roll(q(0), +K). + Cost is O(1) in N; raw rate is p+1, displayed as scheme order p. + cell_shift == 0: T = full period — compare q(T) to q(0). Rate = p. + """ + is_shift = spec.cell_shift > 0 + num_ranks = 1 if is_shift else spec.num_ranks # cell-shift reshape needs single rank + label = "spatial order" if is_shift else "rate" + threshold = spec.expected_order - spec.tol + + errors, hs, nts, history = [], [], [], [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in spec.resolutions: + h = spec.domain_len / N + cell_vol = h**spec.ndim + cell_count = N**spec.ndim + args = ["-N", str(N)] + spec.extra_args + if is_shift: + args += ["--t-end", str(spec.cell_shift * h)] + + cfg, run_dir = _run_mfc(spec.case_path, tmpdir, f"N{N}", args, num_ranks) + Nt = int(cfg["t_step_stop"]) + q0 = _read_field(run_dir, 0, spec.primary_idx, num_ranks, cell_count) + qT = _read_field(run_dir, Nt, spec.primary_idx, num_ranks, cell_count) + ref = q0 + if is_shift: + # Wave moves at v=+1: q(T)[i] = q(0)[i-K], so np.roll shift=+K. + shape = (N,) * spec.ndim + ref = np.roll( + q0.reshape(shape, order="F"), + shift=spec.cell_shift, + axis=tuple(range(spec.ndim)), + ).flatten(order="F") + + errors.append(_l2_norm(qT - ref, cell_vol)) + hs.append(h) + nts.append(Nt) + history.append(_conservation_at_step(run_dir, Nt, cell_vol, spec.cons_vars, num_ranks, cell_count)) + + offset = 1 if is_shift else 0 + widths = [6, 6, 10, 14, 13] + lines = [ + f" (need {label} >= {threshold:.1f})", + "", + _table_line(["N", "Nt", "dx", "L2 error", label.replace(" ", "_")], widths), + _table_line(["-" * w for w in widths], widths), + ] + for i, N in enumerate(spec.resolutions): + if i == 0: + r_str = "---" + else: + r_str = f"{_pairwise_slope(errors[i - 1], errors[i], hs[i - 1], hs[i]) - offset:.2f}" + lines.append(_table_line([str(N), str(nts[i]), f"{hs[i]:.6f}", f"{errors[i]:.6e}", r_str], widths)) + + rate_passed = True + if len(errors) > 1: + fitted = _fit_slope(errors, hs) - offset + lines.append(f"\n Fitted {label}: {fitted:.2f} (need >= {threshold:.1f})") + rate_passed = fitted >= threshold + + cons_passed, cons_lines = _conservation_lines(history, spec.cons_vars) + lines += cons_lines + return rate_passed and cons_passed, "\n".join(lines) + + +def run_dt_sweep(spec: ConvergenceSpec) -> typing.Tuple[bool, str]: + """Fix N, vary CFL — measures temporal scheme order.""" + threshold = spec.expected_order - spec.tol + h = 1.0 / spec.fixed_N + + errors, dts, nts, history = [], [], [], [] + with tempfile.TemporaryDirectory() as tmpdir: + for cfl in spec.cfls: + tag = f"cfl{cfl:.4f}".replace(".", "p") + args = ["-N", str(spec.fixed_N), "--cfl", str(cfl)] + spec.extra_args + cfg, run_dir = _run_mfc(spec.case_path, tmpdir, tag, args, spec.num_ranks) + Nt = int(cfg["t_step_stop"]) + q0 = _read_field(run_dir, 0, spec.primary_idx, spec.num_ranks, spec.fixed_N) + qT = _read_field(run_dir, Nt, spec.primary_idx, spec.num_ranks, spec.fixed_N) + errors.append(_l2_norm(qT - q0, h)) + dts.append(float(cfg["dt"])) + nts.append(Nt) + history.append(_conservation_at_step(run_dir, Nt, h, spec.cons_vars, spec.num_ranks, spec.fixed_N)) + + widths = [7, 12, 6, 14, 8] + lines = [ + f" N={spec.fixed_N} (need rate >= {threshold:.1f})", + "", + _table_line(["CFL", "dt", "Nt", "L2 error", "rate"], widths), + _table_line(["-" * w for w in widths], widths), + ] + for i, cfl in enumerate(spec.cfls): + if i == 0: + r_str = "---" + else: + r_str = f"{_pairwise_slope(errors[i - 1], errors[i], dts[i - 1], dts[i]):.2f}" + lines.append(_table_line([f"{cfl:.3f}", f"{dts[i]:.6e}", str(nts[i]), f"{errors[i]:.6e}", r_str], widths)) + + rate_passed = True + if len(errors) > 1: + fitted = _fit_slope(errors, dts) + lines.append(f"\n Fitted rate: {fitted:.2f} (need >= {threshold:.1f})") + rate_passed = fitted >= threshold + + cons_passed, cons_lines = _conservation_lines(history, spec.cons_vars) + lines += cons_lines + return rate_passed and cons_passed, "\n".join(lines) + + +def run_sod_l1(spec: ConvergenceSpec) -> typing.Tuple[bool, str]: + """1D L1 self-convergence: compare each N against 2N (fine cell-averaged 2:1).""" + threshold = spec.expected_order - spec.tol + + nts, run_dirs = [], [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in spec.resolutions: + cfg, run_dir = _run_mfc(spec.case_path, tmpdir, f"N{N}", ["-N", str(N)] + spec.extra_args, spec.num_ranks) + nts.append(int(cfg["t_step_stop"])) + run_dirs.append(run_dir) + + errors, error_resolutions = [], [] + for i in range(len(spec.resolutions) - 1): + N_c, N_f = spec.resolutions[i], spec.resolutions[i + 1] + if N_f != 2 * N_c: + continue + rho_c = _read_field(run_dirs[i], nts[i], 1, spec.num_ranks, N_c) + rho_f = _read_field(run_dirs[i + 1], nts[i + 1], 1, spec.num_ranks, N_f) + avg = 0.5 * (rho_f[0::2] + rho_f[1::2]) + errors.append(float(np.sum(np.abs(rho_c - avg)) * (1.0 / N_c))) + error_resolutions.append(N_c) + + hs = [1.0 / N for N in error_resolutions] + widths = [6, 6, 14, 8] + lines = [ + f" (need L1 rate >= {threshold:.1f})", + "", + _table_line(["N", "Nt", "L1 self-err", "rate"], widths), + _table_line(["-" * w for w in widths], widths), + ] + for i, N in enumerate(error_resolutions): + if i == 0: + r_str = "---" + else: + r_str = f"{_pairwise_slope(errors[i - 1], errors[i], hs[i - 1], hs[i]):.2f}" + lines.append(_table_line([str(N), str(nts[i]), f"{errors[i]:.6e}", r_str], widths)) + + if len(errors) >= 2: + fitted = _fit_slope(errors, hs) + lines.append(f"\n Fitted rate: {fitted:.2f} (need >= {threshold:.1f})") + return fitted >= threshold, "\n".join(lines) + if len(errors) == 1: + rate = _pairwise_slope(errors[0], errors[0], hs[0], hs[0]) + lines.append(f"\n Single pair rate (need >= {threshold:.1f})") + return rate >= threshold, "\n".join(lines) + lines.append("\n ERROR: need >= 2 consecutive 2x-apart resolutions") + return False, "\n".join(lines) + + +# Entry point used by test.py. + + +def run_case(spec: ConvergenceSpec) -> typing.Tuple[bool, str]: + """Run one convergence case. Returns (passed, full output report).""" + try: + return spec.runner(spec) + except Exception as exc: + return False, f" ERROR: {exc}\n" diff --git a/toolchain/mfc/test/test.py b/toolchain/mfc/test/test.py index 5e0d90d244..2d8ede6952 100644 --- a/toolchain/mfc/test/test.py +++ b/toolchain/mfc/test/test.py @@ -111,6 +111,18 @@ def __filter(cases_) -> typing.Tuple[typing.List[TestCase], typing.List[TestCase if not cases: raise MFCException(f"--only filter matched zero test cases. Specified: {ARG('only')}. Check that UUIDs/names are valid.") + # Convergence cases are slow (multiple resolutions × MPI ranks). Skip + # by default unless the user explicitly opted in via --only "Convergence" + # or a convergence UUID. _filter_only above has already narrowed cases + # to the user's selection, so any convergence case still present here + # was selected on purpose. Listing (`-l`) shows all cases regardless. + if not ARG("only") and not ARG("list"): + convergence_cases = [c for c in cases if getattr(c, "kind", "golden") == "convergence"] + if convergence_cases: + for c in convergence_cases: + cases.remove(c) + skipped_cases.append(c) + # --only-changes: filter based on file-level gcov coverage if ARG("only_changes"): from .coverage import ( @@ -435,10 +447,40 @@ def _process_silo_file(silo_filepath: str, case: TestCase, out_filepath: str): raise MFCException(f"Test {case}: Post Process has detected an Infinity. You can find the run's output in {out_filepath}, and the case dictionary in {case.get_filepath()}.") +def _handle_convergence_case(case: TestCase, start_time: float): + """Dispatch convergence/order-of-accuracy cases through convergence.py.""" + from .convergence import run_case + + if ARG("dry_run"): + trace_display = case.trace if len(case.trace) <= 50 else case.trace[:47] + "..." + cons.print(f" (dry-run) {trace_display:50s} SKIP [magenta]{case.get_uuid()}[/magenta]") + return + + passed, output = run_case(case.convergence_spec) + + log_dir = os.path.join(common.MFC_TEST_DIR, case.get_uuid()) + common.create_directory(log_dir) + common.file_write(os.path.join(log_dir, "convergence.log"), output) + + duration = time.time() - start_time + global current_test_number # noqa: PLW0603 + current_test_number += 1 + progress_str = f"({current_test_number:3d}/{total_test_count:3d})" + trace_display = case.trace if len(case.trace) <= 50 else case.trace[:47] + "..." + cons.print(f" {progress_str} {trace_display:50s} {duration:6.2f} [magenta]{case.get_uuid()}[/magenta]") + + if not passed: + raise MFCException(f"Test {case}: convergence rate check failed (see {log_dir}/convergence.log)") + + def _handle_case(case: TestCase, devices: typing.Set[int]): global current_test_number # noqa: PLW0603 start_time = time.time() + if getattr(case, "kind", "golden") == "convergence": + _handle_convergence_case(case, start_time) + return + # Set timeout using threading.Timer (works in worker threads) # Note: we intentionally do not use signal.alarm() here because signals # only work in the main thread; sched.sched runs tests in worker threads.