diff --git a/.github/workflows/run-dumux-benchmark.yml b/.github/workflows/run-dumux-benchmark.yml index 04dddf5..d7eacaa 100644 --- a/.github/workflows/run-dumux-benchmark.yml +++ b/.github/workflows/run-dumux-benchmark.yml @@ -11,7 +11,7 @@ env: CACHE_NUMBER: 1 # increase to reset cache manually jobs: - tests: + dumux-tests: runs-on: ubuntu-latest steps: @@ -41,11 +41,10 @@ jobs: run: | cd $GITHUB_WORKSPACE/benchmarks/rotating-cylinders/dumux python3 run_benchmark.py - python3 convergence_test.py - name: Upload results uses: actions/upload-artifact@v4 with: - name: convergence-results - path: benchmarks/rotating-cylinders/dumux/rotating-cylinders/results - retention-days: 5 \ No newline at end of file + name: results + path: benchmarks/rotating-cylinders/dumux/results + retention-days: 5 diff --git a/.github/workflows/run-openfoam-benchmark.yml b/.github/workflows/run-openfoam-benchmark.yml index ee7dac3..6f0db02 100644 --- a/.github/workflows/run-openfoam-benchmark.yml +++ b/.github/workflows/run-openfoam-benchmark.yml @@ -11,7 +11,7 @@ env: CACHE_NUMBER: 1 # increase to reset cache manually jobs: - openfoam: + openfoam-tests: runs-on: ubuntu-latest steps: @@ -45,6 +45,6 @@ jobs: - name: Upload OpenFOAM results uses: actions/upload-artifact@v4 with: - name: openfoam-results + name: results path: benchmarks/rotating-cylinders/openfoam/results retention-days: 5 diff --git a/benchmarks/rotating-cylinders/dumux/Snakefile b/benchmarks/rotating-cylinders/dumux/Snakefile index a0a8fe9..cc6cdf9 100644 --- a/benchmarks/rotating-cylinders/dumux/Snakefile +++ b/benchmarks/rotating-cylinders/dumux/Snakefile @@ -11,15 +11,13 @@ rule all: input: "solution_metrics.json", "test_rotatingcylinders.zip" - # "test_rotatingcylinders-00000.vtu" rule run_simulation: input: rc_parameters_file = "parameters.json" output: zip = "test_rotatingcylinders.zip", - metrics = "solution_metrics.json", - # "test_rotatingcylinders-00000.vtu" + metrics = "solution_metrics.json" resources: serial_run=1 singularity: diff --git a/benchmarks/rotating-cylinders/dumux/convergence_test.py b/benchmarks/rotating-cylinders/dumux/convergence_test.py deleted file mode 100644 index b181cf6..0000000 --- a/benchmarks/rotating-cylinders/dumux/convergence_test.py +++ /dev/null @@ -1,77 +0,0 @@ -import json -import zipfile -import pandas as pd -import matplotlib.pyplot as plt -from pathlib import Path - -# Setup paths -root_dir = Path(__file__).resolve().parent -results_base_dir = root_dir / "rotating-cylinders" / "results" - -data = [] - -if not results_base_dir.exists(): - print(f"Error: Base directory {results_base_dir} does not exist.") - exit(1) - -# Iterate through every configuration folder -for conf_dir in results_base_dir.iterdir(): - if conf_dir.is_dir(): - zip_path = conf_dir / "test_rotatingcylinders.zip" - - # 1. Extraction: Unzip everything in the folder - if zip_path.exists(): - try: - with zipfile.ZipFile(zip_path, 'r') as zip_ref: - zip_ref.extractall(conf_dir) - print(f"Extracted zip for: {conf_dir.name}") - except zipfile.BadZipFile: - print(f"Error: {zip_path.name} is corrupt.") - - # 2. Data Collection: Look for the newly extracted JSON - metrics_path = conf_dir / "solution_metrics.json" - if metrics_path.exists(): - try: - with open(metrics_path, 'r') as f: - metrics = json.load(f) - - data.append({ - "conf": conf_dir.name, - # "pressure_error": metrics.get("l2_error_pressure_rel"), - "velocity_error": metrics.get("l2_error_velocity_rel") - }) - except (json.JSONDecodeError, IOError) as e: - print(f"Error reading metrics in {conf_dir.name}: {e}") - -# 3. Processing and Plotting -if not data: - print("No metrics found. Ensure the zip files contain solution_metrics.json.") - exit(1) - -df = pd.DataFrame(data) - -# Sort numerically if folder names are numbers, otherwise alphabetical -try: - df['conf_sort'] = pd.to_numeric(df['conf']) - df = df.sort_values('conf_sort') -except ValueError: - df = df.sort_values('conf') - -plt.figure(figsize=(10, 6)) -# plt.plot(df['conf'], df['pressure_error'], marker='o', label='$L^2$ Rel. Pressure Error') -plt.plot(df['conf'], df['velocity_error'], marker='s', label='Relative $L^2$ Error for Velocity') - -plt.yscale('log') -plt.xlabel('Configuration') -plt.ylabel('Relative Error') -plt.title('Convergence Summary') -plt.legend() -plt.grid(True, which="both", ls="-", alpha=0.3) -plt.xticks(rotation=45) -plt.tight_layout() - -# 4. Save the plot specifically into the results directory -plot_output = results_base_dir / "solution_metrics_plot.png" -plt.savefig(plot_output) - -print(f"\nSummary plot saved to: {plot_output}") \ No newline at end of file diff --git a/benchmarks/rotating-cylinders/dumux/plot_results.py b/benchmarks/rotating-cylinders/dumux/plot_results.py new file mode 100644 index 0000000..83a27af --- /dev/null +++ b/benchmarks/rotating-cylinders/dumux/plot_results.py @@ -0,0 +1,653 @@ +#!/usr/bin/env pvpython +""" +plot_results.py +================ + +Unified post-processing suite for rotating-cylinder benchmark cases. + +Features +-------- +1. Global convergence plot +2. Optional ParaView 2D field rendering +3. Optional analytical comparison plots +4. Single-configuration filtering +5. Custom results directory support + +Examples +-------- +# Process ALL cases (convergence plot only) +pvpython3 plot_results.py + +# Process a SPECIFIC configuration +pvpython3 plot_results.py --conf 100_0_80_320 + +# Also generate 2D field maps +pvpython3 plot_results.py --conf 100_0_80_320 --fields + +# Also generate analytical comparison plot +pvpython3 plot_results.py --conf 100_0_80_320 --analytical + +# Generate everything +pvpython3 plot_results.py --conf 100_0_80_320 --fields --analytical + +# Explicit results directory +pvpython3 plot_results.py --conf 100_0_80_320 --fields --analytical /path/to/results +""" + +from __future__ import annotations + +import argparse +import re +import tempfile +import zipfile +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +# ParaView +from paraview import servermanager +from paraview.simple import OpenDataFile, Delete + +# ============================================================================= +# Configuration +# ============================================================================= + +ROOT_DIR = Path(__file__).resolve().parent +DEFAULT_RESULTS_DIR = ROOT_DIR / "results" + +DPI = 150 +VIEW_SIZE = [2000, 1600] + +# ============================================================================= +# Case Discovery +# ============================================================================= + + +def discover_cases(results_dir: Path, conf: str | None = None): + """ + Discover benchmark case directories. + """ + + cases = [] + + for d in sorted(results_dir.iterdir()): + + if not d.is_dir(): + continue + + if conf and d.name != conf: + continue + + zip_file = d / "test_rotatingcylinders.zip" + + if not zip_file.exists(): + continue + + tokens = d.name.split("_") + + if len(tokens) != 4: + continue + + omega1, omega2, nr, ntheta = map(float, tokens) + + cases.append({ + "path": d, + "zip": zip_file, + "omega1": omega1, + "omega2": omega2, + "nr": int(nr), + "ntheta": int(ntheta), + }) + + return cases + + +# ============================================================================= +# Extraction Utilities +# ============================================================================= + + +def extract_pvd_and_vtus(zip_path: Path, tmpdir: Path): + """ + Extract PVD and VTU files from archive and clean paths inside PVD. + """ + + with zipfile.ZipFile(zip_path, "r") as zf: + + for item in zf.namelist(): + + filename = Path(item).name + + if item.endswith(".pvd") or item.endswith(".vtu"): + + with zf.open(item) as source, open(tmpdir / filename, "wb") as target: + target.write(source.read()) + + pvd_files = list(tmpdir.glob("*.pvd")) + + if not pvd_files: + raise RuntimeError(f"No PVD file found in: {zip_path}") + + pvd_path = pvd_files[0] + + pvd_content = pvd_path.read_text() + + cleaned_content = re.sub( + r'file="[^"]*/([^"]+\.vtu)"', + r'file="\1"', + pvd_content, + ) + + pvd_path.write_text(cleaned_content) + + return pvd_path + + +# ============================================================================= +# Data Extraction +# ============================================================================= + + +def parse_data_via_vtk(pvd_filepath: Path): + """ + Extract pressure, velocity, and cell centers using ParaView. + """ + + data_reader = OpenDataFile(str(pvd_filepath)) + data_reader.UpdatePipeline(1.0) + + client_data = servermanager.Fetch(data_reader) + + vtk_p_array = client_data.GetCellData().GetArray("p") + vtk_v_array = client_data.GetCellData().GetArray("velocity_liq (m/s)") + + num_cells = client_data.GetNumberOfCells() + + pressure = np.array([ + vtk_p_array.GetValue(i) + for i in range(num_cells) + ]) + + velocity = np.array([ + vtk_v_array.GetTuple(i) + for i in range(num_cells) + ]) + + cell_centers = [] + + for i in range(num_cells): + + cell = client_data.GetCell(i) + pts = cell.GetPoints() + + num_pts = pts.GetNumberOfPoints() + + coords = np.array([ + pts.GetPoint(j) + for j in range(num_pts) + ]) + + cell_centers.append(coords.mean(axis=0)) + + cell_centers = np.array(cell_centers) + + x = cell_centers[:, 0] + y = cell_centers[:, 1] + + Delete(data_reader) + + return x, y, pressure, velocity + + +# ============================================================================= +# Analytical Solution +# ============================================================================= + + +def analytical_solution(r, omega1, omega2, r1=1.0, r2=2.0): + """ + Taylor-Couette analytical solution. + """ + + A = ( + (omega2 * r2**2 - omega1 * r1**2) + / (r2**2 - r1**2) + ) + + B = ( + ((omega1 - omega2) * r1**2 * r2**2) + / (r2**2 - r1**2) + ) + + u = A * r + B / r + + p = ( + A**2 * r**2 / 2 + + 2 * A * B * np.log(r) + - B**2 / (2 * r**2) + ) + + p -= p[0] + + return u, p + + +# ============================================================================= +# Radial Profiles +# ============================================================================= + + +def compute_radial_profiles( + cx, + cy, + pressure, + velocity, + r1=1.0, + r2=2.0, + n_bins=80, +): + """ + Compute radial averages. + """ + + vx = velocity[:, 0] + vy = velocity[:, 1] + + speed = np.sqrt(vx**2 + vy**2) + + r_cells = np.sqrt(cx**2 + cy**2) + + mask = (r_cells >= r1) & (r_cells <= r2) + + r_filtered = r_cells[mask] + p_filtered = pressure[mask] + s_filtered = speed[mask] + + bins = np.linspace(r1, r2, n_bins + 1) + + bc = 0.5 * (bins[:-1] + bins[1:]) + + idx = np.clip( + np.digitize(r_filtered, bins) - 1, + 0, + n_bins - 1, + ) + + p_bin = np.array([ + p_filtered[idx == i].mean() + if (idx == i).any() + else np.nan + for i in range(n_bins) + ]) + + s_bin = np.array([ + s_filtered[idx == i].mean() + if (idx == i).any() + else np.nan + for i in range(n_bins) + ]) + + if not np.isnan(p_bin).all(): + p_bin -= p_bin[~np.isnan(p_bin)][0] + + return bc, s_bin, p_bin + + +# ============================================================================= +# Analytical Comparison Plot +# ============================================================================= + + +def generate_analytical_plot(case, r_num, u_num, p_num): + """ + Generate analytical comparison plot. + """ + + u_ref, p_ref = analytical_solution( + r_num, + case["omega1"], + case["omega2"], + ) + + fig, (ax1, ax2) = plt.subplots( + 1, + 2, + figsize=(12, 5), + dpi=DPI, + ) + + # Velocity + ax1.plot( + r_num, + u_ref, + lw=2, + color="black", + label="Analytical", + ) + + ax1.plot( + r_num, + u_num, + "o", + ms=4, + mfc="none", + color="red", + label="DuMux", + ) + + ax1.set_title("Velocity Profile") + ax1.set_xlabel("r") + ax1.set_ylabel("u_theta") + ax1.legend() + ax1.grid(True, alpha=0.3) + + # Pressure + ax2.plot( + r_num, + p_ref, + lw=2, + color="black", + label="Analytical", + ) + + ax2.plot( + r_num, + p_num, + "s", + ms=4, + mfc="none", + color="tab:blue", + label="DuMux", + ) + + ax2.set_title("Pressure Profile") + ax2.set_xlabel("r") + ax2.set_ylabel("p") + ax2.legend() + ax2.grid(True, alpha=0.3) + + plt.tight_layout() + + outpath = case["path"] / "analytical_comparison.png" + + fig.savefig(outpath, bbox_inches="tight") + + plt.close(fig) + + print(f" Saved analytical comparison: {outpath}") + + +# ============================================================================= +# ParaView Rendering +# ============================================================================= + + +import matplotlib.tri as mtri + + +def render_matplotlib_fields( + cx, + cy, + pressure, + velocity, + outdir: Path, + r1=1.0, + r2=2.0, +): + """ + Generate pressure_field.png and velocity_field.png + using pure matplotlib triangulation. + """ + + vx = velocity[:, 0] + vy = velocity[:, 1] + + speed = np.sqrt(vx**2 + vy**2) + + triang = mtri.Triangulation(cx, cy) + + # Mask triangles inside inner cylinder + tri_cx = cx[triang.triangles].mean(axis=1) + tri_cy = cy[triang.triangles].mean(axis=1) + + tri_r = np.sqrt(tri_cx**2 + tri_cy**2) + + triang.set_mask(tri_r < r1) + + theta = np.linspace(0, 2*np.pi, 400) + + fields = [ + ( + pressure, + "p", + "coolwarm", + "pressure_field.png", + ), + ( + speed, + "|u|", + "viridis", + "velocity_field.png", + ), + ] + + for field_vals, label, cmap, filename in fields: + + fig, ax = plt.subplots( + figsize=(7, 6), + dpi=DPI, + ) + + tc = ax.tricontourf( + triang, + field_vals, + levels=64, + cmap=cmap, + ) + + fig.colorbar( + tc, + ax=ax, + label=label, + fraction=0.046, + pad=0.04, + ) + + # Inner cylinder + ax.plot( + r1*np.cos(theta), + r1*np.sin(theta), + "k-", + lw=1.2, + ) + + # Outer cylinder + ax.plot( + r2*np.cos(theta), + r2*np.sin(theta), + "k-", + lw=1.2, + ) + + # White fill inside hole + from matplotlib.patches import Circle + + ax.add_patch( + Circle( + (0, 0), + r1, + color="white", + zorder=3, + ) + ) + + ax.plot( + r1*np.cos(theta), + r1*np.sin(theta), + "k-", + lw=1.2, + zorder=4, + ) + + ax.set_aspect("equal") + + ax.set_xlabel("x") + ax.set_ylabel("y") + + ax.set_title(label) + + ax.tick_params( + which="both", + direction="in", + ) + + plt.tight_layout() + + outpath = outdir / filename + + fig.savefig( + outpath, + bbox_inches="tight", + ) + + plt.close(fig) + + print(f" Saved: {outpath}") + +# ============================================================================= +# Main +# ============================================================================= + + +def main(): + + parser = argparse.ArgumentParser( + description="Unified rotating-cylinder post-processing suite" + ) + + parser.add_argument( + "results_dir", + nargs="?", + default=DEFAULT_RESULTS_DIR, + type=Path, + help="Path to results directory", + ) + + parser.add_argument( + "--conf", + type=str, + help="Process only a specific configuration", + ) + + parser.add_argument( + "--fields", + action="store_true", + help="Generate pressure and velocity field maps", + ) + + parser.add_argument( + "--analytical", + action="store_true", + help="Generate analytical comparison plots", + ) + + args = parser.parse_args() + + results_dir = args.results_dir.resolve() + + if not results_dir.exists(): + raise RuntimeError( + f"Results directory not found: {results_dir}" + ) + + cases = discover_cases(results_dir, args.conf) + + if not cases: + raise RuntimeError( + f"No matching cases found inside: {results_dir}" + ) + + metrics_data = [] + + print("=" * 80) + print(f"Results directory: {results_dir}") + + if args.conf: + print(f"Selected configuration: {args.conf}") + + print("=" * 80) + + # ========================================================================= + # Process Cases + # ========================================================================= + + for case in cases: + + print(f"\nProcessing Case: {case['path'].name}") + + with tempfile.TemporaryDirectory() as tmp: + + tmpdir = Path(tmp) + + pvd_path = extract_pvd_and_vtus( + case["zip"], + tmpdir, + ) + + # ----------------------------------------------------------------- + # Extract data + # ----------------------------------------------------------------- + + cx, cy, pressure, velocity = parse_data_via_vtk( + pvd_path + ) + + # ----------------------------------------------------------------- + # Radial profiles + # ----------------------------------------------------------------- + + r_num, u_num, p_num = compute_radial_profiles( + cx, + cy, + pressure, + velocity, + ) + + # ----------------------------------------------------------------- + # Analytical comparison + # ----------------------------------------------------------------- + + if args.analytical: + + generate_analytical_plot( + case, + r_num, + u_num, + p_num, + ) + + # ----------------------------------------------------------------- + # ParaView field rendering + # ----------------------------------------------------------------- + + if args.fields: + + print(" Rendering pressure field...") + + render_matplotlib_fields( + cx, + cy, + pressure, + velocity, + case["path"], + ) + + +# ============================================================================= +# Entry Point +# ============================================================================= + +if __name__ == "__main__": + main() + diff --git a/benchmarks/rotating-cylinders/dumux/run_benchmark.py b/benchmarks/rotating-cylinders/dumux/run_benchmark.py index d4f37dc..8531710 100644 --- a/benchmarks/rotating-cylinders/dumux/run_benchmark.py +++ b/benchmarks/rotating-cylinders/dumux/run_benchmark.py @@ -1,5 +1,4 @@ import json -import shutil import subprocess import zipfile import os @@ -8,20 +7,19 @@ # We look one level up for the zip and stay in the current folder for extraction. root_dir = Path(__file__).resolve().parent zip_path = root_dir.parent / "rotating-cylinders.zip" -benchmark_dir = root_dir / "rotating-cylinders" snakefile_path = root_dir / "Snakefile" # Extraction if zip_path.exists(): with zipfile.ZipFile(zip_path, 'r') as zip_ref: - zip_ref.extractall(benchmark_dir) - print(f"Successfully extracted benchmark to: {benchmark_dir}") + zip_ref.extractall(root_dir) + print(f"Successfully extracted benchmark to: {root_dir}") else: print(f"Error: Could not find {zip_path} at {root_dir.parent}") exit(1) # Iterate through all parameter files -for param_file in benchmark_dir.glob("parameters_*.json"): +for param_file in root_dir.glob("parameters_*.json"): with open(param_file, "r") as f: data = json.load(f) config_name = data.get("configuration") @@ -31,21 +29,13 @@ continue # Create output directory for the configuration - output_dir = benchmark_dir / "results" / data.get("configuration") + output_dir = root_dir / "results" / data.get("configuration") output_dir.mkdir(parents=True, exist_ok=True) # Copy the selected parameter file to the output directory with a standardised name with open(output_dir / "parameters.json", "w") as outfile: json.dump(data, outfile, indent=2) - - # Copy files from benchmark_dir to output_dir, excluding non-matching parameter files. - for item in benchmark_dir.iterdir(): - if item.is_file(): - if item.name.startswith("parameters_") and item.suffix == ".json": - continue - else: - shutil.copy(item, output_dir / item.name) - + # Run the Snakemake workflow for the configuration subprocess.run([ "snakemake", @@ -53,19 +43,19 @@ "--use-singularity", "--cores", "all", "--resources", "serial_run=1", - "--singularity-args", f"--bind {benchmark_dir}:/dumux/shared", + "--singularity-args", f"--bind {root_dir}:/dumux/shared", "--config", f'conf_name="{config_name}"', "--force" ], check=True, cwd=output_dir) print(f"Workflow executed successfully for {config_name}.") -print("\nAll configurations processed.") +print("\nAll configurations processed.") # --- CLEANUP SECTION --- print("\nStarting cleanup...") # Find all JSON files starting with 'parameters_' in the benchmark directory -for param_file in benchmark_dir.glob("parameters_*.json"): +for param_file in root_dir.glob("parameters_*.json"): try: os.remove(param_file) print(f"Deleted: {param_file.name}") diff --git a/benchmarks/rotating-cylinders/openfoam/0/U b/benchmarks/rotating-cylinders/openfoam/0/U deleted file mode 100644 index eb38d67..0000000 --- a/benchmarks/rotating-cylinders/openfoam/0/U +++ /dev/null @@ -1,39 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class volVectorField; - object U; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -dimensions [0 1 -1 0 0 0 0]; - -internalField uniform (0 0 0); - -boundaryField -{ - innerWall - { - type noSlip; - } - - outerWall - { - type noSlip; - } - - frontAndBack - { - type empty; - } -} - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/0/p b/benchmarks/rotating-cylinders/openfoam/0/p deleted file mode 100644 index 8799498..0000000 --- a/benchmarks/rotating-cylinders/openfoam/0/p +++ /dev/null @@ -1,34 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class volScalarField; - object p; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -dimensions [0 2 -2 0 0 0 0]; - -internalField uniform 0; - -boundaryField -{ - "innerWall|outerWall" - { - type zeroGradient; - } - - frontAndBack - { - type empty; - } -} - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/Allclean b/benchmarks/rotating-cylinders/openfoam/Allclean deleted file mode 100755 index ad9bdea..0000000 --- a/benchmarks/rotating-cylinders/openfoam/Allclean +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/sh -cd "${0%/*}" || exit # Run from this directory -. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions # Tutorial clean functions -#------------------------------------------------------------------------------ - -cleanCase - -rm -f *.png - -#------------------------------------------------------------------------------ diff --git a/benchmarks/rotating-cylinders/openfoam/Allrun b/benchmarks/rotating-cylinders/openfoam/Allrun deleted file mode 100755 index 8cc4e5f..0000000 --- a/benchmarks/rotating-cylinders/openfoam/Allrun +++ /dev/null @@ -1,12 +0,0 @@ -#!/bin/sh -cd "${0%/*}" || exit # Run from this directory -. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions -#------------------------------------------------------------------------------ - -runApplication blockMesh || exit 1 -runApplication $(getApplication) || exit 1 - -echo "Simulation completed. You could now post-process (e.g. visualize) the results." -# ./plot # disabled for benchmark/snakemake-setup - -#------------------------------------------------------------------------------ diff --git a/benchmarks/rotating-cylinders/openfoam/Snakefile b/benchmarks/rotating-cylinders/openfoam/Snakefile index 27ab7dd..768e2b1 100644 --- a/benchmarks/rotating-cylinders/openfoam/Snakefile +++ b/benchmarks/rotating-cylinders/openfoam/Snakefile @@ -4,16 +4,17 @@ from pathlib import Path with open("parameters.json") as f: config = json.load(f) -rule all: - input: "postProcessing/done" +_SNAKEFILE_DIR = Path(workflow.snakefile).resolve().parent +rule all: + input: "solution_metrics.json" rule run_simulation_container: input: template="system/blockMeshDict.template" output: directory("10000"), - touch("postProcessing/done") + "log.simpleFoam" params: r0=config["domain"]["r1"], r1=config["domain"]["r2"], @@ -42,6 +43,13 @@ rule run_simulation_container: ls -la tail log.simpleFoam - - touch postProcessing/done """ + +rule compute_l2_error: + input: + sim_done="log.simpleFoam", + script=ancient(str(_SNAKEFILE_DIR / "compute_l2_error.py")) + output: + "solution_metrics.json" + shell: + "python3 {input.script} {_SNAKEFILE_DIR}/results" diff --git a/benchmarks/rotating-cylinders/openfoam/compute_l2_error.py b/benchmarks/rotating-cylinders/openfoam/compute_l2_error.py new file mode 100644 index 0000000..864476e --- /dev/null +++ b/benchmarks/rotating-cylinders/openfoam/compute_l2_error.py @@ -0,0 +1,325 @@ +#!/usr/bin/env python3 +""" +Compute L2 errors for all rotating-cylinders OpenFOAM configurations. + +Matches the DuMux main.cc approach: + - Full 2D domain integral (all cells, not just a sample line) + - Cell-centre quadrature (FVM cell-average = 1-point quadrature, weighted by cell area) + - Pressure offset: subtract (sim - exact) at the first cell, then compute error + - Velocity: both vx and vy components over all cells + +Cell centres are computed from polyMesh points+faces+owner — no writeCellCentres needed. +Fields are parsed directly from the ASCII OpenFOAM field files. + +Usage: + python3 compute_l2_error.py # auto-finds ./results/ + python3 compute_l2_error.py # explicit path +""" + +import json +import re +import sys +import numpy as np +from pathlib import Path + + +# --------------------------------------------------------------------------- +# OpenFOAM ASCII field reader +# --------------------------------------------------------------------------- + +def _parse_scalar_list(text: str) -> np.ndarray: + """Extract the internalField nonuniform scalar list.""" + m = re.search(r'internalField\s+nonuniform\s+List\s+\d+\s*\(([^)]+)\)', text, re.S) + if not m: + # uniform scalar + m = re.search(r'internalField\s+uniform\s+([\d.eE+\-]+)', text) + if m: + return np.array([float(m.group(1))]) + raise ValueError("Could not parse scalar internalField") + return np.fromstring(m.group(1), sep='\n') + + +def _parse_vector_list(text: str) -> np.ndarray: + """Extract the internalField nonuniform vector list → shape (N, 3).""" + m = re.search(r'internalField\s+nonuniform\s+List\s+\d+\s*\((.+?)\)\s*;', text, re.S) + if not m: + m2 = re.search(r'internalField\s+uniform\s+\(([\d.eE+\-\s]+)\)', text) + if m2: + vals = list(map(float, m2.group(1).split())) + return np.array([vals]) + raise ValueError("Could not parse vector internalField") + rows = re.findall(r'\(([\d.eE+\-\s]+)\)', m.group(1)) + return np.array([list(map(float, r.split())) for r in rows]) + + +def read_scalar_field(path: Path) -> np.ndarray: + return _parse_scalar_list(path.read_text()) + + +def read_vector_field(path: Path) -> np.ndarray: + """Returns shape (N, 3).""" + return _parse_vector_list(path.read_text()) + + +# --------------------------------------------------------------------------- +# polyMesh reader — compute cell centres and areas from scratch +# --------------------------------------------------------------------------- + +def _read_foam_points(text: str) -> np.ndarray: + """Parse polyMesh/points -> (N, 3) float array.""" + m = re.search(r'\n(\d+)\s*\n\((.+?)\n\)', text, re.S) + if not m: + raise ValueError("Could not parse points list") + tuples = re.findall(r'\(\s*([\d.eE+\-]+)\s+([\d.eE+\-]+)\s+([\d.eE+\-]+)\s*\)', m.group(2)) + return np.array([[float(a), float(b), float(c)] for a, b, c in tuples]) + + +def _read_foam_faces(text: str) -> list: + """Parse polyMesh/faces -> list of int-index lists.""" + m = re.search(r'\n(\d+)\s*\n\((.+?)\n\)', text, re.S) + if not m: + raise ValueError("Could not parse faces list") + faces = [] + for fm in re.finditer(r'\d+\(([^)]+)\)', m.group(2)): + faces.append(list(map(int, fm.group(1).split()))) + return faces + + +def _read_foam_labels(text: str) -> np.ndarray: + """Parse polyMesh/owner or neighbour -> 1-D int array.""" + m = re.search(r'\n(\d+)\s*\n\((.+?)\n\)', text, re.S) + if not m: + raise ValueError("Could not parse label list") + return np.array(list(map(int, m.group(2).split()))) + + +def compute_cell_centres_and_areas(mesh_dir: Path): + """ + Compute cell centres and (2-D) cell areas from polyMesh. + The mesh is 2-D extruded (z from 0 to some dz); we ignore z. + + Returns: + centres : (nCells, 2) x,y cell centres + areas : (nCells,) cell face areas projected to xy-plane + """ + pts = _read_foam_points((mesh_dir / "points").read_text()) # (nPoints, 3) + face_list = _read_foam_faces((mesh_dir / "faces").read_text()) # list of lists + owner = _read_foam_labels((mesh_dir / "owner").read_text()) # (nFaces,) + n_cells = int(owner.max()) + 1 + + # accumulate face-centre contributions per cell + sum_fc = np.zeros((n_cells, 2)) + count = np.zeros(n_cells, dtype=int) + + for fi, face in enumerate(face_list): + fc = pts[face, :2].mean(axis=0) # face centre (x,y) + c = owner[fi] + sum_fc[c] += fc + count[c] += 1 + + centres = sum_fc / count[:, None] + + # cell area: sum of triangle areas formed by face edges (shoelace on 2-D projection) + # For a structured annular mesh each cell is a quad; use face ownership to collect + # the 4 face-centre vertices and compute area via shoelace. + # Simpler: use face-centre average as above and compute area from cell vertex polygon. + # We build cell→vertex list then apply shoelace. + cell_verts = [set() for _ in range(n_cells)] + for fi, face in enumerate(face_list): + c = owner[fi] + cell_verts[c].update(face) + + # also need neighbour to add their faces + neighbour_text = (mesh_dir / "neighbour").read_text() + neighbour = _read_foam_labels(neighbour_text) + + for fi, face in enumerate(face_list): + if fi < len(neighbour): + c = neighbour[fi] + cell_verts[c].update(face) + + areas = np.zeros(n_cells) + for ci in range(n_cells): + vids = list(cell_verts[ci]) + vxy = pts[vids, :2] + # sort by angle around centroid so shoelace works + cx, cy = vxy.mean(axis=0) + angles = np.arctan2(vxy[:, 1] - cy, vxy[:, 0] - cx) + vxy = vxy[np.argsort(angles)] + x, y = vxy[:, 0], vxy[:, 1] + areas[ci] = 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1))) + + return centres, areas + + +# --------------------------------------------------------------------------- +# Analytical solution (Taylor-Couette) +# --------------------------------------------------------------------------- + +def analytical_v_theta(r, r1, r2, omega1, omega2=0.0): + A = (omega2 * r2**2 - omega1 * r1**2) / (r2**2 - r1**2) + B = (omega1 - omega2) * r1**2 * r2**2 / (r2**2 - r1**2) + return A * r + B / r + + +def analytical_velocity_xy(x, y, r1, r2, omega1, omega2=0.0): + """Returns (vx, vy) arrays at cell centres.""" + r = np.sqrt(x**2 + y**2) + theta = np.arctan2(y, x) + v_theta = analytical_v_theta(r, r1, r2, omega1, omega2) + vx = -v_theta * np.sin(theta) + vy = v_theta * np.cos(theta) + return vx, vy + + +def analytical_pressure_kinematic(x, y, r1, r2, omega1, omega2=0.0): + """Kinematic pressure p/rho (up to constant C).""" + r = np.sqrt(x**2 + y**2) + A = (omega2 * r2**2 - omega1 * r1**2) / (r2**2 - r1**2) + B = (omega1 - omega2) * r1**2 * r2**2 / (r2**2 - r1**2) + return A**2 * r**2 / 2.0 + 2*A*B * np.log(r) - B**2 / (2*r**2) + + +# --------------------------------------------------------------------------- +# L2 error (area-weighted, matching DuMux quadrature-over-elements) +# --------------------------------------------------------------------------- + +def l2_error_2d(sim, exact, areas): + """ + Area-weighted L2 error — equivalent to DuMux's element quadrature. + sim, exact : (N,) scalar or (N, dim) vector + areas : (N,) cell areas (integration weights) + """ + diff = sim - exact + if diff.ndim == 1: + err2 = np.sum(diff**2 * areas) + norm2 = np.sum(exact**2 * areas) + else: + err2 = np.sum(np.sum(diff**2, axis=1) * areas) + norm2 = np.sum(np.sum(exact**2, axis=1) * areas) + return np.sqrt(err2), np.sqrt(norm2) + + +# --------------------------------------------------------------------------- +# Per-config processing +# --------------------------------------------------------------------------- + +def process_config(case_dir: Path) -> dict | None: + params_file = case_dir / "parameters.json" + if not params_file.exists(): + print(f" [skip] no parameters.json in {case_dir.name}") + return None + + params = json.loads(params_file.read_text()) + r1 = params["domain"]["r1"] + r2 = params["domain"]["r2"] + omega1 = params["problem"]["omega1"] + omega2 = params["problem"].get("omega2", 0.0) + + # locate latest time directory + time_dirs = sorted( + [d for d in case_dir.iterdir() + if d.is_dir() and re.fullmatch(r'[\d.]+', d.name)], + key=lambda d: float(d.name) + ) + if not time_dirs: + print(f" [skip] no numeric time directories in {case_dir.name}") + return None + time_dir = time_dirs[-1] + print(f" Time directory: {time_dir.name}") + + # read fields + try: + U_all = read_vector_field(time_dir / "U") # (N, 3) + p_sim = read_scalar_field(time_dir / "p") # (N,) + except Exception as e: + print(f" [skip] field read error: {e}") + return None + + U_sim = U_all[:, :2] # drop z + + # compute cell centres and areas from polyMesh + mesh_dir = case_dir / "constant" / "polyMesh" + try: + centres, areas = compute_cell_centres_and_areas(mesh_dir) + except Exception as e: + print(f" [skip] mesh read error: {e}") + return None + + n = len(p_sim) + if len(centres) != n: + print(f" [skip] cell count mismatch: fields={n}, mesh={len(centres)}") + return None + + cx, cy = centres[:, 0], centres[:, 1] + + # --- analytical fields at cell centres --- + vx_ex, vy_ex = analytical_velocity_xy(cx, cy, r1, r2, omega1, omega2) + U_exact = np.stack([vx_ex, vy_ex], axis=1) + + p_exact = analytical_pressure_kinematic(cx, cy, r1, r2, omega1, omega2) + + # --- pressure offset correction (mirrors DuMux: subtract sim-exact at cell 0) --- + # DuMux subtracts (sim_p[0] - exact_p[0]) from the entire pressure vector + p_offset = p_sim[0] - p_exact[0] + p_sim_corrected = p_sim - p_offset + + # --- area-weighted L2 errors --- + l2_vel, norm_vel = l2_error_2d(U_sim, U_exact, areas) + l2_p, norm_p = l2_error_2d(p_sim_corrected, p_exact, areas) + + rel_vel = l2_vel / norm_vel if norm_vel > 1e-18 else l2_vel + rel_p = l2_p / norm_p if norm_p > 1e-18 else l2_p + + print(f" velocity L2 rel: {rel_vel:.4e} | pressure L2 rel: {rel_p:.4e}") + + metrics = { + "l2_error_velocity_rel": float(rel_vel), + "l2_error_pressure_rel": float(rel_p), + } + out_path = case_dir / "solution_metrics.json" + out_path.write_text(json.dumps(metrics, indent=2)) + print(f" Saved: {out_path.relative_to(case_dir.parent)}") + return metrics + +# --------------------------------------------------------------------------- +# Entry point +# --------------------------------------------------------------------------- + +def main(): + if len(sys.argv) > 1: + results_dir = Path(sys.argv[1]).resolve() + # else: + # results_dir = Path(__file__).resolve().parent / "results" + + if not results_dir.exists(): + print(f"Error: results directory not found: {results_dir}") + sys.exit(1) + + print(f"Scanning: {results_dir}\n") + + all_results = [] + for case_dir in sorted(results_dir.iterdir()): + if not case_dir.is_dir(): + continue + print(f"Processing: {case_dir.name}") + metrics = process_config(case_dir) + if metrics is None: + continue + params = json.loads((case_dir / "parameters.json").read_text()) + all_results.append({ + "conf": case_dir.name, + "cells_radial": params["grid"]["cells_radial"], + **metrics, + }) + + if not all_results: + print("\nNo valid configurations found.") + sys.exit(1) + + print(f"\nProcessed {len(all_results)} configuration(s).") + # plot_convergence(all_results, results_dir) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/benchmarks/rotating-cylinders/openfoam/constant/MRFProperties b/benchmarks/rotating-cylinders/openfoam/constant/MRFProperties deleted file mode 100644 index 18e2c3b..0000000 --- a/benchmarks/rotating-cylinders/openfoam/constant/MRFProperties +++ /dev/null @@ -1,30 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - location "constant"; - object MRFProperties; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -MRF1 -{ - cellZone all; - active yes; - - nonRotatingPatches (outerWall); - - origin (0 0 0); - axis (0 0 1); - omega {omega}; -} - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/constant/transportProperties b/benchmarks/rotating-cylinders/openfoam/constant/transportProperties deleted file mode 100644 index 2c6cb08..0000000 --- a/benchmarks/rotating-cylinders/openfoam/constant/transportProperties +++ /dev/null @@ -1,22 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - location "constant"; - object transportProperties; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -transportModel Newtonian; - -nu 1; - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/constant/turbulenceProperties b/benchmarks/rotating-cylinders/openfoam/constant/turbulenceProperties deleted file mode 100644 index 3ee1e9e..0000000 --- a/benchmarks/rotating-cylinders/openfoam/constant/turbulenceProperties +++ /dev/null @@ -1,20 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - location "constant"; - object turbulenceProperties; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -simulationType laminar; - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/output_template.zip b/benchmarks/rotating-cylinders/openfoam/output_template.zip new file mode 100644 index 0000000..3df79c1 Binary files /dev/null and b/benchmarks/rotating-cylinders/openfoam/output_template.zip differ diff --git a/benchmarks/rotating-cylinders/openfoam/plot b/benchmarks/rotating-cylinders/openfoam/plot deleted file mode 100755 index 8e001bd..0000000 --- a/benchmarks/rotating-cylinders/openfoam/plot +++ /dev/null @@ -1,46 +0,0 @@ -#!/bin/sh -cd "${0%/*}" || exit # Run from this directory -# Require gnuplot -command -v gnuplot >/dev/null || { - echo "FOAM FATAL ERROR: gnuplot not found - skipping graph creation" 1>&2 - exit 1 -} - - -gnuplot< np.ndarray: + m = re.search(r'internalField\s+nonuniform\s+List\s+\d+\s*\(([^)]+)\)', text, re.S) + if not m: + m = re.search(r'internalField\s+uniform\s+([\d.eE+\-]+)', text) + if m: + return np.array([float(m.group(1))]) + raise ValueError("Could not parse scalar internalField") + return np.fromstring(m.group(1), sep='\n') + + +def _parse_vector_list(text: str) -> np.ndarray: + m = re.search(r'internalField\s+nonuniform\s+List\s+\d+\s*\((.+?)\)\s*;', text, re.S) + if not m: + m2 = re.search(r'internalField\s+uniform\s+\(([\d.eE+\-\s]+)\)', text) + if m2: + return np.array([list(map(float, m2.group(1).split()))]) + raise ValueError("Could not parse vector internalField") + rows = re.findall(r'\(([\d.eE+\-\s]+)\)', m.group(1)) + return np.array([list(map(float, r.split())) for r in rows]) + + +def read_scalar_field(path: Path) -> np.ndarray: + return _parse_scalar_list(path.read_text()) + + +def read_vector_field(path: Path) -> np.ndarray: + return _parse_vector_list(path.read_text()) + + +# ============================================================================= +# polyMesh reader — cell centres and areas +# ============================================================================= + +def _read_foam_points(text: str) -> np.ndarray: + m = re.search(r'\n(\d+)\s*\n\((.+?)\n\)', text, re.S) + if not m: + raise ValueError("Could not parse points list") + tuples = re.findall(r'\(\s*([\d.eE+\-]+)\s+([\d.eE+\-]+)\s+([\d.eE+\-]+)\s*\)', m.group(2)) + return np.array([[float(a), float(b), float(c)] for a, b, c in tuples]) + + +def _read_foam_faces(text: str) -> list: + m = re.search(r'\n(\d+)\s*\n\((.+?)\n\)', text, re.S) + if not m: + raise ValueError("Could not parse faces list") + faces = [] + for fm in re.finditer(r'\d+\(([^)]+)\)', m.group(2)): + faces.append(list(map(int, fm.group(1).split()))) + return faces + + +def _read_foam_labels(text: str) -> np.ndarray: + m = re.search(r'\n(\d+)\s*\n\((.+?)\n\)', text, re.S) + if not m: + raise ValueError("Could not parse label list") + return np.array(list(map(int, m.group(2).split()))) + + +def compute_cell_centres_and_areas(mesh_dir: Path): + pts = _read_foam_points((mesh_dir / "points").read_text()) + face_list = _read_foam_faces((mesh_dir / "faces").read_text()) + owner = _read_foam_labels((mesh_dir / "owner").read_text()) + n_cells = int(owner.max()) + 1 + + sum_fc = np.zeros((n_cells, 2)) + count = np.zeros(n_cells, dtype=int) + for fi, face in enumerate(face_list): + fc = pts[face, :2].mean(axis=0) + sum_fc[owner[fi]] += fc + count[owner[fi]] += 1 + centres = sum_fc / count[:, None] + + cell_verts = [set() for _ in range(n_cells)] + for fi, face in enumerate(face_list): + cell_verts[owner[fi]].update(face) + neighbour = _read_foam_labels((mesh_dir / "neighbour").read_text()) + for fi, face in enumerate(face_list): + if fi < len(neighbour): + cell_verts[neighbour[fi]].update(face) + + areas = np.zeros(n_cells) + for ci in range(n_cells): + vxy = pts[list(cell_verts[ci]), :2] + cx, cy = vxy.mean(axis=0) + angles = np.arctan2(vxy[:, 1] - cy, vxy[:, 0] - cx) + vxy = vxy[np.argsort(angles)] + x, y = vxy[:, 0], vxy[:, 1] + areas[ci] = 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1))) + + return centres, areas + + +# ============================================================================= +# Analytical solution (Taylor-Couette) +# ============================================================================= + +def analytical_solution(r, omega1, omega2, r1, r2): + A = (omega2 * r2**2 - omega1 * r1**2) / (r2**2 - r1**2) + B = (omega1 - omega2) * r1**2 * r2**2 / (r2**2 - r1**2) + u = A * r + B / r + p = A**2 * r**2 / 2.0 + 2*A*B * np.log(r) - B**2 / (2*r**2) + p -= p[0] + return u, p + + +def analytical_velocity_xy(x, y, omega1, omega2, r1, r2): + r = np.sqrt(x**2 + y**2) + theta = np.arctan2(y, x) + A = (omega2 * r2**2 - omega1 * r1**2) / (r2**2 - r1**2) + B = (omega1 - omega2) * r1**2 * r2**2 / (r2**2 - r1**2) + v_theta = A * r + B / r + return -v_theta * np.sin(theta), v_theta * np.cos(theta) + + +def analytical_pressure_kinematic(x, y, omega1, omega2, r1, r2): + r = np.sqrt(x**2 + y**2) + A = (omega2 * r2**2 - omega1 * r1**2) / (r2**2 - r1**2) + B = (omega1 - omega2) * r1**2 * r2**2 / (r2**2 - r1**2) + return A**2 * r**2 / 2.0 + 2*A*B * np.log(r) - B**2 / (2*r**2) + + +# ============================================================================= +# Radial profile binning (mirrors DuMux compute_radial_profiles) +# ============================================================================= + +def compute_radial_profiles(cx, cy, pressure, velocity, r1, r2, n_bins=80): + vx, vy = velocity[:, 0], velocity[:, 1] + speed = np.sqrt(vx**2 + vy**2) + r_cells = np.sqrt(cx**2 + cy**2) + + mask = (r_cells >= r1) & (r_cells <= r2) + r_f, p_f, s_f = r_cells[mask], pressure[mask], speed[mask] + + bins = np.linspace(r1, r2, n_bins + 1) + bc = 0.5 * (bins[:-1] + bins[1:]) + idx = np.clip(np.digitize(r_f, bins) - 1, 0, n_bins - 1) + + p_bin = np.array([p_f[idx == i].mean() if (idx == i).any() else np.nan for i in range(n_bins)]) + s_bin = np.array([s_f[idx == i].mean() if (idx == i).any() else np.nan for i in range(n_bins)]) + + if not np.isnan(p_bin).all(): + p_bin -= p_bin[~np.isnan(p_bin)][0] + + return bc, s_bin, p_bin + + +# ============================================================================= +# Case discovery +# ============================================================================= + +def discover_cases(results_dir: Path): + cases = [] + for d in sorted(results_dir.iterdir()): + if not d.is_dir(): + continue + params_file = d / "parameters.json" + if not params_file.exists(): + continue + # find latest time directory — may be absent if outputs were zipped + time_dirs = sorted( + [t for t in d.iterdir() if t.is_dir() and re.fullmatch(r'[\d.]+', t.name)], + key=lambda t: float(t.name) + ) + cases.append({ + "path": d, + "params": json.loads(params_file.read_text()), + "time_dir": time_dirs[-1] if time_dirs else None, + "mesh_dir": d / "constant" / "polyMesh", + }) + return cases + + +# ============================================================================= +# Per-case plots +# ============================================================================= + +def plot_profiles(case: dict): + """Radial profiles — tangential velocity and pressure vs analytical.""" + params = case["params"] + r1 = params["domain"]["r1"] + r2 = params["domain"]["r2"] + omega1 = params["problem"]["omega1"] + omega2 = params["problem"].get("omega2", 0.0) + + try: + U_all = read_vector_field(case["time_dir"] / "U") + p_sim = read_scalar_field(case["time_dir"] / "p") + centres, _ = compute_cell_centres_and_areas(case["mesh_dir"]) + except Exception as e: + print(f" [skip profiles] {e}") + return + + cx, cy = centres[:, 0], centres[:, 1] + velocity = U_all[:, :2] + + # pressure offset correction + p_exact0 = analytical_pressure_kinematic(cx, cy, omega1, omega2, r1, r2) + p_sim = p_sim - (p_sim[0] - p_exact0[0]) + + r_num, u_num, p_num = compute_radial_profiles(cx, cy, p_sim, velocity, r1, r2) + u_ref, p_ref = analytical_solution(r_num, omega1, omega2, r1, r2) + + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=FIG_SIZE, dpi=DPI) + + ax1.plot(r_num, u_ref, lw=2, color="black", label="Analytical") + ax1.plot(r_num, u_num, "o", ms=4, mfc="none", color="red", label="OpenFOAM") + ax1.set_title("Tangential Velocity Profile") + ax1.set_xlabel("r") + ax1.set_ylabel("|u|") + ax1.legend() + ax1.tick_params(which="both", direction="in") + + ax2.plot(r_num, p_ref, lw=2, color="black", label="Analytical") + ax2.plot(r_num, p_num, "o", ms=4, mfc="none", color="red", label="OpenFOAM") + ax2.set_title("Pressure Profile") + ax2.set_xlabel("r") + ax2.set_ylabel("p") + ax2.legend() + ax2.tick_params(which="both", direction="in") + + plt.suptitle(case["path"].name, fontsize=11) + plt.tight_layout() + out = case["path"] / "profiles.png" + fig.savefig(out, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {out.name}") + + +def plot_fields(case: dict): + """2D colour maps for pressure and velocity magnitude over the annular domain.""" + params = case["params"] + r1 = params["domain"]["r1"] + r2 = params["domain"]["r2"] + omega1 = params["problem"]["omega1"] + omega2 = params["problem"].get("omega2", 0.0) + + try: + U_all = read_vector_field(case["time_dir"] / "U") + p_sim = read_scalar_field(case["time_dir"] / "p") + centres, _ = compute_cell_centres_and_areas(case["mesh_dir"]) + except Exception as e: + print(f" [skip fields] {e}") + return + + cx, cy = centres[:, 0], centres[:, 1] + + # pressure offset correction + p_exact0 = analytical_pressure_kinematic(cx, cy, omega1, omega2, r1, r2) + p_sim = p_sim - (p_sim[0] - p_exact0[0]) + + speed = np.sqrt(U_all[:, 0]**2 + U_all[:, 1]**2) + + # Triangulate cell centres, then mask any triangle whose centroid + # falls inside the inner cylinder (r < r1) — these are spurious + # triangles that the Delaunay algorithm draws across the hole. + triang = mtri.Triangulation(cx, cy) + tri_cx = cx[triang.triangles].mean(axis=1) + tri_cy = cy[triang.triangles].mean(axis=1) + tri_r = np.sqrt(tri_cx**2 + tri_cy**2) + triang.set_mask(tri_r < r1) + + theta = np.linspace(0, 2*np.pi, 360) + + for field_vals, label, cmap, fname in [ + (p_sim, "p", "coolwarm", "pressure_field.png"), + (speed, "|U|", "viridis", "velocity_field.png"), + ]: + fig, ax = plt.subplots(figsize=(7, 6), dpi=FIELD_DPI) + tc = ax.tricontourf(triang, field_vals, levels=64, cmap=cmap) + fig.colorbar(tc, ax=ax, label=label, fraction=0.046, pad=0.04) + + # cylinder boundary lines + ax.plot(r1*np.cos(theta), r1*np.sin(theta), "k-", lw=1.2) + ax.plot(r2*np.cos(theta), r2*np.sin(theta), "k-", lw=1.2) + + # fill the inner-cylinder hole with a solid white patch so no + # colour bleeds inside even at boundaries + from matplotlib.patches import Circle + ax.add_patch(Circle((0, 0), r1, color="white", zorder=3)) + ax.plot(r1*np.cos(theta), r1*np.sin(theta), "k-", lw=1.2, zorder=4) + + ax.set_aspect("equal") + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_title(f"{label} — {case['path'].name}") + # axis ticks and labels only — no grid boxes + ax.tick_params(which="both", direction="in") + plt.tight_layout() + + out = case["path"] / fname + fig.savefig(out, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {out.name}") + + +# ============================================================================= +# Global convergence plot +# ============================================================================= + +def plot_convergence(results: list, results_dir: Path): + try: + results.sort(key=lambda d: int(d["cells_radial"])) + except Exception: + results.sort(key=lambda d: d["conf"]) + + confs = [d["conf"] for d in results] + p_errors = [d["l2_error_pressure_rel"] for d in results] + v_errors = [d["l2_error_velocity_rel"] for d in results] + + fig, ax = plt.subplots(figsize=(10, 6), dpi=DPI) + ax.plot(confs, p_errors, marker="o", linewidth=1.5, + label=r"Relative $L^2$ Error — Pressure") + ax.plot(confs, v_errors, marker="s", linewidth=1.5, + label=r"Relative $L^2$ Error — Velocity") + + ax.set_yscale("log") + ax.set_xlabel("Configuration") + ax.set_ylabel(r"Relative $L^2$ Error") + ax.set_title("OpenFOAM Rotating Cylinders — Convergence Summary") + ax.legend() + ax.grid(True, which="both", linestyle="-", alpha=0.3) + plt.xticks(rotation=45) + plt.tight_layout() + + out = results_dir / "solution_metrics_plot.png" + fig.savefig(out, dpi=DPI) + plt.close(fig) + print(f"\nConvergence plot saved to: {out}") + + + +# ============================================================================= +# Load raw field data — from disk or by unzipping on the fly +# ============================================================================= + +def load_case_fields(case: dict) -> dict | None: + """ + Return a dict with time_dir, mesh_dir pointing at real or temp-extracted + paths, plus a cleanup() callable. Returns None if data cannot be found. + + The caller must call result["cleanup"]() when done to remove any temp dir. + """ + # Raw files already on disk + if case["time_dir"] is not None: + return { + "time_dir": case["time_dir"], + "mesh_dir": case["mesh_dir"], + "cleanup": lambda: None, + } + + # Try unzipping + zip_files = sorted(case["path"].glob("*.zip")) + if not zip_files: + return None + + tmp = tempfile.mkdtemp() + tmp_path = Path(tmp) + try: + with zipfile.ZipFile(zip_files[0], "r") as zf: + zf.extractall(tmp_path) + + time_dirs = sorted( + [t for t in tmp_path.iterdir() + if t.is_dir() and re.fullmatch(r"[\d.]+", t.name)], + key=lambda t: float(t.name), + ) + mesh_candidates = list((tmp_path / "constant").glob("polyMesh")) + + if not time_dirs or not mesh_candidates: + shutil.rmtree(tmp, ignore_errors=True) + return None + + return { + "time_dir": time_dirs[-1], + "mesh_dir": mesh_candidates[0], + "cleanup": lambda: shutil.rmtree(tmp, ignore_errors=True), + } + except Exception: + shutil.rmtree(tmp, ignore_errors=True) + return None + + +# ============================================================================= +# Analytical comparison plot +# ============================================================================= + +def plot_analytical_comparison(case: dict, fields: dict, results_dir: Path): + """ + Side-by-side radial profiles comparing OpenFOAM results with the + analytical Taylor-Couette solution. + Saved to /analytical_comparison.png + """ + params = case["params"] + r1 = params["domain"]["r1"] + r2 = params["domain"]["r2"] + omega1 = params["problem"]["omega1"] + omega2 = params["problem"].get("omega2", 0.0) + + try: + U_all = read_vector_field(fields["time_dir"] / "U") + p_sim = read_scalar_field(fields["time_dir"] / "p") + centres, _ = compute_cell_centres_and_areas(fields["mesh_dir"]) + except Exception as e: + print(f" [skip analytical comparison] {e}") + return + + cx, cy = centres[:, 0], centres[:, 1] + velocity = U_all[:, :2] + + # Align pressure reference to analytical at first cell + p_exact0 = analytical_pressure_kinematic(cx, cy, omega1, omega2, r1, r2) + p_sim = p_sim - (p_sim[0] - p_exact0[0]) + + # Bin numerical data + r_num, u_num, p_num = compute_radial_profiles(cx, cy, p_sim, velocity, r1, r2) + + # Dense analytical curves + r_ref = np.linspace(r1, r2, 400) + u_ref, p_ref = analytical_solution(r_ref, omega1, omega2, r1, r2) + + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5), dpi=DPI) + + # Velocity + ax1.plot(r_ref, u_ref, lw=2, color="black", label="Analytical") + ax1.plot(r_num, u_num, "o", ms=4, mfc="none", color="tab:red", + label=f"OpenFOAM") + ax1.set_xlabel("r", fontsize=12) + ax1.set_ylabel(r"$u_\theta$", fontsize=12) + ax1.set_title("Velocity", fontsize=12) + ax1.legend(fontsize=11) + ax1.grid(True, linestyle="--", alpha=0.4) + ax1.tick_params(which="both", direction="in") + + # Pressure + ax2.plot(r_ref, p_ref, lw=2, color="black", label="Analytical") + ax2.plot(r_num, p_num, "s", ms=4, mfc="none", color="tab:blue", + label=f"OpenFOAM") + ax2.set_xlabel("r", fontsize=12) + ax2.set_ylabel("p", fontsize=12) + ax2.set_title("Pressure", fontsize=12) + ax2.legend(fontsize=11) + ax2.grid(True, linestyle="--", alpha=0.4) + ax2.tick_params(which="both", direction="in") + + fig.suptitle( + f"OpenFOAM vs Analytical — {case['path'].name}", + fontsize=12, fontweight="bold", + ) + plt.tight_layout() + + out = case["path"] / "analytical_comparison.png" + fig.savefig(out, dpi=DPI, bbox_inches="tight") + plt.close(fig) + print(f" Saved: {out}") + + +# ============================================================================= +# Entry point +# ============================================================================= + +def main(): + parser = argparse.ArgumentParser( + description="Post-process OpenFOAM rotating-cylinders benchmark results." + ) + parser.add_argument( + "results_dir", nargs="?", default=None, + help="Path to results directory (default: ./results/)", + ) + parser.add_argument( + "--conf", metavar="CONF_NAME", default=None, + help=( + "Name of a single configuration to process, e.g. 100_0_80_320. " + "When omitted, all configurations are processed for the convergence plot." + ), + ) + parser.add_argument( + "--fields", action="store_true", + help="Generate 2D colour field maps (pressure_field.png, velocity_field.png). " + "Requires --conf.", + ) + parser.add_argument( + "--analytical", action="store_true", + help="Generate analytical comparison plot (analytical_comparison.png). " + "Requires --conf.", + ) + args = parser.parse_args() + + # Validate flag combinations + if (args.fields or args.analytical) and not args.conf: + parser.error("--fields and --analytical require --conf ") + + results_dir = Path(args.results_dir).resolve() if args.results_dir else RESULTS_DIR + if not results_dir.exists(): + print(f"Error: results directory not found: {results_dir}") + sys.exit(1) + + print(f"Scanning: {results_dir}\n") + cases = discover_cases(results_dir) + if not cases: + print("No valid cases found.") + sys.exit(1) + + # ── Single-configuration mode ───────────────────────────────────────────── + if args.conf: + matched = [c for c in cases if c["path"].name == args.conf] + if not matched: + available = " \n".join(c["path"].name for c in cases) + print(f"Error: configuration '{args.conf}' not found.\n" + f"Available configurations:\n {available}") + sys.exit(1) + + case = matched[0] + print("=" * 60) + print(f"Configuration: {case['path'].name}") + + if args.fields or args.analytical: + print(f" Loading field data ...") + fields = load_case_fields(case) + if fields is None: + print(" [error] Could not load field data — no time directory or zip found.") + sys.exit(1) + try: + if args.fields: + print(" Rendering 2D field maps...") + # plot_fields expects case["time_dir"] / case["mesh_dir"] at top level + plot_fields({**case, "time_dir": fields["time_dir"], + "mesh_dir": fields["mesh_dir"]}) + if args.analytical: + print(" Generating analytical comparison plot...") + plot_analytical_comparison(case, fields, results_dir) + finally: + fields["cleanup"]() + else: + # No extra flags — just run profiles (needs field data too) + print(" Loading field data for profiles...") + fields = load_case_fields(case) + if fields is None: + print(" [skip profiles] no field data available") + else: + try: + plot_profiles({**case, "time_dir": fields["time_dir"], + "mesh_dir": fields["mesh_dir"]}) + finally: + fields["cleanup"]() + + print("=" * 60) + print("Done.") + return + + # ── All-configurations mode (convergence plot) ──────────────────────────── + metrics_data = [] + + for case in cases: + print("=" * 60) + print(f"Processing: {case['path'].name}") + + # Profiles where field data is available + fields = load_case_fields(case) + if fields is not None: + try: + plot_profiles({**case, "time_dir": fields["time_dir"], + "mesh_dir": fields["mesh_dir"]}) + finally: + fields["cleanup"]() + else: + print(" [skip profiles] no field data available") + + # Collect metrics + metrics_file = case["path"] / "solution_metrics.json" + if metrics_file.exists(): + try: + metrics = json.loads(metrics_file.read_text()) + metrics_data.append({ + "conf": case["path"].name, + "cells_radial": case["params"]["grid"]["cells_radial"], + **metrics, + }) + except Exception as e: + print(f" Warning: could not read solution_metrics.json: {e}") + else: + print(f" Warning: no solution_metrics.json in {case['path'].name}") + + if metrics_data: + plot_convergence(metrics_data, results_dir) + else: + print("\nNo metrics found — skipping convergence plot.") + + print("=" * 60) + print("All processing complete.") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/benchmarks/rotating-cylinders/openfoam/run_benchmark.py b/benchmarks/rotating-cylinders/openfoam/run_benchmark.py index 8e90bf9..0740ea7 100644 --- a/benchmarks/rotating-cylinders/openfoam/run_benchmark.py +++ b/benchmarks/rotating-cylinders/openfoam/run_benchmark.py @@ -1,5 +1,4 @@ import json -import shutil import subprocess import zipfile import os @@ -7,20 +6,25 @@ root_dir = Path(__file__).resolve().parent zip_path = root_dir.parent / "rotating-cylinders.zip" -benchmark_dir = root_dir / "rotating-cylinders" +template_zip_path = root_dir / "output_template.zip" snakefile_path = root_dir / "Snakefile" # Extraction if zip_path.exists(): with zipfile.ZipFile(zip_path, 'r') as zip_ref: - zip_ref.extractall(benchmark_dir) - print(f"Successfully extracted benchmark to: {benchmark_dir}") + zip_ref.extractall(root_dir) + print(f"Successfully extracted benchmark to: {root_dir}") else: print(f"Error: Could not find {zip_path} at {root_dir.parent}") exit(1) +# Validate template zip exists +if not template_zip_path.exists(): + print(f"Error: Could not find output_template.zip at {template_zip_path}") + exit(1) + # Iterate through all parameter files -for param_file in benchmark_dir.glob("parameters_*.json"): +for param_file in root_dir.glob("parameters_*.json"): with open(param_file, "r") as f: data = json.load(f) config_name = data.get("configuration") @@ -37,24 +41,10 @@ with open(output_dir / "parameters.json", "w") as outfile: json.dump(data, outfile, indent=2) - # Copy OpenFOAM template files to the output directory - for item in ["0", "system", "constant", "Allrun", "Allclean", "plot"]: - src = root_dir / item - dst = output_dir / item - if src.is_dir(): - if dst.exists(): - shutil.rmtree(dst) - shutil.copytree(src, dst) - elif src.is_file(): - shutil.copy(src, dst) - - # Copy any additional files from the zip extract dir (excluding parameter files) - for item in benchmark_dir.iterdir(): - if item.is_file(): - if item.name.startswith("parameters_") and item.suffix == ".json": - continue - else: - shutil.copy(item, output_dir / item.name) + # Extract OpenFOAM template files from output_template.zip into the output directory + with zipfile.ZipFile(template_zip_path, 'r') as zip_ref: + zip_ref.extractall(output_dir) + print(f"Extracted output_template.zip to: {output_dir}") # Run the Snakemake workflow for the configuration subprocess.run([ @@ -67,12 +57,29 @@ ], check=True, cwd=output_dir) print(f"Workflow executed successfully for {config_name}.") + # Zip all output files except solution_metrics.json, then remove the originals + output_zip_path = output_dir / f"{config_name}_files.zip" + with zipfile.ZipFile(output_zip_path, 'w', zipfile.ZIP_DEFLATED) as zf: + for item in output_dir.rglob("*"): + if item.is_file() and item.name not in {"solution_metrics.json", "parameters.json"} and item != output_zip_path: + zf.write(item, item.relative_to(output_dir)) + item.unlink() + + # Remove any empty subdirectories left after zipping + for item in sorted(output_dir.rglob("*"), reverse=True): + if item.is_dir(): + try: + item.rmdir() + except OSError: + pass + print(f"Outputs zipped to: {output_zip_path.name}") + print("\nAll configurations processed.") # --- CLEANUP SECTION --- print("\nStarting cleanup...") -for param_file in benchmark_dir.glob("parameters_*.json"): +for param_file in root_dir.glob("parameters_*.json"): try: os.remove(param_file) print(f"Deleted: {param_file.name}") diff --git a/benchmarks/rotating-cylinders/openfoam/system/blockMeshDict.template b/benchmarks/rotating-cylinders/openfoam/system/blockMeshDict.template deleted file mode 100644 index 9a82a6c..0000000 --- a/benchmarks/rotating-cylinders/openfoam/system/blockMeshDict.template +++ /dev/null @@ -1,130 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - object blockMeshDict; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -scale 1; - -geom -{ - r0 {r0}; - r1 {r1}; - - nr {nr}; - ntheta {ntheta}; - - g {g}; - - mr0 #eval{ - $r0 }; - mr1 #eval{ - $r1 }; - invG #eval{ 1./$g }; -} - -vertices -( - ( $/geom/r0 0 0) - ( 0 $/geom/r0 0) - ($/geom/mr0 0 0) - ( 0 $/geom/mr0 0) - ( $/geom/r1 0 0) - ( 0 $/geom/r1 0) - ($/geom/mr1 0 0) - ( 0 $/geom/mr1 0) - - ( $/geom/r0 0 1) - ( 0 $/geom/r0 1) - ($/geom/mr0 0 1) - ( 0 $/geom/mr0 1) - ( $/geom/r1 0 1) - ( 0 $/geom/r1 1) - ($/geom/mr1 0 1) - ( 0 $/geom/mr1 1) -); - -blockInfo -all -($/geom/ntheta $/geom/nr 1) -simpleGrading (1 ((0.5 0.5 $/geom/g)(0.5 0.5 $/geom/invG)) 1); - -blocks -( - hex (1 0 4 5 9 8 12 13) $blockInfo - hex (2 1 5 6 10 9 13 14) $blockInfo - hex (3 2 6 7 11 10 14 15) $blockInfo - hex (0 3 7 4 8 11 15 12) $blockInfo -); - -edges -( - arc 0 1 origin (0 0 0) - arc 1 2 origin (0 0 0) - arc 2 3 origin (0 0 0) - arc 3 0 origin (0 0 0) - arc 8 9 origin (0 0 1) - arc 9 10 origin (0 0 1) - arc 10 11 origin (0 0 1) - arc 11 8 origin (0 0 1) - - arc 4 5 origin (0 0 0) - arc 5 6 origin (0 0 0) - arc 6 7 origin (0 0 0) - arc 7 4 origin (0 0 0) - arc 12 13 origin (0 0 1) - arc 13 14 origin (0 0 1) - arc 14 15 origin (0 0 1) - arc 15 12 origin (0 0 1) -); - -defaultPatch -{ - name frontAndBack; - type empty; -} - -boundary -( - innerWall - { - type wall; - faces - ( - (1 0 8 9) - (2 1 9 10) - (3 2 10 11) - (0 3 11 8) - ); - } - - outerWall - { - type wall; - faces - ( - (5 13 12 4) - (6 14 13 5) - (7 15 14 6) - (4 12 15 7) - ); - } -); - -mergePatchPairs -( -); - - -// Cleanup -#remove ( geom blockInfo ) - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/system/controlDict b/benchmarks/rotating-cylinders/openfoam/system/controlDict deleted file mode 100644 index a40242d..0000000 --- a/benchmarks/rotating-cylinders/openfoam/system/controlDict +++ /dev/null @@ -1,74 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - location "system"; - object controlDict; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -application simpleFoam; - -startFrom startTime; - -startTime 0; - -stopAt endTime; - -endTime 10000; - -deltaT 1; - -writeControl timeStep; - -writeInterval 10000; - -purgeWrite 0; - -writeFormat ascii; - -writePrecision 6; - -writeCompression off; - -timeFormat general; - -timePrecision 6; - -runTimeModifiable true; - -functions -{ - sample1 - { - type sets; - libs (sampling); - writeControl writeTime; - fields (U p); - interpolationScheme cellPoint; - setFormat raw; - - sets - ( - centreLine - { - type uniform; - axis x; - start (1 0 0); - end (2 0 0); - nPoints 20; - } - ); - } -} - - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/system/fvSchemes b/benchmarks/rotating-cylinders/openfoam/system/fvSchemes deleted file mode 100644 index 0a0bd18..0000000 --- a/benchmarks/rotating-cylinders/openfoam/system/fvSchemes +++ /dev/null @@ -1,51 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - location "system"; - object fvSchemes; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -ddtSchemes -{ - default steadyState; -} - -gradSchemes -{ - default Gauss linear; -} - -divSchemes -{ - default none; - div(phi,U) bounded Gauss linearUpwind grad(U); - div((nuEff*dev2(T(grad(U))))) Gauss linear; -} - -laplacianSchemes -{ - default Gauss linear corrected; -} - -interpolationSchemes -{ - default linear; -} - -snGradSchemes -{ - default corrected; -} - - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/openfoam/system/fvSolution b/benchmarks/rotating-cylinders/openfoam/system/fvSolution deleted file mode 100644 index f137672..0000000 --- a/benchmarks/rotating-cylinders/openfoam/system/fvSolution +++ /dev/null @@ -1,57 +0,0 @@ -/*--------------------------------*- C++ -*----------------------------------*\ -| ========= | | -| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: v2412 | -| \\ / A nd | Website: www.openfoam.com | -| \\/ M anipulation | | -\*---------------------------------------------------------------------------*/ -FoamFile -{ - version 2.0; - format ascii; - class dictionary; - location "system"; - object fvSolution; -} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -solvers -{ - p - { - solver GAMG; - tolerance 1e-10; - relTol 0.1; - smoother GaussSeidel; - } - - U - { - solver smoothSolver; - smoother GaussSeidel; - tolerance 1e-10; - relTol 0.1; - } -} - -SIMPLE -{ - nNonOrthogonalCorrectors 0; - pRefCell 0; - pRefValue 0; -} - -relaxationFactors -{ - fields - { - p 0.3; - } - equations - { - U 0.7; - } -} - - -// ************************************************************************* // diff --git a/benchmarks/rotating-cylinders/plot_convergence.py b/benchmarks/rotating-cylinders/plot_convergence.py new file mode 100644 index 0000000..b781209 --- /dev/null +++ b/benchmarks/rotating-cylinders/plot_convergence.py @@ -0,0 +1,162 @@ +""" +plot_convergence.py +------------------- +General convergence plotting script for the rotating-cylinders benchmark. + +Place this file at: + benchmarks/rotating-cylinders/plot_convergence.py + +It auto-discovers every solver subfolder that contains a results/ directory, +reads solution_metrics.json and parameters.json from each case, and produces: + + 1. One plot per solver → /results/solution_metrics_plot.png + 2. One combined plot → rotating-cylinders/convergence_comparison.png +""" + +import json +import matplotlib.pyplot as plt +from pathlib import Path + +# --------------------------------------------------------------------------- +# Paths +# --------------------------------------------------------------------------- +root_dir = Path(__file__).resolve().parent # rotating-cylinders/ +results_key = "results" # expected subfolder name + +# --------------------------------------------------------------------------- +# Discover solvers: any direct subdirectory that contains a results/ folder +# --------------------------------------------------------------------------- +solvers = sorted( + d for d in root_dir.iterdir() + if d.is_dir() and (d / results_key).is_dir() +) + +if not solvers: + print("No solver directories with a 'results/' folder found. Exiting.") + exit(1) + +print(f"Found solvers: {[s.name for s in solvers]}\n") + +# --------------------------------------------------------------------------- +# Helper: load all cases for one solver +# --------------------------------------------------------------------------- +def load_solver_results(solver_dir: Path) -> list[dict]: + """Return a list of result dicts sorted by cells_radial (or conf name).""" + results_dir = solver_dir / results_key + all_results = [] + + for case_dir in results_dir.iterdir(): + if not case_dir.is_dir(): + continue + metrics_file = case_dir / "solution_metrics.json" + params_file = case_dir / "parameters.json" + + if not metrics_file.exists(): + print(f" [skip] no solution_metrics.json in {solver_dir.name}/{case_dir.name}") + continue + + metrics = json.loads(metrics_file.read_text()) + params = json.loads(params_file.read_text()) if params_file.exists() else {} + + all_results.append({ + "conf": case_dir.name, + "cells_radial": params.get("grid", {}).get("cells_radial"), + "pressure_error": metrics.get("l2_error_pressure_rel"), + "velocity_error": metrics.get("l2_error_velocity_rel"), + }) + + if not all_results: + return all_results + + # Sort by cells_radial when available, otherwise by conf name + if all(d["cells_radial"] is not None for d in all_results): + all_results.sort(key=lambda d: d["cells_radial"]) + else: + all_results.sort(key=lambda d: d["conf"]) + + return all_results + +# --------------------------------------------------------------------------- +# Helper: draw a single-solver convergence plot +# --------------------------------------------------------------------------- +def plot_single_solver(solver_name: str, data: list[dict], save_path: Path) -> None: + confs = [d["conf"] for d in data] + p_errors = [d["pressure_error"] for d in data] + v_errors = [d["velocity_error"] for d in data] + + fig, ax = plt.subplots(figsize=(10, 6)) + ax.plot(confs, p_errors, marker="o", linewidth=1.5, + label=r"Relative $L^2$ Error — Pressure") + ax.plot(confs, v_errors, marker="s", linewidth=1.5, + label=r"Relative $L^2$ Error — Velocity") + + ax.set_yscale("log") + ax.set_xlabel("Configuration") + ax.set_ylabel(r"Relative $L^2$ Error") + ax.set_title(f"{solver_name.capitalize()} Rotating Cylinders — Convergence Summary") + ax.legend() + ax.grid(True, which="both", linestyle="-", alpha=0.3) + plt.xticks(rotation=45) + plt.tight_layout() + + fig.savefig(save_path, dpi=150) + plt.close(fig) + print(f" Saved: {save_path.relative_to(root_dir)}") + +# --------------------------------------------------------------------------- +# 1. Per-solver plots +# --------------------------------------------------------------------------- +print("--- Per-solver plots ---") +solver_data = {} # name → list[dict] (kept for the combined plot) + +for solver_dir in solvers: + data = load_solver_results(solver_dir) + if not data: + print(f" [skip] no valid cases for {solver_dir.name}") + continue + + solver_data[solver_dir.name] = data + save_path = solver_dir / results_key / "solution_metrics_plot.png" + plot_single_solver(solver_dir.name, data, save_path) + +# --------------------------------------------------------------------------- +# 2. Combined comparison plot +# --------------------------------------------------------------------------- +print("\n--- Combined comparison plot ---") + +if len(solver_data) < 2: + print(" Only one solver with data — skipping combined plot.") +else: + # Use a distinct marker per solver so lines stay distinguishable in B&W too + markers = ["o", "s", "^", "D", "v", "P", "X", "*"] + + fig, axes = plt.subplots(1, 2, figsize=(16, 6), sharey=False) + + for ax, error_key, field in zip( + axes, + ["pressure_error", "velocity_error"], + ["Pressure", "Velocity"], + ): + for idx, (solver_name, data) in enumerate(solver_data.items()): + confs = [d["conf"] for d in data] + errors = [d[error_key] for d in data] + marker = markers[idx % len(markers)] + ax.plot(confs, errors, marker=marker, linewidth=1.5, label=solver_name.capitalize()) + + ax.set_yscale("log") + ax.set_xlabel("Configuration") + ax.set_ylabel(r"Relative $L^2$ Error") + ax.set_title(f"Convergence Comparison — {field}") + ax.legend() + ax.grid(True, which="both", linestyle="-", alpha=0.3) + plt.setp(ax.get_xticklabels(), rotation=45) + + fig.suptitle("Rotating Cylinders — Solver Convergence Comparison", fontsize=14) + plt.tight_layout() + + combined_path = root_dir / "convergence_comparison.png" + fig.savefig(combined_path, dpi=150) + plt.close(fig) + print(f" Saved: {combined_path.relative_to(root_dir)}") + +print("\nDone.") diff --git a/docs/benchmarks/rotating_cylinders/convergence_comparison.png b/docs/benchmarks/rotating_cylinders/convergence_comparison.png new file mode 100644 index 0000000..bc84093 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/convergence_comparison.png differ diff --git a/docs/benchmarks/rotating_cylinders/dumux_plots/analytical_comparison.png b/docs/benchmarks/rotating_cylinders/dumux_plots/analytical_comparison.png new file mode 100644 index 0000000..8c7ae84 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/dumux_plots/analytical_comparison.png differ diff --git a/docs/benchmarks/rotating_cylinders/dumux_plots/pressure_field.png b/docs/benchmarks/rotating_cylinders/dumux_plots/pressure_field.png new file mode 100644 index 0000000..038b002 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/dumux_plots/pressure_field.png differ diff --git a/docs/benchmarks/rotating_cylinders/dumux_plots/solution_metrics_plot.png b/docs/benchmarks/rotating_cylinders/dumux_plots/solution_metrics_plot.png new file mode 100644 index 0000000..21f9339 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/dumux_plots/solution_metrics_plot.png differ diff --git a/docs/benchmarks/rotating_cylinders/dumux_plots/velocity_field.png b/docs/benchmarks/rotating_cylinders/dumux_plots/velocity_field.png new file mode 100644 index 0000000..8eea08e Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/dumux_plots/velocity_field.png differ diff --git a/docs/benchmarks/rotating_cylinders/openfoam_plots/analytical_comparison.png b/docs/benchmarks/rotating_cylinders/openfoam_plots/analytical_comparison.png new file mode 100644 index 0000000..6d775e6 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/openfoam_plots/analytical_comparison.png differ diff --git a/docs/benchmarks/rotating_cylinders/openfoam_plots/pressure_field.png b/docs/benchmarks/rotating_cylinders/openfoam_plots/pressure_field.png new file mode 100644 index 0000000..623190e Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/openfoam_plots/pressure_field.png differ diff --git a/docs/benchmarks/rotating_cylinders/openfoam_plots/solution_metrics_plot.png b/docs/benchmarks/rotating_cylinders/openfoam_plots/solution_metrics_plot.png new file mode 100644 index 0000000..cf1a3e7 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/openfoam_plots/solution_metrics_plot.png differ diff --git a/docs/benchmarks/rotating_cylinders/openfoam_plots/velocity_field.png b/docs/benchmarks/rotating_cylinders/openfoam_plots/velocity_field.png new file mode 100644 index 0000000..5648fe3 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/openfoam_plots/velocity_field.png differ diff --git a/docs/benchmarks/rotating_cylinders/rotating-cylinders-2d.md b/docs/benchmarks/rotating_cylinders/rotating-cylinders-2d.md new file mode 100644 index 0000000..cd5c728 --- /dev/null +++ b/docs/benchmarks/rotating_cylinders/rotating-cylinders-2d.md @@ -0,0 +1,200 @@ +# Rotating cylinders + +## Problem description + +We consider the case of two concentric cylinders, with the inner cylinder of radius $R_1$ rotating at angular velocity $\Omega_1$ and the outer cylinder of radius $R_2$ rotating at angular velocity $\Omega_2$. The fluid between the cylinders is incompressible and the flow is steady and laminar. + +![Case description](rotating-cylinders-boundary-conditions.png) + +The laminar case corresponds to a Reynolds number of $Re = 100$, where the Reynolds number is defined as: + +$$ +Re = \frac{|u_0| \, d}{\nu} +$$ + +where $u_0 = \Omega_1 R_1$ is the surface velocity of the inner cylinder, $d = R_2 - R_1$ is the gap width between the cylinders, and $\nu$ is the kinematic viscosity of the fluid. Using inner and outer radii of $R_1 = 1$ and $R_2 = 2$, respectively, and kinematic viscosity $\nu = 1$, the angular velocity of the inner cylinder is $\Omega_1 = 100$ rad/s. + +## Analytical solution + +Taylor [@taylor_stability_1922] shows that the rotational velocity $u_{\theta}$ in the annular gap is described by: + +$$ +u_{\theta}(r) = A r + \frac{B}{r} +$$ + +where the constants $A$ and $B$ are defined as: + +$$ +A = \frac{\Omega_2 R_2^2 - \Omega_1 R_1^2}{R_2^2 - R_1^2} +$$ + +$$ +B = \frac{(\Omega_1 - \Omega_2) R_1^2 R_2^2}{R_2^2 - R_1^2} +$$ + +Here, $\Omega_1$ and $\Omega_2$ are the rotational speeds of the inner and outer cylinders, and the ratio $\mu = \Omega_2 / \Omega_1$ can be used to express these alternatively as: + +$$ +A = \frac{\Omega_1 \left(1 - \mu \dfrac{R_2^2}{R_1^2}\right)}{1 - \dfrac{R_2^2}{R_1^2}}, \qquad +B = \frac{R_1^2 \,\Omega_1 (1 - \mu)}{1 - \dfrac{R_1^2}{R_2^2}} +$$ + +The steady flow equations in cylindrical coordinates reduce to a radial pressure balance: + +$$ +\frac{1}{\rho}\frac{\partial p}{\partial r} = \frac{u_{\theta}^2}{r} +$$ + +Integrating with respect to $r$ yields the analytical pressure field: + +$$ +p(r) = \frac{A^2 r^2}{2} + 2AB \ln(r) - \frac{B^2}{2r^2} + C +$$ + +where $C$ is an integration constant determined by a reference pressure condition. + +## Mesh + +The domain is the annular region between the two cylinders, discretized with a 2D structured mesh created using blockMesh. + +![Mesh](rotating-cylinders-mesh.png) + +## Boundary conditions + +The boundary conditions are: + +$$ +\begin{aligned} +u_{\theta} &= \Omega_1 R_1 & \text{on } r = R_1 & \quad \text{(inner cylinder, rotating)} \\ +u_{\theta} &= \Omega_2 R_2 & \text{on } r = R_2 & \quad \text{(outer cylinder, fixed or rotating)} +\end{aligned} +$$ + +In the reference configuration, the outer cylinder is fixed ($\Omega_2 = 0$) and the inner cylinder rotates at $\Omega_1 = 100$ rad/s. + +## Comparison of numerical solution with analytical solution + +The numerical solution $u_{\theta,h}$ and analytical solution $u_{\theta}$ can be compared using the $L_2$ error norm: + +$$ +e_{L_2} = \left\| u_{\theta} - u_{\theta,h} \right\|_{L_2} = \sqrt{\int_\Omega \left(u_{\theta}(\boldsymbol{x}) - u_{\theta,h}(\boldsymbol{x})\right)^2 \mathrm{d}\boldsymbol{x}} +$$ + +and the relative $L_2$ error: + +$$ +e_{L_2,\mathrm{rel}} = \frac{\left\| u_{\theta} - u_{\theta,h} \right\|_{L_2}}{\left\| u_{\theta} \right\|_{L_2}} +$$ + +Analogously, the pressure field $p_h$ can be compared to the analytical pressure $p$. Plotting the relative error against mesh size $h$ in a log–log plot allows determination of the convergence order of the numerical scheme. + +## Table of parameters + +### Model parameters + +| Parameter | Description | +|-----------|-------------| +| $R_1$ [m] | Radius of the inner cylinder. | +| $R_2$ [m] | Radius of the outer cylinder. | +| $\Omega_1$ [rad/s] | Angular velocity of the inner cylinder. | +| $\Omega_2$ [rad/s] | Angular velocity of the outer cylinder. | +| $\nu$ [m²/s] | Kinematic viscosity of the fluid. | +| $Re$ [-] | Reynolds number. | + +### Numerical parameters + +| Parameter | Description | +|-----------|-------------| +| $h$ [m] | Element size. | +| $q$ [-] | Element order (geometry interpolation order). | +| $p$ [-] | Degree of the ansatz functions. | +| $r$ [-] | Degree of the quadrature rule. | +| $\mathcal{Q}$ [-] | Quadrature rule (e.g. Gauss or Gauss–Lobatto). | + +## Numerical results + +### DuMux + +Results obtained with DuMux show good agreement with the analytical Taylor–Couette solution, with errors decreasing consistently across all mesh refinement levels. Three configurations were tested, with the mesh refined by doubling the number of radial and angular cells at each level. + +| Configuration | Radial cells $N_r$ | Angular cells $N_\theta$ | $e_{L^2,\mathrm{rel}}$ velocity | $e_{L^2,\mathrm{rel}}$ pressure | +|---|---|---|---|---| +| `100_0_20_80` | 20 | 80 | 1.33 × 10⁻¹ | 3.79 × 10⁻² | +| `100_0_40_160` | 40 | 160 | 2.76 × 10⁻² | 5.00 × 10⁻³ | +| `100_0_80_320` | 80 | 320 | 9.92 × 10⁻³ | 2.24 × 10⁻³ | + +#### Convergence + +The relative $L^2$ errors for velocity and pressure are plotted against mesh configuration on a log–log scale. Errors decrease consistently with mesh refinement, confirming convergence of the numerical scheme. + +![Convergence of velocity and pressure errors — DuMux](dumux_plots/solution_metrics_plot.png) + +#### Radial profiles vs analytical + +The rotational velocity $u_\theta(r)$ and pressure $p(r)$ profiles for the finest mesh (`100_0_80_320`) are compared against the analytical solution. The numerical results (open red circles) closely follow the analytical curves (black line) across the full gap $r \in [1, 2]$ m. + +![Radial profiles — DuMux vs analytical, mesh 80×320](dumux_plots/analytical_comparison.png) + +#### 2D field plots + +The 2D pressure and velocity magnitude fields for the finest configuration confirm the expected annular structure. The pressure increases centrifugally from the inner to the outer cylinder, and the velocity magnitude decays monotonically from $|u| \approx 100$ m/s at $r = R_1$ to near zero at $r = R_2$. + +
+
+ 2D velocity field — DuMux, mesh 80×320 +
Velocity magnitude u (m/s) — DuMux, mesh 80×320
+
+
+ 2D pressure field — DuMux, mesh 80×320 +
Pressure p (Pa) — DuMux, mesh 80×320
+
+
+ +### OpenFOAM + +Results obtained with OpenFOAM (simpleFoam, laminar, MRF) show good agreement with the analytical Taylor–Couette solution, with errors decreasing consistently across all mesh refinement levels. Three configurations were tested, with the mesh refined by doubling the number of radial and angular cells at each level. + +| Configuration | Radial cells $N_r$ | Angular cells $N_\theta$ | $e_{L^2,\mathrm{rel}}$ velocity | $e_{L^2,\mathrm{rel}}$ pressure | +|---|---|---|---|---| +| `100_0_20_80` | 20 | 80 | 4.69 × 10⁻³ | 1.74 × 10⁻² | +| `100_0_40_160` | 40 | 160 | 8.54 × 10⁻⁴ | 7.33 × 10⁻³ | +| `100_0_80_320` | 80 | 320 | 1.94 × 10⁻⁴ | 3.46 × 10⁻³ | + +#### Convergence + +The relative $L^2$ errors for velocity and pressure are plotted against mesh configuration on a log–log scale. Errors decrease consistently with mesh refinement, confirming convergence of the numerical scheme. + +![Convergence of velocity and pressure errors — OpenFOAM](openfoam_plots/solution_metrics_plot.png) + +#### Radial profiles vs analytical + +The rotational velocity $u_\theta(r)$ and pressure $p(r)$ profiles for the finest mesh (`100_0_80_320`) are compared against the analytical solution. The numerical results closely follow the analytical curves across the full gap $r \in [1, 2]$ m. + +![Radial profiles — OpenFOAM vs analytical, mesh 80×320](openfoam_plots/analytical_comparison.png) + +#### 2D field plots + +The 2D pressure and velocity magnitude fields for the finest configuration confirm the expected annular structure. The pressure increases centrifugally from the inner to the outer cylinder, and the velocity magnitude decays monotonically from $|u| \approx 100$ m/s at $r = R_1$ to near zero at $r = R_2$. + +
+
+ 2D velocity field — OpenFOAM, mesh 80×320 +
Velocity magnitude u (m/s) — OpenFOAM, mesh 80×320
+
+
+ 2D pressure field — OpenFOAM, mesh 80×320 +
Pressure p (Pa) — OpenFOAM, mesh 80×320
+
+
+ +### Solver convergence comparison + +The relative $L^2$ errors for both velocity and pressure are compared across solvers (DuMux and OpenFOAM) on a log–log scale. OpenFOAM achieves lower absolute errors at the coarsest mesh level, while DuMux converges more steeply and reaches comparable or lower errors at finer resolutions, particularly for pressure. + +![Solver convergence comparison — DuMux vs OpenFOAM](convergence_comparison.png) + +## References + +- Taylor, G. I. (1923). Stability of a viscous liquid contained between two rotating cylinders. *Philosophical Transactions of the Royal Society A*, 223, 289–343. +- OpenFOAM Documentation. Rotating Cylinders 2D — Verification & Validation. +- DuMux Benchmark Case. Rotating Cylinders. diff --git a/docs/benchmarks/rotating_cylinders/rotating-cylinders-boundary-conditions.png b/docs/benchmarks/rotating_cylinders/rotating-cylinders-boundary-conditions.png new file mode 100644 index 0000000..85ab6a7 Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/rotating-cylinders-boundary-conditions.png differ diff --git a/docs/benchmarks/rotating_cylinders/rotating-cylinders-mesh.png b/docs/benchmarks/rotating_cylinders/rotating-cylinders-mesh.png new file mode 100644 index 0000000..dc6fbcc Binary files /dev/null and b/docs/benchmarks/rotating_cylinders/rotating-cylinders-mesh.png differ