diff --git a/.github/workflows/convergence.yml b/.github/workflows/convergence.yml new file mode 100644 index 0000000000..d77c5786b5 --- /dev/null +++ b/.github/workflows/convergence.yml @@ -0,0 +1,154 @@ +name: Convergence + +on: + push: + branches: [master] + pull_request: + workflow_dispatch: + +env: + OMPI_MCA_rmaps_base_oversubscribe: 1 + +jobs: + convergence-1d: + name: "1D Advection Convergence" + runs-on: ubuntu-latest + timeout-minutes: 60 + + 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: Initialize MFC + run: ./mfc.sh init + + - name: Run 1D convergence tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_convergence_1d.py \ + --resolutions 64 128 256 512 1024 \ + --num-ranks 4 + + convergence-2d: + name: "2D Isentropic Vortex Convergence" + runs-on: ubuntu-latest + timeout-minutes: 60 + + 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: Initialize MFC + run: ./mfc.sh init + + - name: Run 2D convergence tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_convergence.py \ + --resolutions 32 64 128 \ + --num-ranks 4 + + convergence-sod: + name: "1D Sod Shock Tube L1 Convergence" + runs-on: ubuntu-latest + timeout-minutes: 90 + + 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: Initialize MFC + run: ./mfc.sh init + + - name: Run Sod shock tube L1 convergence tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_sod.py \ + --resolutions 64 128 256 512 \ + --num-ranks 4 + + convergence-2d-axisym: + name: "2D Axisymmetric WENO5 Convergence" + runs-on: ubuntu-latest + timeout-minutes: 60 + + 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: Initialize MFC + run: ./mfc.sh init + + - name: Run 2D axisymmetric convergence test + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_convergence_axisym.py \ + --resolutions 64 128 256 + + convergence-temporal: + name: "RK3 Temporal Order" + runs-on: ubuntu-latest + timeout-minutes: 60 + + 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: Initialize MFC + run: ./mfc.sh init + + - name: Run temporal order tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_temporal_order.py \ + --num-ranks 4 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..671355780f --- /dev/null +++ b/examples/1D_euler_convergence/case.py @@ -0,0 +1,105 @@ +#!/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)") +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 = 1.0 +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.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_axisym_convergence/case.py b/examples/2D_axisym_convergence/case.py new file mode 100644 index 0000000000..8dff936d46 --- /dev/null +++ b/examples/2D_axisym_convergence/case.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +""" +2D axisymmetric convergence case: density sine wave in x (axial/z), cyl_coord=T. + +In MFC 2D cylindrical: x=axial(z), y=radial(r). +rho = 1 + 0.2*sin(2*pi*x), u_x=1, p=1, u_y=0. +Exact solution at time T: rho(x,T) = 1 + 0.2*sin(2*pi*(x-T)). +Domain [0,5] = 5 full periods of the sine IC; periodic x-BCs make the exact +solution valid everywhere (no upstream contamination from ghost cells). +""" + +import argparse +import json +import math +import sys + +parser = argparse.ArgumentParser() +parser.add_argument("--mfc", type=str, default=None) +parser.add_argument("-N", type=int, default=64, help="Grid cells in x (axial) direction") +parser.add_argument("--order", type=int, default=5) +parser.add_argument("--cfl", type=float, default=0.02) +parser.add_argument("--nr", type=int, default=6, help="Radial cells (fixed)") +parser.add_argument("--Tfinal", type=float, default=0.1, help="Final simulation time") +args = parser.parse_args() + +N = args.N # axial cells (x-direction, refined) +Nr = args.nr # radial cells (y-direction, fixed) +# Domain [0,5] = 5 periods of sin(2*pi*x); periodic x-BCs make it exact everywhere +Lx = 5.0 +Lr = 0.5 # radial domain length + +dx = Lx / N +dt = args.cfl * dx +Nt = math.ceil(args.Tfinal / dt) +dt = args.Tfinal / Nt # adjust to hit Tfinal exactly + +case = { + # x=axial, y=radial + "m": N - 1, + "n": Nr - 1, + "p": 0, + "x_domain%beg": 0.0, + "x_domain%end": Lx, + "y_domain%beg": 0.0, + "y_domain%end": Lr, + "cyl_coord": "T", + "dt": dt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Nt, + "weno_order": args.order, + "weno_eps": 1e-16, + "mapped_weno": "F", + "wenoz": "F", + "teno": "F", + "mp_weno": "F", + "time_stepper": 3, + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "model_eqns": 2, + "num_fluids": 1, + "fluid_pp(1)%gamma": 1.4, + "fluid_pp(1)%pi_inf": 0.0, + # x: periodic (domain [0,5] = 5 full periods of the sine IC; wave advects cleanly) + # y: symmetry axis at r=0, extrapolation at outer r (same as CI tests) + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -2, + "bc_y%end": -3, + "num_patches": 1, + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": Lx / 2, + "patch_icpp(1)%y_centroid": Lr / 2, + "patch_icpp(1)%length_x": Lx, + "patch_icpp(1)%length_y": Lr, + "patch_icpp(1)%vel(1)": 1.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": "1.0 + 0.2*sin(2.0*acos(-1.0_wp)*x)", + "patch_icpp(1)%alpha(1)": 1.0, +} + +if args.mfc is not None: + print(json.dumps(case)) diff --git a/examples/2D_isentropicvortex_convergence/case.py b/examples/2D_isentropicvortex_convergence/case.py new file mode 100644 index 0000000000..865908d98e --- /dev/null +++ b/examples/2D_isentropicvortex_convergence/case.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python3 +# 2D isentropic vortex convergence case (weak-vortex, small-epsilon formulation). +# +# Uses hcid=283: 3-pt Gauss-Legendre cell averages of conserved variables as IC. +# Vortex strength eps=0.01 (weak vortex) so the primitive→conserved covariance +# error O(eps^3 h^2) stays well below the WENO5 scheme error O(eps^2 h^5) at +# resolutions N=16..64. With periodic BCs and a stationary vortex the comparison +# L2(rho(T) - rho(0)) isolates the scheme's accumulated spatial truncation error. +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="2D isentropic vortex 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 (default: 5)") +parser.add_argument("--muscl", action="store_true", help="Use MUSCL 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("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)") +parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4; use 0.005 for WENO7/TENO7 so temporal error is negligible)") +args = parser.parse_args() + +gamma = 1.4 +eps_vortex = 0.01 # vortex strength: small enough that prim->cons floor is negligible +N = args.N +m = N - 1 +dx = 10.0 / N + +# Max wave speed: c_sound at ambient + max rotational velocity (at r~0.7 for exp(1-r^2)) +c_max = math.sqrt(gamma) + eps_vortex / (2.0 * math.pi) +dt = args.cfl * dx / c_max +T_end = 2.0 +Nt = max(4, math.ceil(T_end / dt)) +dt = T_end / Nt # adjust 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-40, + "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": -5.0, + "x_domain%end": 5.0, + "y_domain%beg": -5.0, + "y_domain%end": 5.0, + "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": 3, + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -4, + "bc_x%end": -4, + "bc_y%beg": -4, + "bc_y%end": -4, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%length_x": 10.0, + "patch_icpp(1)%length_y": 10.0, + "patch_icpp(1)%hcid": 283, + "patch_icpp(1)%epsilon": eps_vortex, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": 1.0, + "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 dd151210c6..214f575857 100644 --- a/src/simulation/m_muscl.fpp +++ b/src/simulation/m_muscl.fpp @@ -170,7 +170,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/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 33661b54d7..78daecd625 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -772,7 +772,7 @@ contains ! Recalculating the radial momentum geometric source flux flux_gsrc_rs${XYZ}$_vf(j, k, l, eqn_idx%cont%end + 2) = flux_rs${XYZ}$_vf(j, k, l, & & eqn_idx%cont%end + 2) - (s_M*pres_R - s_P*pres_L)/(s_M - s_P) - ! Geometrical source of the void fraction(s) is zero + ! Cylindrical correction for HLL/LF numerical diffusion of void fraction(s) $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%adv%beg, eqn_idx%adv%end flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i) @@ -1382,7 +1382,7 @@ contains ! Recalculating the radial momentum geometric source flux flux_gsrc_rs${XYZ}$_vf(j, k, l, eqn_idx%cont%end + 2) = flux_rs${XYZ}$_vf(j, k, l, & & eqn_idx%cont%end + 2) - (s_M*pres_R - s_P*pres_L)/(s_M - s_P) - ! Geometrical source of the void fraction(s) is zero + ! Cylindrical correction for HLL/LF numerical diffusion of void fraction(s) $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%adv%beg, eqn_idx%adv%end flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index ce0ae48c4b..91b083b61e 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -548,6 +548,12 @@ contains if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) + ! Single-fluid cyl_coord: alpha is trivially 1 but drifts due to varying HLLC contact velocity across faces. Reset to + ! prevent pressure NaN. + if (num_fluids == 1 .and. cyl_coord .and. .not. bubbles_euler) then + call s_reset_single_fluid_alpha(q_cons_ts(1)%vf) + end if + if (ib) then ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms if (moving_immersed_boundary_flag) then @@ -703,6 +709,25 @@ contains end subroutine s_apply_bodyforces + !> Reset the single-fluid volume-fraction field to 1, preventing per-stage drift caused by varying contact-wave speed across + !! cylindrical faces. + subroutine s_reset_single_fluid_alpha(q_cons_vf) + + type(scalar_field), dimension(1:sys_size), intent(inout) :: q_cons_vf + integer :: j, k, l + + $:GPU_PARALLEL_LOOP(collapse=3) + do l = 0, p + do k = 0, n + do j = 0, m + q_cons_vf(eqn_idx%adv%beg)%sf(j, k, l) = 1._wp + end do + end do + end do + $:END_GPU_PARALLEL_LOOP() + + end subroutine s_reset_single_fluid_alpha + !> Update immersed boundary positions and velocities at the current Runge-Kutta stage subroutine s_propagate_immersed_boundaries(s) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 687431f54b..73b6b4ce66 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", @@ -803,7 +803,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 66166f43c8..a25bc33c65 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/cases.py b/toolchain/mfc/test/cases.py index 097aac4e69..cc89c8c72c 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -1569,6 +1569,10 @@ def foreach_example(): "1D_titarevtorro_analytical", "2D_acoustic_pulse_analytical", "2D_isentropicvortex_analytical", + "2D_isentropicvortex_convergence", + "1D_euler_convergence", + "1D_advection_convergence", + "1D_sod_convergence", "2D_zero_circ_vortex_analytical", "3D_TaylorGreenVortex_analytical", "3D_IGR_TaylorGreenVortex_nvidia", @@ -1580,6 +1584,7 @@ def foreach_example(): "2D_ibm_stl_MFCCharacter", "1D_qbmm", "2D_Thermal_Flatplate", # formatted I/O field overflow on gfortran 12 + "2D_axisym_convergence", # WENO5 stencil requires n+1 >= 25; case uses nr=6 (run via run_convergence_axisym.py) ] if path in casesToSkip: continue diff --git a/toolchain/mfc/test/run_convergence.py b/toolchain/mfc/test/run_convergence.py new file mode 100644 index 0000000000..101ed6644b --- /dev/null +++ b/toolchain/mfc/test/run_convergence.py @@ -0,0 +1,254 @@ +#!/usr/bin/env python3 +""" +Convergence-rate verification for MFC's 2D isentropic vortex problem. + +Uses hcid=283: 3-pt Gauss-Legendre cell averages of conserved variables as IC. +The vortex strength eps=0.01 (set in case.py) is chosen so that the dominant +error source is the WENO spatial truncation error O(eps^2 * h^p), not the +primitive-to-conserved covariance floor O(eps^3 * h^2). For h > eps^(1/3)=0.22 +(i.e., N < 46 per dimension), the p-th order scheme shows rate p. + +L2(rho(T) - rho(0)) measures accumulated scheme error; the comparison to rho(0) +(the numerical IC) eliminates IC discretisation error, isolating the scheme error. + +WENO7/TENO7 are NOT tested here. For the isentropic vortex, the IC +primitive→conserved covariance error is O(eps^3 * h^2). The WENO7 scheme +error is O(eps^2 * h^7). Scheme error dominates only when h > eps^(1/5); +with eps=0.01 that requires h > 0.40, i.e., N < 25. At N=64-128 the +covariance floor dominates and the measured rate is ~2, not 7. +WENO7/TENO7 7th-order convergence is verified by the 1D test (run_convergence_1d.py) +which uses a pure advection problem that avoids this nonlinear floor. + +Usage: + python toolchain/mfc/test/run_convergence.py [--resolutions 32 64 128] +""" + +import argparse +import json +import math +import os +import shutil +import struct +import subprocess +import sys +import tempfile + +import numpy as np + +CASE = "examples/2D_isentropicvortex_convergence/case.py" +MFC = "./mfc.sh" + +# (label, extra_args, expected_order, tolerance, min_N, max_N) +# With eps=0.01 and N=32..128 the prim->cons covariance error O(eps^3 h^2) is +# well below the scheme's spatial error O(eps^2 h^p), so each scheme shows its +# nominal rate. The tolerance is the allowable shortfall from the nominal order. +# +# WENO3: at N=32-128 the rate is ~2.0-2.2 (pre-asymptotic; approaches 3 at +# finer grids). Threshold 1.8. +# WENO7/TENO7 are omitted — see module docstring for why. +SCHEMES = [ + # WENO5/TENO5: need (N/2) >= num_stcls_min*recon_order = 25 cells/rank with 4 + # ranks in 2x2; min_N=64 ensures 32 cells/rank which satisfies this constraint. + ("WENO5", ["--order", "5"], 5, 1.0, 64, None), + ("WENO3", ["--order", "3"], 3, 1.2, 32, None), + ("WENO1", ["--order", "1"], 1, 0.4, 32, None), + ("MUSCL2", ["--muscl"], 2, 0.5, 32, None), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6"], 5, 1.0, 64, None), +] + + +def read_cons_var(run_dir: str, step: int, var_idx: int, N: int, num_ranks: int = 1) -> np.ndarray: + """Read q_cons_vf{var_idx} from all MPI ranks and return as a flat array. + + Spatial ordering is not preserved across ranks but that is fine for L2 + norm and sum computations, which are invariant to permutation of elements. + """ + 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] + data = np.frombuffer(f.read(rec_len), dtype=np.float64) + f.read(4) + chunks.append(data.copy()) + combined = np.concatenate(chunks) + if combined.size != N * N: + raise ValueError(f"Expected {N * N} values across {num_ranks} ranks, got {combined.size}") + return combined + + +# 2D single-fluid Euler (model_eqns=2, num_fluids=1): vf1=ρ, vf2=ρu, vf3=ρv, vf4=E +# Momentum is excluded: the isentropic vortex has zero net linear momentum, making +# the relative error formula ill-conditioned (denominator ≈ 0). Density and energy +# have large nonzero integrals and are the meaningful conserved quantities to verify. +CONS_VARS_2D = [("density", 1), ("energy", 4)] +CONS_TOL = 1e-10 + + +def conservation_errors(run_dir: str, Nt: int, N: int, cell_vol: float, var_list: list, num_ranks: int) -> dict: + """Return relative conservation error |Σq(T) - Σq(0)| / |Σq(0)| for each variable.""" + errs = {} + for name, idx in var_list: + q0 = read_cons_var(run_dir, 0, idx, N, num_ranks) + qT = read_cons_var(run_dir, Nt, idx, N, num_ranks) + s0 = float(np.sum(q0)) * cell_vol + sT = float(np.sum(qT)) * cell_vol + errs[name] = abs(sT - s0) / (abs(s0) + 1e-300) + return errs + + +def l2_error(rho_final: np.ndarray, rho_init: np.ndarray, dx: float) -> float: + """L2 error: sqrt(sum((f-g)^2 * dx^2)).""" + diff = rho_final - rho_init + return float(np.sqrt(np.sum(diff**2) * dx**2)) + + +def convergence_rate(errors: list, resolutions: list) -> float: + """Least-squares slope of log(error) vs log(dx), dx = 10/N.""" + log_dx = np.log(10.0 / np.array(resolutions, dtype=float)) + log_err = np.log(np.array(errors, dtype=float)) + rate, _ = np.polyfit(log_dx, log_err, 1) + return float(rate) + + +def run_case(tmpdir: str, N: int, extra_args: list, num_ranks: int = 1): + """Run the vortex case at resolution N. Returns (Nt, run_dir).""" + # Query case parameters to find t_step_stop + result = subprocess.run( + [sys.executable, CASE, "--mfc", "{}", "-N", str(N)] + extra_args, + capture_output=True, + text=True, + check=False, + ) + if result.returncode != 0: + raise RuntimeError(f"case.py failed:\n{result.stderr}") + cfg = json.loads(result.stdout) + Nt = int(cfg["t_step_stop"]) + + # Run only pre_process and simulation (post_process not needed for p_all) + cmd = [ + MFC, + "run", + CASE, + "-t", + "pre_process", + "simulation", + "-n", + str(num_ranks), + "--", + "-N", + str(N), + ] + extra_args + result = subprocess.run(cmd, capture_output=True, text=True, cwd=os.getcwd(), check=False) + if result.returncode != 0: + print(result.stdout[-2000:]) + print(result.stderr) + raise RuntimeError(f"./mfc.sh run failed for N={N}") + + # Copy p_all to temp dir, then clean the case directory for next run + case_dir = os.path.dirname(CASE) + src = os.path.join(case_dir, "p_all") + dst = os.path.join(tmpdir, f"N{N}", "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 Nt, os.path.join(tmpdir, f"N{N}") + + +def test_scheme(label, extra_args, expected_order, tol, resolutions, min_N=None, max_N=None, num_ranks=1): + if min_N is not None: + resolutions = [N for N in resolutions if N >= min_N] + if max_N is not None: + resolutions = [N for N in resolutions if N <= max_N] + print(f"\n{'=' * 60}") + print(f" {label} (need rate >= {expected_order - tol:.1f})") + print(f"{'=' * 60}") + + errors = [] + nts = [] + all_cons_errs = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in resolutions: + dx = 10.0 / N + Nt, run_dir = run_case(tmpdir, N, extra_args, num_ranks) + nts.append(Nt) + rho0 = read_cons_var(run_dir, 0, 1, N, num_ranks) + rhoT = read_cons_var(run_dir, Nt, 1, N, num_ranks) + err = l2_error(rhoT, rho0, dx) + errors.append(err) + all_cons_errs.append(conservation_errors(run_dir, Nt, N, dx**2, CONS_VARS_2D, num_ranks)) + + # Compute pairwise rates + rates = [None] + for i in range(1, len(resolutions)): + log_dx0 = math.log(10.0 / resolutions[i - 1]) + log_dx1 = math.log(10.0 / resolutions[i]) + rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (log_dx1 - log_dx0)) + + print(f" {'N':>6} {'Nt':>5} {'dx':>10} {'L2 error':>14} {'rate':>8}") + print(f" {'-' * 6} {'-' * 5} {'-' * 10} {'-' * 14} {'-' * 8}") + for i, N in enumerate(resolutions): + dx = 10.0 / N + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {N:>6} {nts[i]:>5} {dx:>10.5f} {errors[i]:>14.6e} {r_str}") + + if len(resolutions) > 1: + overall = convergence_rate(errors, resolutions) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + rate_passed = overall >= expected_order - tol + else: + print("\n (need >= 2 resolutions to compute rate)") + rate_passed = True # can't fail with a single resolution + + print(f"\n Conservation (need rel. error < {CONS_TOL:.0e}):") + cons_passed = True + for name, _ in CONS_VARS_2D: + max_err = max(ce[name] for ce in all_cons_errs) + ok = max_err < CONS_TOL + print(f" {name:<14}: max = {max_err:.2e} {'OK' if ok else 'FAIL'}") + if not ok: + cons_passed = False + + passed = rate_passed and cons_passed + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC convergence-rate verification") + parser.add_argument("--resolutions", type=int, nargs="+", default=[32, 64, 128], help="Grid resolutions (default: 32 64 128; N<32 unsupported for WENO5)") + parser.add_argument("--schemes", nargs="+", default=["WENO5", "WENO3", "WENO1", "MUSCL2", "TENO5"], help="Schemes to test (default: all)") + parser.add_argument("--num-ranks", type=int, default=1, help="MPI ranks per simulation (default: 1)") + args = parser.parse_args() + + results = {} + for label, extra_args, expected_order, tol, min_N, max_N in SCHEMES: + if label not in args.schemes: + continue + try: + passed = test_scheme(label, extra_args, expected_order, tol, args.resolutions, min_N, max_N, args.num_ranks) + except Exception as e: + import traceback + + print(f" ERROR: {e}") + traceback.print_exc() + passed = False + results[label] = passed + + print(f"\n{'=' * 60}") + print(" Summary") + print(f"{'=' * 60}") + all_pass = True + for label, passed in results.items(): + print(f" {label:<12} {'PASS' if passed else 'FAIL'}") + if not passed: + all_pass = False + + sys.exit(0 if all_pass else 1) + + +if __name__ == "__main__": + main() diff --git a/toolchain/mfc/test/run_convergence_1d.py b/toolchain/mfc/test/run_convergence_1d.py new file mode 100644 index 0000000000..679819328f --- /dev/null +++ b/toolchain/mfc/test/run_convergence_1d.py @@ -0,0 +1,269 @@ +#!/usr/bin/env python3 +""" +Convergence-rate verification for MFC's 1D single-fluid Euler equations. + +Single fluid with a density sine wave: rho = 1 + 0.2*sin(2*pi*x), u=1, p=1. +After exactly one period (T=1, u=1, L=1), the exact solution equals the IC. +L2(rho(T) - rho(0)) measures the accumulated scheme spatial truncation error. +No non-conservative alpha equation — clean benchmark for all schemes. + +WENO5/TENO5 use CFL=0.02: RK3 temporal error O(dt^3) is then negligible +relative to the O(h^5) spatial error at N=128-512. + +WENO7/TENO7 use CFL=0.005: at CFL=0.02 the RK3 temporal error (~3.4e-12 at +N=128) is comparable to the spatial error (~4.4e-12), giving a spurious rate +of ~3.7. With CFL=0.005 the temporal error drops by (0.005/0.02)^3 = 1/64 +to ~5.3e-14, well below spatial, and the measured rate approaches 7. +N is capped at 256 — the machine-precision floor is reached near N=512. + +WENO3-JS degrades to 2nd order at smooth extrema (Henrick et al. 2005). +The expected rate for WENO3 here is therefore 2, not 3; the 2D isentropic +vortex test (run_convergence.py) verifies WENO3 rate 3. + +MUSCL2 uses muscl_lim=0 (unlimited central-difference) by default. TVD +limiters clip slopes to zero at smooth extrema and stall at 1st order on the +sine wave; the unlimited limiter preserves 2nd-order convergence everywhere. + +Usage: + python toolchain/mfc/test/run_convergence_1d.py [--resolutions 128 256 512 1024] +""" + +import argparse +import json +import math +import os +import shutil +import struct +import subprocess +import sys +import tempfile + +import numpy as np + +CASE = "examples/1D_euler_convergence/case.py" +MFC = "./mfc.sh" + +# (label, extra_args, expected_order, tolerance, min_N, max_N) +# CFL is baked into each scheme's extra_args so that WENO7/TENO7 can use a +# smaller CFL independently of all other schemes. +# +# Per-scheme resolution bounds let each scheme run over the range where its +# asymptotic order is cleanly visible: +# WENO5 : cap at N=512 — double-precision floor kills the rate at N=1024 +# (error ~2.6e-12, rate collapses to 0.69); [128,512] gives 4.99. +# WENO3 : start at N=256 — skips the coarsest pre-asymptotic points; +# WENO3-JS degrades to 2nd order at smooth extrema (Henrick 2005), +# asymptote confirmed 1.99 at N=4096; [256,1024] gives ~1.87. +# WENO1 : full range [128,1024]; rate 0.97. +# MUSCL2 : full range [128,1024]; unlimited slope, rate exactly 2.00. +# TENO5 : same range as WENO5; CT=1e-6; rate matches WENO5 on smooth problems. +# WENO7 : CFL=0.005, range [64,128] — at N=256 the spatial error (~1.7e-14) +# falls below the round-off accumulation floor (~2.5e-12 for ~28M +# cell-steps), so only N=64 and N=128 give a clean rate ≥6.5. +# TENO7 : same range and CFL as WENO7; CT=1e-9. +SCHEMES = [ + ("WENO5", ["--order", "5", "--cfl", "0.02"], 5, 0.2, 128, 512), + ("WENO3", ["--order", "3", "--cfl", "0.02"], 2, 0.2, 256, None), + ("WENO1", ["--order", "1", "--cfl", "0.02"], 1, 0.05, 128, None), + ("MUSCL2", ["--muscl", "--cfl", "0.02"], 2, 0.1, 128, None), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 5, 0.2, 128, 512), + ("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), +] + + +def read_cons_var(run_dir: str, step: int, var_idx: int, num_ranks: int = 1, expected_size: int = None) -> np.ndarray: + """Read q_cons_vf{var_idx} from all MPI ranks and concatenate into one 1D array.""" + 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] + data = np.frombuffer(f.read(rec_len), dtype=np.float64) + f.read(4) + chunks.append(data.copy()) + combined = np.concatenate(chunks) + if expected_size is not None and combined.size != expected_size: + raise ValueError(f"Expected {expected_size} values across {num_ranks} ranks, got {combined.size}") + return combined + + +# 1D single-fluid Euler (model_eqns=2, num_fluids=1): vf1=ρ, vf2=ρu, vf3=E +CONS_VARS_1D = [("density", 1), ("x-momentum", 2), ("energy", 3)] +CONS_TOL = 1e-10 + + +def conservation_errors(run_dir: str, Nt: int, cell_vol: float, var_list: list, num_ranks: int, expected_size: int = None) -> dict: + """Return relative conservation error |Σq(T) - Σq(0)| / |Σq(0)| for each variable.""" + errs = {} + for name, idx in var_list: + q0 = read_cons_var(run_dir, 0, idx, num_ranks, expected_size) + qT = read_cons_var(run_dir, Nt, idx, num_ranks, expected_size) + s0 = float(np.sum(q0)) * cell_vol + sT = float(np.sum(qT)) * cell_vol + errs[name] = abs(sT - s0) / (abs(s0) + 1e-300) + return errs + + +def l2_error(a: np.ndarray, b: np.ndarray, dx: float) -> float: + return float(np.sqrt(np.sum((a - b) ** 2) * dx)) + + +def convergence_rate(errors: list, resolutions: list) -> float: + log_dx = np.log(1.0 / np.array(resolutions, dtype=float)) + log_err = np.log(np.array(errors, dtype=float)) + rate, _ = np.polyfit(log_dx, log_err, 1) + return float(rate) + + +def run_case(tmpdir: str, N: int, extra_args: list, num_ranks: int = 1): + """Run the 1D advection case at resolution N. Returns (Nt, run_dir).""" + result = subprocess.run( + [sys.executable, CASE, "--mfc", "{}", "-N", str(N)] + extra_args, + capture_output=True, + text=True, + check=False, + ) + if result.returncode != 0: + raise RuntimeError(f"case.py failed:\n{result.stderr}") + cfg = json.loads(result.stdout) + Nt = int(cfg["t_step_stop"]) + + cmd = [ + MFC, + "run", + CASE, + "-t", + "pre_process", + "simulation", + "-n", + str(num_ranks), + "--", + "-N", + str(N), + ] + extra_args + result = subprocess.run(cmd, capture_output=True, text=True, cwd=os.getcwd(), check=False) + if result.returncode != 0: + print(result.stdout[-3000:]) + print(result.stderr) + raise RuntimeError(f"./mfc.sh run failed for N={N}") + + case_dir = os.path.dirname(CASE) + src = os.path.join(case_dir, "p_all") + dst = os.path.join(tmpdir, f"N{N}", "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 Nt, os.path.join(tmpdir, f"N{N}") + + +def test_scheme(label, extra_args, expected_order, tol, resolutions, min_N=None, max_N=None, num_ranks=1): + if min_N is not None: + resolutions = [N for N in resolutions if N >= min_N] + if max_N is not None: + resolutions = [N for N in resolutions if N <= max_N] + print(f"\n{'=' * 60}") + print(f" {label} (need rate >= {expected_order - tol:.1f})") + print(f"{'=' * 60}") + + errors = [] + nts = [] + all_cons_errs = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in resolutions: + dx = 1.0 / N + Nt, run_dir = run_case(tmpdir, N, extra_args, num_ranks) + nts.append(Nt) + vf0 = read_cons_var(run_dir, 0, 1, num_ranks, expected_size=N) + vfT = read_cons_var(run_dir, Nt, 1, num_ranks, expected_size=N) + err = l2_error(vfT, vf0, dx) + errors.append(err) + all_cons_errs.append(conservation_errors(run_dir, Nt, dx, CONS_VARS_1D, num_ranks, expected_size=N)) + print(f" N={N}: Nt={Nt}, |vf0|={len(vf0)}, err={err:.4e}") + + rates = [None] + for i in range(1, len(resolutions)): + log_dx0 = math.log(1.0 / resolutions[i - 1]) + log_dx1 = math.log(1.0 / resolutions[i]) + rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (log_dx1 - log_dx0)) + + print() + print(f" {'N':>6} {'Nt':>6} {'dx':>10} {'L2 error':>14} {'rate':>8}") + print(f" {'-' * 6} {'-' * 6} {'-' * 10} {'-' * 14} {'-' * 8}") + for i, N in enumerate(resolutions): + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {N:>6} {nts[i]:>6} {1.0 / N:>10.6f} {errors[i]:>14.6e} {r_str}") + + if len(resolutions) > 1: + overall = convergence_rate(errors, resolutions) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + rate_passed = overall >= expected_order - tol + else: + rate_passed = True + + print(f"\n Conservation (need rel. error < {CONS_TOL:.0e}):") + cons_passed = True + for name, _ in CONS_VARS_1D: + max_err = max(ce[name] for ce in all_cons_errs) + ok = max_err < CONS_TOL + print(f" {name:<14}: max = {max_err:.2e} {'OK' if ok else 'FAIL'}") + if not ok: + cons_passed = False + + passed = rate_passed and cons_passed + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC 1D advection convergence-rate verification") + parser.add_argument( + "--resolutions", + type=int, + nargs="+", + default=[64, 128, 256, 512, 1024], + help="Grid resolutions to test (default: 64 128 256 512 1024)", + ) + parser.add_argument( + "--schemes", + nargs="+", + default=["WENO5", "WENO3", "WENO1", "MUSCL2", "TENO5", "WENO7", "TENO7"], + help="Schemes to test", + ) + parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter (0=unlimited 1=minmod ...; default: 0)") + parser.add_argument("--num-ranks", type=int, default=1, help="MPI ranks per simulation (default: 1)") + args = parser.parse_args() + + muscl_extra = ["--muscl-lim", str(args.muscl_lim)] + + results = {} + for label, extra_args, expected_order, tol, min_N, max_N in SCHEMES: + if label not in args.schemes: + continue + try: + passed = test_scheme(label, extra_args + muscl_extra, expected_order, tol, args.resolutions, min_N, max_N, args.num_ranks) + except Exception as e: + import traceback + + print(f" ERROR: {e}") + traceback.print_exc() + passed = False + results[label] = passed + + print(f"\n{'=' * 60}") + print(" Summary") + print(f"{'=' * 60}") + all_pass = True + for label, passed in results.items(): + print(f" {label:<12} {'PASS' if passed else 'FAIL'}") + if not passed: + all_pass = False + + sys.exit(0 if all_pass else 1) + + +if __name__ == "__main__": + main() diff --git a/toolchain/mfc/test/run_convergence_axisym.py b/toolchain/mfc/test/run_convergence_axisym.py new file mode 100644 index 0000000000..a4bf6d4f41 --- /dev/null +++ b/toolchain/mfc/test/run_convergence_axisym.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python3 +""" +Convergence-rate verification for WENO5 on a 2D axisymmetric (cyl_coord=T) grid. + +Density sine wave in z: rho = 1 + 0.2*sin(2*pi*z), u_z=1, p=1, u_r=0. +Exact solution at time T: rho_exact(z,T) = 1 + 0.2*sin(2*pi*(z-T)). +Nr is held fixed; Nz is refined. +L2 error is computed by averaging density over r then comparing to exact solution. +""" + +import argparse +import json +import math +import os +import shutil +import struct +import subprocess +import sys +import tempfile + +import numpy as np + +CASE = "examples/2D_axisym_convergence/case.py" +MFC = "./mfc.sh" +NR = 32 # fixed radial cells (n=31; WENO5 needs n+1 >= num_stcls_min*5 = 25) +TFINAL = 0.1 # short final time avoids long-time cylindrical instability + + +def read_field(run_dir: str, step: int, var_idx: int, nz: int, nr: int) -> np.ndarray: + """Read 2D field. In MFC: x=axial(Nz), y=radial(Nr). Fortran col-major -> shape (Nz, Nr).""" + path = os.path.join(run_dir, "p_all", "p0", str(step), f"q_cons_vf{var_idx}.dat") + with open(path, "rb") as f: + rec_len = struct.unpack("i", f.read(4))[0] + data = np.frombuffer(f.read(rec_len), dtype=np.float64).copy() + f.read(4) + if data.size != nr * nz: + raise ValueError(f"Expected {nr * nz} values, got {data.size}") + # Fortran stores (x, y) = (axial, radial) in column-major: first index varies fastest + return data.reshape(nz, nr, order="F") + + +def run_case(tmpdir: str, nz: int, extra_args: list) -> tuple: + """Run the 2D axisym case at resolution nz. Returns (Nt, run_dir).""" + result = subprocess.run( + [sys.executable, CASE, "--mfc", "{}", "-N", str(nz), "--nr", str(NR), "--Tfinal", str(TFINAL)] + extra_args, + capture_output=True, + text=True, + check=False, + ) + if result.returncode != 0: + raise RuntimeError(f"case.py failed:\n{result.stderr}") + cfg = json.loads(result.stdout) + Nt = int(cfg["t_step_stop"]) + dt = float(cfg["dt"]) + + cmd = [MFC, "run", CASE, "-t", "pre_process", "simulation", "-n", "1", "--", "-N", str(nz), "--nr", str(NR), "--Tfinal", str(TFINAL)] + extra_args + result = subprocess.run(cmd, capture_output=True, text=True, cwd=os.getcwd(), check=False) + if result.returncode != 0: + print(result.stdout[-3000:]) + print(result.stderr) + raise RuntimeError(f"./mfc.sh run failed for Nz={nz}") + + case_dir = os.path.dirname(CASE) + src = os.path.join(case_dir, "p_all") + dst = os.path.join(tmpdir, f"N{nz}", "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 Nt, dt, os.path.join(tmpdir, f"N{nz}") + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--resolutions", type=int, nargs="+", default=[64, 128, 256]) + parser.add_argument("--cfl", type=float, default=0.02) + parser.add_argument("--order", type=int, default=5) + args = parser.parse_args() + + extra = ["--cfl", str(args.cfl), "--order", str(args.order)] + + print(f"\n{'=' * 60}") + print(f" WENO{args.order} on 2D axisymmetric (cyl_coord=T) grid") + print(f" Sine wave in z, Nr={NR} fixed, Nz refined, T={TFINAL}") + print(f"{'=' * 60}") + + errors, nts = [], [] + with tempfile.TemporaryDirectory() as tmpdir: + for nz in args.resolutions: + Lx = 5.0 # must match case.py + dz = Lx / nz + Nt, dt, run_dir = run_case(tmpdir, nz, extra) + nts.append(Nt) + T_actual = Nt * dt # actual final time + + # Read final density field (Nz x Nr), average over r + rhoT = read_field(run_dir, Nt, 1, nz, NR) + rhoT_z = rhoT.mean(axis=1) # shape (nz,) + + # Exact solution: rho(z, T) = 1 + 0.2*sin(2*pi*(z - T)) + x_cc = (np.arange(nz) + 0.5) * dz + rho_exact = 1.0 + 0.2 * np.sin(2.0 * np.pi * (x_cc - T_actual)) + + err = float(np.sqrt(np.sum((rhoT_z - rho_exact) ** 2) * dz)) + errors.append(err) + print(f" Nz={nz:4d}: Nt={Nt}, err={err:.4e}") + + Lx = 5.0 + rates = [None] + for i in range(1, len(args.resolutions)): + nz0, nz1 = args.resolutions[i - 1], args.resolutions[i] + rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (math.log(Lx / nz1) - math.log(Lx / nz0))) + + Lx = 5.0 + print(f"\n {'Nz':>6} {'Nt':>6} {'dz':>10} {'L2 error':>14} {'rate':>8}") + print(f" {'-' * 6} {'-' * 6} {'-' * 10} {'-' * 14} {'-' * 8}") + for i, nz in enumerate(args.resolutions): + r = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {nz:>6} {nts[i]:>6} {Lx / nz:>10.6f} {errors[i]:>14.6e} {r}") + + if len(args.resolutions) > 1: + log_dz = np.log(Lx / np.array(args.resolutions, dtype=float)) + log_err = np.log(np.array(errors, dtype=float)) + rate, _ = np.polyfit(log_dz, log_err, 1) + expected = args.order - 0.2 + passed = rate >= expected + print(f"\n Fitted rate: {rate:.2f} (need >= {expected:.1f})") + print(f" {'PASS' if passed else 'FAIL'}") + sys.exit(0 if passed else 1) + + +if __name__ == "__main__": + main() diff --git a/toolchain/mfc/test/run_sod.py b/toolchain/mfc/test/run_sod.py new file mode 100644 index 0000000000..774900422d --- /dev/null +++ b/toolchain/mfc/test/run_sod.py @@ -0,0 +1,229 @@ +#!/usr/bin/env python3 +""" +L1 self-convergence study for the 1D Sod shock tube across all MFC schemes. + +Sod problem: rho_L=1, u_L=0, p_L=1; rho_R=0.125, u_R=0, p_R=0.1; T=0.2. +Contains a shock, contact discontinuity, and rarefaction fan. + +By Godunov's theorem, any conservative monotone scheme converges at 1st order +in L1 for problems with shocks. Higher-order schemes (WENO5, TENO7, ...) also +achieve L1 rate ~1 globally because the shock contributes an O(h) error that +dominates the smooth-region high-order accuracy. + +Self-convergence method: run at N and 2N, cell-average the finer solution to +the coarse grid, compute L1(rho_N - avg(rho_{2N})). No exact solution needed. +Ranks are read in rank order, which equals spatial order for 1D decomposition. + +Usage: + python toolchain/mfc/test/run_sod.py + python toolchain/mfc/test/run_sod.py --resolutions 64 128 256 512 --schemes WENO5 TENO5 +""" + +import argparse +import json +import math +import os +import shutil +import struct +import subprocess +import sys +import tempfile + +import numpy as np + +CASE = "examples/1D_sod_convergence/case.py" +MFC = "./mfc.sh" + +# (label, extra_args, expected_order, tolerance, min_N) +# All schemes achieve L1 rate ~1 on shock problems (Godunov's theorem). +SCHEMES = [ + # WENO1 (1st-order upwind): contact discontinuity smears over O(sqrt(h*T)) width, + # giving L1 contribution O(h^0.5), so fitted rate ~0.6-0.7 is physically expected. + ("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), + # SUPERBEE is over-compressive near contacts; at N=64 the rate is pre-asymptotic + # (~0.40), so min_N=128 skips the pre-asymptotic point and gives a reliable fit. + ("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), +] + + +def read_cons_var(run_dir: str, step: int, var_idx: int, num_ranks: int = 1, expected_size: int = None) -> np.ndarray: + """Read q_cons_vf{var_idx} from all ranks in rank order (= spatial order for 1D).""" + 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] + data = np.frombuffer(f.read(rec_len), dtype=np.float64) + f.read(4) + chunks.append(data.copy()) + combined = np.concatenate(chunks) + if expected_size is not None and combined.size != expected_size: + raise ValueError(f"Expected {expected_size} values across {num_ranks} ranks, got {combined.size}") + return combined + + +def l1_self_error(coarse: np.ndarray, fine: np.ndarray, dx_coarse: float) -> float: + """L1 diff between coarse solution and cell-averaged fine solution.""" + assert len(fine) == 2 * len(coarse), f"Expected 2:1 ratio, got {len(fine)}:{len(coarse)}" + fine_avg = (fine[0::2] + fine[1::2]) / 2.0 + return float(np.sum(np.abs(coarse - fine_avg)) * dx_coarse) + + +def run_case(tmpdir: str, N: int, extra_args: list, num_ranks: int = 1): + """Run the Sod case at resolution N. Returns (Nt, run_dir).""" + result = subprocess.run( + [sys.executable, CASE, "--mfc", "{}", "-N", str(N)] + extra_args, + capture_output=True, + text=True, + check=False, + ) + if result.returncode != 0: + raise RuntimeError(f"case.py failed:\n{result.stderr}") + cfg = json.loads(result.stdout) + Nt = int(cfg["t_step_stop"]) + + cmd = [ + MFC, + "run", + CASE, + "-t", + "pre_process", + "simulation", + "-n", + str(num_ranks), + "--", + "-N", + str(N), + ] + extra_args + result = subprocess.run(cmd, capture_output=True, text=True, cwd=os.getcwd(), check=False) + if result.returncode != 0: + print(result.stdout[-3000:]) + print(result.stderr) + raise RuntimeError(f"./mfc.sh run failed for N={N}") + + case_dir = os.path.dirname(CASE) + src = os.path.join(case_dir, "p_all") + dst = os.path.join(tmpdir, f"N{N}", "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 Nt, os.path.join(tmpdir, f"N{N}") + + +def test_scheme(label, extra_args, expected_order, tol, resolutions, min_N=None, num_ranks=1): + if min_N is not None: + resolutions = [N for N in resolutions if N >= min_N] + print(f"\n{'=' * 60}") + print(f" {label} (need L1 rate >= {expected_order - tol:.1f})") + print(f"{'=' * 60}") + + # Need consecutive pairs — resolutions must be factors of 2 apart + nts = [] + run_dirs = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in resolutions: + Nt, run_dir = run_case(tmpdir, N, extra_args, num_ranks) + nts.append(Nt) + run_dirs.append(run_dir) + + # Compute L1 self-errors: compare each N against 2N + errors = [] + error_resolutions = [] + for i in range(len(resolutions) - 1): + N_c = resolutions[i] + N_f = resolutions[i + 1] + if N_f != 2 * N_c: + continue # skip non-2x pairs + dx_c = 1.0 / N_c + rho_c = read_cons_var(run_dirs[i], nts[i], 1, num_ranks, expected_size=N_c) + rho_f = read_cons_var(run_dirs[i + 1], nts[i + 1], 1, num_ranks, expected_size=N_f) + err = l1_self_error(rho_c, rho_f, dx_c) + errors.append(err) + error_resolutions.append(N_c) + + rates = [None] + for i in range(1, len(errors)): + log_h0 = math.log(1.0 / error_resolutions[i - 1]) + log_h1 = math.log(1.0 / error_resolutions[i]) + rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (log_h1 - log_h0)) + + print(f" {'N':>6} {'Nt':>6} {'L1 self-err':>14} {'rate':>8}") + print(f" {'-' * 6} {'-' * 6} {'-' * 14} {'-' * 8}") + for i, N in enumerate(error_resolutions): + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {N:>6} {nts[i]:>6} {errors[i]:>14.6e} {r_str}") + + if len(errors) >= 2: + log_h = np.log(1.0 / np.array(error_resolutions, dtype=float)) + log_e = np.log(np.array(errors, dtype=float)) + overall, _ = np.polyfit(log_h, log_e, 1) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + passed = overall >= expected_order - tol + elif len(errors) == 1: + print(f"\n Single pair rate: {rates[-1]:.2f} (need >= {expected_order - tol:.1f})") + passed = rates[-1] >= expected_order - tol + else: + print("\n ERROR: need >= 2 consecutive 2x-apart resolutions to compute a rate") + passed = False + + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC Sod shock tube L1 convergence") + parser.add_argument( + "--resolutions", + type=int, + nargs="+", + default=[128, 256, 512, 1024], + help="Grid resolutions — must be consecutive factors of 2 (default: 128 256 512 1024)", + ) + parser.add_argument( + "--schemes", + nargs="+", + default=[s[0] for s in SCHEMES], # label is always first element + help="Schemes to test (default: all)", + ) + parser.add_argument("--num-ranks", type=int, default=1, help="MPI ranks per simulation (default: 1)") + args = parser.parse_args() + + results = {} + for label, extra_args, expected_order, tol, min_N in SCHEMES: + if label not in args.schemes: + continue + try: + passed = test_scheme(label, extra_args, expected_order, tol, args.resolutions, min_N, args.num_ranks) + except Exception as e: + import traceback + + print(f" ERROR: {e}") + traceback.print_exc() + passed = False + results[label] = passed + + print(f"\n{'=' * 60}") + print(" Summary") + print(f"{'=' * 60}") + all_pass = True + for label, passed in results.items(): + print(f" {label:<18} {'PASS' if passed else 'FAIL'}") + if not passed: + all_pass = False + + sys.exit(0 if all_pass else 1) + + +if __name__ == "__main__": + main() diff --git a/toolchain/mfc/test/run_temporal_order.py b/toolchain/mfc/test/run_temporal_order.py new file mode 100644 index 0000000000..0df4aa97d7 --- /dev/null +++ b/toolchain/mfc/test/run_temporal_order.py @@ -0,0 +1,247 @@ +#!/usr/bin/env python3 +""" +Time-integration order verification for MFC's RK1, RK2, and RK3 time steppers. + +Uses the 1D single-fluid Euler advection problem (rho = 1 + 0.2*sin(2*pi*x), +u=1, p=1, L=1, T=1) with a fine spatial grid (N=512, WENO5) so the spatial +error (~4e-12) is negligible compared to the temporal error at the CFLs tested. + +L2(rho(T) - rho(0)) measures total accumulated error. By fixing N and varying +CFL (and hence dt), the spatial contribution is constant and the measured rate +reflects the time integration order. + +CFL ranges are chosen to be within each stepper's stability region and keep +temporal errors well above the ~4e-12 spatial floor: + RK1 (Euler, 1st order): CFL=[0.10, 0.05] — stable limit ~0.1 with WENO5+LF + (nearly-imaginary eigenvalues constrain Euler more than TVD RK); + error ~2.5e-4 and ~1.2e-4 (rate ≈ 1.0) + RK2 (TVD Heun, 2nd order): CFL=[0.50, 0.25]; + error ~1.2e-6 and ~2.9e-7 (rate ≈ 2.0) + RK3 (TVD Shu-Osher, 3rd order): CFL=[0.50, 0.25]; + error ~8.3e-10 and ~1.1e-10 (rate ≈ 3.0) + +Usage: + python toolchain/mfc/test/run_temporal_order.py + python toolchain/mfc/test/run_temporal_order.py --schemes RK3/WENO5 --cfls 0.5 0.25 0.125 +""" + +import argparse +import json +import math +import os +import shutil +import struct +import subprocess +import sys +import tempfile + +import numpy as np + +CASE = "examples/1D_euler_convergence/case.py" +MFC = "./mfc.sh" + +# (label, extra_args, expected_order, tolerance, cfls) +# N=512 is fixed; WENO5 keeps spatial error ~4e-12 (negligible at CFL>=0.25). +# RK1/RK2 temporal errors at CFL=0.5 are ~2e-4 and ~2e-7, both >> spatial floor. +SCHEMES = [ + # RK1 (Forward Euler): nearly-imaginary WENO5+LF eigenvalues constrain stability + # to CFL < ~0.1; use [0.10, 0.05] which are provably stable and temporal-dominated + ("RK1/WENO5", ["--order", "5", "--time-stepper", "1"], 1, 0.1, [0.10, 0.05]), + # RK2/RK3 (TVD): stable to CFL~1; use [0.50, 0.25] for clean temporal dominance + ("RK2/WENO5", ["--order", "5", "--time-stepper", "2"], 2, 0.2, [0.50, 0.25]), + ("RK3/WENO5", ["--order", "5", "--time-stepper", "3"], 3, 0.3, [0.50, 0.25]), +] + +N_SPATIAL = 512 # fixed spatial resolution + + +def read_cons_var(run_dir: str, step: int, var_idx: int, num_ranks: int = 1, expected_size: int = None) -> np.ndarray: + """Read q_cons_vf{var_idx} from all MPI ranks and concatenate into one 1D array.""" + 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] + data = np.frombuffer(f.read(rec_len), dtype=np.float64) + f.read(4) + chunks.append(data.copy()) + combined = np.concatenate(chunks) + if expected_size is not None and combined.size != expected_size: + raise ValueError(f"Expected {expected_size} values across {num_ranks} ranks, got {combined.size}") + return combined + + +# 1D single-fluid Euler (model_eqns=2, num_fluids=1): vf1=ρ, vf2=ρu, vf3=E +CONS_VARS_1D = [("density", 1), ("x-momentum", 2), ("energy", 3)] +CONS_TOL = 1e-10 + + +def conservation_errors(run_dir: str, Nt: int, cell_vol: float, var_list: list, num_ranks: int, expected_size: int = None) -> dict: + """Return relative conservation error |Σq(T) - Σq(0)| / |Σq(0)| for each variable.""" + errs = {} + for name, idx in var_list: + q0 = read_cons_var(run_dir, 0, idx, num_ranks, expected_size) + qT = read_cons_var(run_dir, Nt, idx, num_ranks, expected_size) + s0 = float(np.sum(q0)) * cell_vol + sT = float(np.sum(qT)) * cell_vol + errs[name] = abs(sT - s0) / (abs(s0) + 1e-300) + return errs + + +def l2_error(a: np.ndarray, b: np.ndarray, dx: float) -> float: + return float(np.sqrt(np.sum((a - b) ** 2) * dx)) + + +def run_case(tmpdir: str, cfl: float, extra_args: list, num_ranks: int = 1): + """Run the 1D advection case at fixed N=512 with given CFL. Returns (dt, Nt, run_dir).""" + result = subprocess.run( + [sys.executable, CASE, "--mfc", "{}", "-N", str(N_SPATIAL), "--cfl", str(cfl)] + extra_args, + capture_output=True, + text=True, + check=False, + ) + if result.returncode != 0: + raise RuntimeError(f"case.py failed:\n{result.stderr}") + cfg = json.loads(result.stdout) + Nt = int(cfg["t_step_stop"]) + dt = float(cfg["dt"]) + + cmd = [ + MFC, + "run", + CASE, + "-t", + "pre_process", + "simulation", + "-n", + str(num_ranks), + "--", + "-N", + str(N_SPATIAL), + "--cfl", + str(cfl), + ] + extra_args + result = subprocess.run(cmd, capture_output=True, text=True, cwd=os.getcwd(), check=False) + if result.returncode != 0: + print(result.stdout[-3000:]) + print(result.stderr) + raise RuntimeError(f"./mfc.sh run failed for CFL={cfl}") + + case_dir = os.path.dirname(CASE) + src = os.path.join(case_dir, "p_all") + cfl_tag = f"cfl{cfl:.4f}".replace(".", "p") + dst = os.path.join(tmpdir, cfl_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 dt, Nt, os.path.join(tmpdir, cfl_tag) + + +def test_scheme(label, extra_args, expected_order, tol, cfls, num_ranks=1): + print(f"\n{'=' * 60}") + print(f" {label} N={N_SPATIAL} (need rate >= {expected_order - tol:.1f})") + print(f"{'=' * 60}") + + errors = [] + dts = [] + nts = [] + dx = 1.0 / N_SPATIAL + all_cons_errs = [] + + with tempfile.TemporaryDirectory() as tmpdir: + for cfl in cfls: + dt, Nt, run_dir = run_case(tmpdir, cfl, extra_args, num_ranks) + dts.append(dt) + nts.append(Nt) + vf0 = read_cons_var(run_dir, 0, 1, num_ranks, expected_size=N_SPATIAL) + vfT = read_cons_var(run_dir, Nt, 1, num_ranks, expected_size=N_SPATIAL) + err = l2_error(vfT, vf0, dx) + errors.append(err) + all_cons_errs.append(conservation_errors(run_dir, Nt, dx, CONS_VARS_1D, num_ranks, expected_size=N_SPATIAL)) + + rates = [None] + for i in range(1, len(cfls)): + log_dt0 = math.log(dts[i - 1]) + log_dt1 = math.log(dts[i]) + rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (log_dt1 - log_dt0)) + + print(f" {'CFL':>7} {'dt':>12} {'Nt':>6} {'L2 error':>14} {'rate':>8}") + print(f" {'-' * 7} {'-' * 12} {'-' * 6} {'-' * 14} {'-' * 8}") + for i, cfl in enumerate(cfls): + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {cfl:>7.3f} {dts[i]:>12.6e} {nts[i]:>6} {errors[i]:>14.6e} {r_str}") + + if len(cfls) > 1: + log_dt = np.log(np.array(dts, dtype=float)) + log_err = np.log(np.array(errors, dtype=float)) + overall, _ = np.polyfit(log_dt, log_err, 1) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + rate_passed = overall >= expected_order - tol + else: + print("\n (need >= 2 CFL values to compute rate)") + rate_passed = True + + print(f"\n Conservation (need rel. error < {CONS_TOL:.0e}):") + cons_passed = True + for name, _ in CONS_VARS_1D: + max_err = max(ce[name] for ce in all_cons_errs) + ok = max_err < CONS_TOL + print(f" {name:<14}: max = {max_err:.2e} {'OK' if ok else 'FAIL'}") + if not ok: + cons_passed = False + + passed = rate_passed and cons_passed + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC RK3 temporal order verification") + parser.add_argument( + "--cfls", + type=float, + nargs="+", + default=None, + help="CFL values to test (default: per-scheme values)", + ) + parser.add_argument( + "--schemes", + nargs="+", + default=["RK1/WENO5", "RK2/WENO5", "RK3/WENO5"], + help="Schemes to test (default: all)", + ) + parser.add_argument("--num-ranks", type=int, default=1, help="MPI ranks per simulation (default: 1)") + args = parser.parse_args() + + results = {} + for label, extra_args, expected_order, tol, default_cfls in SCHEMES: + if label not in args.schemes: + continue + cfls = args.cfls if args.cfls is not None else default_cfls + try: + passed = test_scheme(label, extra_args, expected_order, tol, cfls, args.num_ranks) + except Exception as e: + import traceback + + print(f" ERROR: {e}") + traceback.print_exc() + passed = False + results[label] = passed + + print(f"\n{'=' * 60}") + print(" Summary") + print(f"{'=' * 60}") + all_pass = True + for label, passed in results.items(): + print(f" {label:<14} {'PASS' if passed else 'FAIL'}") + if not passed: + all_pass = False + + sys.exit(0 if all_pass else 1) + + +if __name__ == "__main__": + main()