From aad533c475f837f2cb6de47691ec1badba5f4fe9 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Sun, 9 Nov 2025 14:01:27 -0500 Subject: [PATCH 1/4] add 1D convergence test example --- examples/1D_convergence/README.md | 11 ++++ examples/1D_convergence/case.py | 83 +++++++++++++++++++++++++++++++ toolchain/mfc/case.py | 2 +- 3 files changed, 95 insertions(+), 1 deletion(-) create mode 100644 examples/1D_convergence/README.md create mode 100755 examples/1D_convergence/case.py diff --git a/examples/1D_convergence/README.md b/examples/1D_convergence/README.md new file mode 100644 index 0000000000..29ddf4dbab --- /dev/null +++ b/examples/1D_convergence/README.md @@ -0,0 +1,11 @@ +This example case contains an automated convergence test using a 1D, two-component advection case. +The case can be ran by executing the bash script `./submitJobs.sh` in a terminal after enabling +execution permissions with `chmod +x ./submitJobs.sh` and setting the `ROOT_DIR` and `MFC_DIR` +variables. The default settings the script run the case for 6 different grid resolutions with 1st, +3rd, and 5th, order spatial reconstructions. These settings can be modified by editing the variables +at the top of the script. You can also run different model equations by setting the `ME` variable +and different Riemann solvers by setting the `RS` variable. + +Once the simulations have been run, you can generate convergence plots with matplotlib by running +`python3 plot.py` in a terminal. This will generate plots of the L1, L2, and Linf error norms and +save the results to `errors.csv`. diff --git a/examples/1D_convergence/case.py b/examples/1D_convergence/case.py new file mode 100755 index 0000000000..12b8b341d2 --- /dev/null +++ b/examples/1D_convergence/case.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python3 +import math +import json +import argparse + +# Parsing command line arguments +parser = argparse.ArgumentParser( + description="Generate JSON case configuration for two-fluid convergence simulation." +) +parser.add_argument("--mfc", type=json.loads,default="{}",metavar="DICT",help="MFC's toolchain's internal state.") +parser.add_argument("--order", type=int, default=5, help="WENO order (default: 5)") +parser.add_argument("--meqns", type=int, default=2, help="Model equations (default: 2 (five-equation model))") +parser.add_argument("--rs", type=int, default=2, help="Riemann solver (default: 2 (HLLC))") +parser.add_argument("-N", type=int, default=1024, help="Number of grid points (default: 1024)") + +args = parser.parse_args() + +# Numerical setup +Nx = args.N - 1 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend, Nt = 5.0, 200000 +mydt = Tend / (1.0 * Nt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0.0e00, + "x_domain%end": 1.0e00, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(math.ceil(Nt / 10.0)), + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": args.meqns, + "alt_soundspeed": "F", + "num_fluids": 2, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": args.order, + "weno_eps": dx**2, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": args.rs, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Patch 1 L + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%vel(1)": 1.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": f"0.5 - 0.5*sin(2*pi*x)", + "patch_icpp(1)%alpha(1)": f"0.5 - 0.5*sin(2*pi*x)", + "patch_icpp(1)%alpha_rho(2)": f"0.5 + 0.5*sin(2*pi*x)", + "patch_icpp(1)%alpha(2)": f"0.5 + 0.5*sin(2*pi*x)", + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + "fluid_pp(2)%gamma": 1.0e00 / (1.4 - 1.0e00), + "fluid_pp(2)%pi_inf": 0.0, + } + ) +) diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index fb30a21f61..269c8438bf 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -142,7 +142,7 @@ def rhs_replace(match): 'r': f'patch_icpp({pid})%radius', 'eps': f'patch_icpp({pid})%epsilon', 'beta': f'patch_icpp({pid})%beta', 'tau_e': f'patch_icpp({pid})%tau_e', 'radii': f'patch_icpp({pid})%radii', - 'e' : f'{math.e}', 'pi': f'{math.pi}', + 'e' : f'{math.e}', }.get(match.group(), match.group()) lines = [] From 87e0dcb6b5f7a50d26a3de62129e404fe5c73aff Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Sun, 9 Nov 2025 14:02:07 -0500 Subject: [PATCH 2/4] format --- examples/1D_convergence/case.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/1D_convergence/case.py b/examples/1D_convergence/case.py index 12b8b341d2..fe1ed3c976 100755 --- a/examples/1D_convergence/case.py +++ b/examples/1D_convergence/case.py @@ -4,10 +4,8 @@ import argparse # Parsing command line arguments -parser = argparse.ArgumentParser( - description="Generate JSON case configuration for two-fluid convergence simulation." -) -parser.add_argument("--mfc", type=json.loads,default="{}",metavar="DICT",help="MFC's toolchain's internal state.") +parser = argparse.ArgumentParser(description="Generate JSON case configuration for two-fluid convergence simulation.") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC's toolchain's internal state.") parser.add_argument("--order", type=int, default=5, help="WENO order (default: 5)") parser.add_argument("--meqns", type=int, default=2, help="Model equations (default: 2 (five-equation model))") parser.add_argument("--rs", type=int, default=2, help="Riemann solver (default: 2 (HLLC))") From f1c26f9b38098b0af8de559cb4796d56570c2e6b Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Sun, 9 Nov 2025 14:04:15 -0500 Subject: [PATCH 3/4] add missing files --- examples/1D_convergence/plot.py | 70 +++++++++++++++++++++++++++ examples/1D_convergence/submitJobs.sh | 27 +++++++++++ 2 files changed, 97 insertions(+) create mode 100644 examples/1D_convergence/plot.py create mode 100755 examples/1D_convergence/submitJobs.sh diff --git a/examples/1D_convergence/plot.py b/examples/1D_convergence/plot.py new file mode 100644 index 0000000000..0f239b5fae --- /dev/null +++ b/examples/1D_convergence/plot.py @@ -0,0 +1,70 @@ +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +N = [32, 64, 128, 256, 512, 1024] +Ord = [1, 3, 5] + +errors = np.nan * np.zeros((len(N), len(Ord), 3)) + +TEND = 200000 + +for i in range(len(N)): + for j in range(len(Ord)): + + sim_a1 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.5.00.{TEND}.dat", sep=r"\s+", header=None, names=["x", "y"]) + sim_a2 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.6.00.{TEND}.dat", sep=r"\s+", header=None, names=["x", "y"]) + + exact_a1 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.5.00.000000.dat", sep=r"\s+", header=None, names=["x", "y"]) + exact_a2 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.6.00.000000.dat", sep=r"\s+", header=None, names=["x", "y"]) + + ## 2 norm + errors[i, j, 0] = 1 / N[i] * np.sum(np.sqrt((sim_a1.y - exact_a1.y) ** 2)) + errors[i, j, 0] += 1 / N[i] * np.sum(np.sqrt((sim_a2.y - exact_a2.y) ** 2)) + + ## 1 norm + errors[i, j, 1] = 1 / N[i] * np.sum(np.abs(sim_a1.y - exact_a1.y)) + errors[i, j, 1] += 1 / N[i] * np.sum(np.abs(sim_a2.y - exact_a2.y)) + + ## Inf norm + errors[i, j, 2] = np.nanmax(np.abs(sim_a1.y - exact_a1.y)) + errors[i, j, 2] += np.nanmax(np.abs(sim_a2.y - exact_a2.y)) + +fig, ax = plt.subplots(1, 3, figsize=(12, 8), sharex=True) + +colors = ["blue", "green", "red", "purple"] + +ref = np.nan * np.zeros((len(N), len(Ord))) + +for i in range(3): + ax[i].plot(N, 30 / np.array(N) ** 1, label="Slope = -1", color=colors[0]) + ax[i].plot(N, 3000 / np.array(N) ** 3, label="Slope = -3", color=colors[1]) + ax[i].plot(N, 5000 / np.array(N) ** 5, label="Slope = -5", color=colors[2]) + +for j in range(len(Ord)): + ax[0].plot(N, errors[:, j, 0], "o-", color=colors[j]) + ax[0].set_xscale("log", base=2) + ax[0].set_yscale("log") + ax[0].set_title("||error||_2") + ax[0].legend() + + ax[1].plot(N, errors[:, j, 1], "o-", color=colors[j]) + ax[1].set_xscale("log", base=2) + ax[1].set_yscale("log") + ax[1].set_title("||error||_1") + + ax[2].plot(N, errors[:, j, 2], "o-", color=colors[j]) + ax[2].set_xscale("log", base=2) + ax[2].set_yscale("log") + ax[2].set_title("||error||_inf") + +plt.tight_layout() +plt.show() + +errors = np.column_stack((N, errors[:, :, 0])) +errors = np.column_stack((errors, 7 / np.array(N) ** 1)) +errors = np.column_stack((errors, 700 / np.array(N) ** 3)) +errors = np.column_stack((errors, 2000 / np.array(N) ** 5)) + +df = pd.DataFrame(errors, columns=["N", "1", "3", "5", "R1", "R3", "R5"], index=N) +df.to_csv("errors.csv", index=False) diff --git a/examples/1D_convergence/submitJobs.sh b/examples/1D_convergence/submitJobs.sh new file mode 100755 index 0000000000..174c7b8cc2 --- /dev/null +++ b/examples/1D_convergence/submitJobs.sh @@ -0,0 +1,27 @@ +#!/usr/bin/env sh +Nx=(32 64 128 256 512 1024) +Order=(1 3 5) +ME=2 # Model equations = 2 for five-equation model +RS=2 # Riemann solver = 2 for HLLC + +#ROOT_DIR= +#MFC_DIR= +ROOT_DIR="/Users/benwilfong/Documents/software/MFC-Wilfong/examples/1D_convergence" +MFC_DIR="/Users/benwilfong/Documents/software/MFC-Wilfong" + +for i in "${Nx[@]}"; do + for j in "${Order[@]}"; do + rm -rf N${i}_O${j} + mkdir N${i}_O${j} + cp case.py N${i}_O${j}/ + done +done + +cd $MFC_DIR + +for i in "${Nx[@]}"; do + for j in "${Order[@]}"; do + ./mfc.sh run $ROOT_DIR/N${i}_O${j}/case.py --case-optimization --no-debug -- --order $j -N $i --meqns $ME --rs $RS + done +done + From 3dd46e2b404a523b1111dbcece8627bcaff154b2 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Sun, 9 Nov 2025 17:35:36 -0500 Subject: [PATCH 4/4] Github suggestions and bug fix --- examples/1D_convergence/README.md | 4 ++-- examples/1D_convergence/plot.py | 7 +++---- examples/1D_convergence/submitJobs.sh | 9 +++++---- toolchain/mfc/test/cases.py | 3 ++- 4 files changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/1D_convergence/README.md b/examples/1D_convergence/README.md index 29ddf4dbab..3b9963e048 100644 --- a/examples/1D_convergence/README.md +++ b/examples/1D_convergence/README.md @@ -1,7 +1,7 @@ This example case contains an automated convergence test using a 1D, two-component advection case. -The case can be ran by executing the bash script `./submitJobs.sh` in a terminal after enabling +The case can be run by executing the bash script `./submitJobs.sh` in a terminal after enabling execution permissions with `chmod +x ./submitJobs.sh` and setting the `ROOT_DIR` and `MFC_DIR` -variables. The default settings the script run the case for 6 different grid resolutions with 1st, +variables. By default the script runs the case for 6 different grid resolutions with 1st, 3rd, and 5th, order spatial reconstructions. These settings can be modified by editing the variables at the top of the script. You can also run different model equations by setting the `ME` variable and different Riemann solvers by setting the `RS` variable. diff --git a/examples/1D_convergence/plot.py b/examples/1D_convergence/plot.py index 0f239b5fae..60747b6965 100644 --- a/examples/1D_convergence/plot.py +++ b/examples/1D_convergence/plot.py @@ -19,16 +19,15 @@ exact_a2 = pd.read_csv(f"N{N[i]}_O{Ord[j]}/D/cons.6.00.000000.dat", sep=r"\s+", header=None, names=["x", "y"]) ## 2 norm - errors[i, j, 0] = 1 / N[i] * np.sum(np.sqrt((sim_a1.y - exact_a1.y) ** 2)) - errors[i, j, 0] += 1 / N[i] * np.sum(np.sqrt((sim_a2.y - exact_a2.y) ** 2)) + errors[i, j, 0] = np.linalg.norm(sim_a1.y - exact_a1.y) / np.sqrt(N[i]) + errors[i, j, 0] += np.linalg.norm(sim_a2.y - exact_a2.y) / np.sqrt(N[i]) ## 1 norm errors[i, j, 1] = 1 / N[i] * np.sum(np.abs(sim_a1.y - exact_a1.y)) errors[i, j, 1] += 1 / N[i] * np.sum(np.abs(sim_a2.y - exact_a2.y)) ## Inf norm - errors[i, j, 2] = np.nanmax(np.abs(sim_a1.y - exact_a1.y)) - errors[i, j, 2] += np.nanmax(np.abs(sim_a2.y - exact_a2.y)) + errors[i, j, 2] = np.max([np.nanmax(np.abs(sim_a1.y - exact_a1.y)), np.nanmax(np.abs(sim_a2.y - exact_a2.y))]) fig, ax = plt.subplots(1, 3, figsize=(12, 8), sharex=True) diff --git a/examples/1D_convergence/submitJobs.sh b/examples/1D_convergence/submitJobs.sh index 174c7b8cc2..0f5315fa9e 100755 --- a/examples/1D_convergence/submitJobs.sh +++ b/examples/1D_convergence/submitJobs.sh @@ -1,13 +1,14 @@ #!/usr/bin/env sh Nx=(32 64 128 256 512 1024) Order=(1 3 5) + ME=2 # Model equations = 2 for five-equation model RS=2 # Riemann solver = 2 for HLLC -#ROOT_DIR= -#MFC_DIR= -ROOT_DIR="/Users/benwilfong/Documents/software/MFC-Wilfong/examples/1D_convergence" -MFC_DIR="/Users/benwilfong/Documents/software/MFC-Wilfong" +# Get the directory of the script itself +ROOT_DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )" +# Assume the script is in examples/1D_convergence, so MFC_DIR is two levels up +MFC_DIR="$(dirname "$(dirname "$ROOT_DIR")")" for i in "${Nx[@]}"; do for j in "${Order[@]}"; do diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index 6fb00781be..fa23933b69 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -1021,7 +1021,8 @@ def foreach_example(): "3D_TaylorGreenVortex_analytical", "3D_IGR_TaylorGreenVortex_nvidia", "2D_backward_facing_step", - "2D_forward_facing_step"] + "2D_forward_facing_step", + "1D_convergence"] if path in casesToSkip: continue name = f"{path.split('_')[0]} -> Example -> {'_'.join(path.split('_')[1:])}"