From 70632b6e0af083d9a3e3eaa7cbf81b46276159e7 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Mon, 23 Feb 2026 16:52:31 -0500 Subject: [PATCH 01/10] add working coalescence case --- examples/3D_droplet_coalescence/case.py | 188 ++++++++++++++++++++++++ 1 file changed, 188 insertions(+) create mode 100644 examples/3D_droplet_coalescence/case.py diff --git a/examples/3D_droplet_coalescence/case.py b/examples/3D_droplet_coalescence/case.py new file mode 100644 index 0000000000..356cf01368 --- /dev/null +++ b/examples/3D_droplet_coalescence/case.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python3 +# Head-on binary droplet coalescence — We = 2, Re ~ 236. +# +# Two water droplets (D = 306 um) in air collide along the x-axis. +# Ported from an OpenFOAM interFoam simulation (temp/ directory). +# +# Nondimensional parameters follow Meng (2016) convention: +# rho_l = 100, rho_g = 1, sigma = 0.72 +# gamma_l = 4.4, pi_inf_l = 100; gamma_g = 1.4, pi_inf_g = 0 +# D = 1, R = 0.5, p_gas = 1.0, p_liq = 1 + sigma/R = 2.44 +# Ur = sqrt(We * sigma / (rho_l * D)) = 0.12 +# Re = rho_l * Ur * D / mu_l = 236 +# +# 3 patches: gas background + two spherical droplets (no hardcoded IC needed). + +import math +import json + +# --- Nondimensional parameters (Meng 2016) --- +rho_l = 100.0 +rho_g = 1.0 +sigma = 0.72 +gamma_l, pi_inf_l = 4.4, 100.0 +gamma_g, pi_inf_g = 1.4, 0.0 +p_gas = 1.0 + +D = 1.0 +R = D / 2.0 +p_liq = p_gas + sigma / R # Laplace pressure jump + +# Collision parameters +We = 2.0 +Re = 236.0 +Ur = math.sqrt(We * sigma / (rho_l * D)) # ~0.12 relative velocity +Ud = Ur / 2.0 # each droplet moves at half the relative velocity + +# Viscosity (matching Re = 236) +mu_l = rho_l * Ur * D / Re # ~0.05085 +mu_g = mu_l / 48.1 # ~0.001057 (physical viscosity ratio water/air) + +# Domain: 9D x 6D x 6D centered at origin +x0, x1 = -4.5, 4.5 +y0, y1 = -3.0, 3.0 +z0, z1 = -3.0, 3.0 + +# Grid: 10 cells/D +Nx = 179 # 90 cells in x (9D * 10) +Ny = 119 # 60 cells in y (6D * 10) +Nz = 119 # 60 cells in z (6D * 10) + +eps = 1e-9 + +# Droplet placement: centers at (+-sep, 0, 0) +# Minimal gap (~0.3D) to avoid wake-induced deformation before contact +sep = 0.65 * D # center-to-center = 1.3D, surface gap = 0.3D + +# Time stepping +dx = (x1 - x0) / (Nx + 1) +c_l = math.sqrt(gamma_l * (p_gas + pi_inf_l) / rho_l) # ~2.11 +mydt = 0.1 * dx / c_l # CFL = 0.1 + +# Collision timescale = D / Ur ~ 8.33 +# Run for ~11 collision times (matching 5ms OpenFOAM end time) +t_end = 11.0 * D / Ur # ~92 +Nt = int(math.ceil(t_end / mydt)) +Ns = max(1, Nt // 100) # ~100 output frames + +# Configuration case dictionary +data = { + # Logistics + "run_time_info": "T", + # Computational Domain + "x_domain%beg": x0, + "x_domain%end": x1, + "y_domain%beg": y0, + "y_domain%end": y1, + "z_domain%beg": z0, + "z_domain%end": z1, + "m": Nx, + "n": Ny, + "p": Nz, + "cyl_coord": "F", + "dt": mydt, + "t_step_start": 0, + "t_step_stop": Nt, + "t_step_save": Ns, + # Simulation Algorithm + "model_eqns": 3, + "alt_soundspeed": "F", + "mixture_err": "T", + "mpp_lim": "F", + "time_stepper": 3, + "weno_order": 3, + "avg_state": 2, + "weno_eps": 1e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "weno_Re_flux": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "viscous": "T", + "bc_x%beg": -6, + "bc_x%end": -6, + "bc_y%beg": -6, + "bc_y%end": -6, + "bc_z%beg": -6, + "bc_z%end": -6, + "num_patches": 3, + "num_fluids": 2, + "weno_avg": "T", + "surface_tension": "T", + # Database Structure Parameters + "format": 1, + "precision": 2, + "alpha_wrt(1)": "T", + "cf_wrt": "T", + "vel_wrt(1)": "T", + "vel_wrt(2)": "T", + "vel_wrt(3)": "T", + "parallel_io": "T", + "sigma": sigma, + # Fluid Parameters (Liquid — fluid 1) + "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), + "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), + "fluid_pp(1)%Re(1)": 1.0 / mu_l, + # Fluid Parameters (Gas — fluid 2) + "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), + "fluid_pp(2)%pi_inf": 0.0, + "fluid_pp(2)%Re(1)": 1.0 / mu_g, + # ===== Patch 1: Gas background (fills entire domain) ===== + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%z_centroid": 0.0, + "patch_icpp(1)%length_x": x1 - x0, + "patch_icpp(1)%length_y": y1 - y0, + "patch_icpp(1)%length_z": z1 - z0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": p_gas, + "patch_icpp(1)%alpha_rho(1)": eps * rho_l, + "patch_icpp(1)%alpha_rho(2)": (1.0 - eps) * rho_g, + "patch_icpp(1)%alpha(1)": eps, + "patch_icpp(1)%alpha(2)": 1.0 - eps, + "patch_icpp(1)%cf_val": 0, + # ===== Patch 2: Left droplet (center at -sep, moves in +x) ===== + "patch_icpp(2)%geometry": 8, + "patch_icpp(2)%x_centroid": -sep, + "patch_icpp(2)%y_centroid": 0.0, + "patch_icpp(2)%z_centroid": 0.0, + "patch_icpp(2)%radius": R, + "patch_icpp(2)%alter_patch(1)": "T", + "patch_icpp(2)%smoothen": "T", + "patch_icpp(2)%smooth_patch_id": 1, + "patch_icpp(2)%smooth_coeff": 0.95, + "patch_icpp(2)%vel(1)": Ud, + "patch_icpp(2)%vel(2)": 0.0, + "patch_icpp(2)%vel(3)": 0.0, + "patch_icpp(2)%pres": p_liq, + "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_l, + "patch_icpp(2)%alpha_rho(2)": eps * rho_g, + "patch_icpp(2)%alpha(1)": 1.0 - eps, + "patch_icpp(2)%alpha(2)": eps, + "patch_icpp(2)%cf_val": 1, + # ===== Patch 3: Right droplet (center at +sep, moves in -x) ===== + "patch_icpp(3)%geometry": 8, + "patch_icpp(3)%x_centroid": sep, + "patch_icpp(3)%y_centroid": 0.0, + "patch_icpp(3)%z_centroid": 0.0, + "patch_icpp(3)%radius": R, + "patch_icpp(3)%alter_patch(1)": "T", + "patch_icpp(3)%smoothen": "T", + "patch_icpp(3)%smooth_patch_id": 1, + "patch_icpp(3)%smooth_coeff": 0.95, + "patch_icpp(3)%vel(1)": -Ud, + "patch_icpp(3)%vel(2)": 0.0, + "patch_icpp(3)%vel(3)": 0.0, + "patch_icpp(3)%pres": p_liq, + "patch_icpp(3)%alpha_rho(1)": (1.0 - eps) * rho_l, + "patch_icpp(3)%alpha_rho(2)": eps * rho_g, + "patch_icpp(3)%alpha(1)": 1.0 - eps, + "patch_icpp(3)%alpha(2)": eps, + "patch_icpp(3)%cf_val": 1, +} + +print(json.dumps(data)) From 56682d83240d2a7c68fabebebbe1bf2927d966ae Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Fri, 27 Feb 2026 11:23:24 -0500 Subject: [PATCH 02/10] Tee interactive run output to .out file Batch runs get a .out log file from the scheduler, but interactive runs had no equivalent. Wrap the interactive command with tee so output streams to both the terminal and {name}.out. Co-Authored-By: Claude Opus 4.6 --- toolchain/mfc/run/run.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/toolchain/mfc/run/run.py b/toolchain/mfc/run/run.py index 95747c3500..a7aa58e839 100644 --- a/toolchain/mfc/run/run.py +++ b/toolchain/mfc/run/run.py @@ -143,6 +143,7 @@ def __execute_job_script(qsystem: queues.QueueSystem): # any files the queue system generates (like .err and .out) are created # in the correct directory. cmd = qsystem.gen_submit_cmd(__job_script_filepath()) + cwd = os.path.dirname(ARG("input")) verbosity = ARG('verbose') @@ -151,11 +152,18 @@ def __execute_job_script(qsystem: queues.QueueSystem): cons.print(f" [dim]$ {' '.join(str(c) for c in cmd)}[/dim]") cons.print() + # For interactive runs on Unix, tee output to a .out file so users get + # the same log file that batch runs produce via the scheduler. + if isinstance(qsystem, queues.InteractiveSystem) and os.name != 'nt': + out_filepath = os.path.abspath(os.path.join(cwd, f"{ARG('name')}.out")) + cmd = ["/bin/bash", "-c", + f"set -o pipefail; {shlex.join([str(x) for x in cmd])} 2>&1 | tee {shlex.quote(out_filepath)}"] + # Execute the job script with appropriate output handling # At verbosity >= 2, show print_cmd=True for system() calls print_cmd = verbosity >= 2 - if system(cmd, cwd=os.path.dirname(ARG("input")), print_cmd=print_cmd).returncode != 0: + if system(cmd, cwd=cwd, print_cmd=print_cmd).returncode != 0: raise MFCException(f"Submitting batch file for {qsystem.name} failed. It can be found here: {__job_script_filepath()}. Please check the file for errors.") From aeecc0c9dc3902543380e3c94e0bc23bb2349c53 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Tue, 21 Apr 2026 21:11:53 -0400 Subject: [PATCH 03/10] Add droplet coalescence case setups for acoustic demulsification project - New 2D axisymmetric drainage validation case (Phase 1) - New 3D MUSCL+THINC variant (a-g parameter sweep) - Update existing 3D WENO-Z case - Ignore *.md at repo root --- .gitignore | 1 + .../2D_axisym_droplet_coalescence/case.py | 338 ++++++++++++++++ examples/3D_droplet_coalescence/case.py | 364 ++++++++++-------- examples/3D_droplet_coalescence/case_muscl.py | 193 ++++++++++ 4 files changed, 743 insertions(+), 153 deletions(-) create mode 100644 examples/2D_axisym_droplet_coalescence/case.py create mode 100644 examples/3D_droplet_coalescence/case_muscl.py diff --git a/.gitignore b/.gitignore index aba54411e1..f17f581ede 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +*.md node_modules/ package.json yarn.lock diff --git a/examples/2D_axisym_droplet_coalescence/case.py b/examples/2D_axisym_droplet_coalescence/case.py new file mode 100644 index 0000000000..5bb8c5a7f4 --- /dev/null +++ b/examples/2D_axisym_droplet_coalescence/case.py @@ -0,0 +1,338 @@ +#!/usr/bin/env python3 +# Phase 1 drainage validation case — acoustic_demulsification project +# Geometry: 2D axisymmetric, head-on collision of two equal droplets +# Fluids : tetradecane drops in nitrogen ambient +# Target : resolve gas-film drainage; compare h_min(t) to Klaseboer 2000 +# Stop : when h_min reaches ~100–200 nm (before numerical rupture) +# Roadmap reference: +Projects/acoustic_demulsification/simulation_roadmap.md +# (Phase 1 — Drainage validation) +# Run: +# ./mfc.sh run case.py -n # single node +# ./mfc.sh run case.py --gpu -n # GPU build +# Knobs at top of file: +# SMOKE_TEST — True: fast end-to-end smoke test (no stretch, fixed dt, 100 steps) +# False: production (stretched grid, CFL-adaptive dt, t_stop = 200 us) +# SCHEME — 'wenoz' or 'muscl_thinc' +# DX_FINE — cell size in drainage film (m) +# DX_BULK — cell size in bulk (m) +# FILM_HALF — axial half-width of the uniform fine zone (m) +# WE — target Weber number (default 30) +# Known issues: +# - With SMOKE_TEST = False and stretching enabled, simulation hits NaN at +# corner cell (m, n) after ~10 steps. Cause under investigation. The +# smoke-test path (uniform grid) runs end-to-end cleanly. + +import json +import math + +# User-facing knobs + +SMOKE_TEST = True # True: coarse + uniform + 100 fixed-dt steps +SCHEME = "wenoz" # 'wenoz' | 'muscl_thinc' +FILM_HALF = 2e-6 # ± half-width of the uniform fine zone on the axis (m) +WE = 30.0 # Weber number (coalescence regime) + +# Droplet diameter defined here so DX_FINE can reference it; also used below. +D = 300e-6 # droplet diameter, 300 um + +# DX_FINE = D/1000 = 300 nm for both modes. Sets the UNIFORM PRE-STRETCH cell +# size domain/N (see Nx/Ny derivation below). MFC's tanh stretch keeps cells +# inside [x_a, x_b] at ~this size and grows cells outward; effective bulk +# cell size is emergent, not a user knob. +DX_FINE = D / 1000.0 + +# Optional explicit overrides for (Nx, Ny). If set, DX_FINE becomes diagnostic. +NX_OVERRIDE = None +NY_OVERRIDE = None + +# Physical constants (tetradecane / nitrogen at ~25 C) + +# Tetradecane (liquid) +rho_l = 762.0 # kg/m^3 +mu_l = 2.1e-3 # Pa*s +sigma = 0.027 # N/m (≈ 27 mN/m) + +# Stiffened-gas fit for tetradecane: gamma = 4.4, p_inf = 3.0e8 Pa +# => c_l = sqrt(gamma*(p+p_inf)/rho) ≈ 1316 m/s at 1 atm — matches bulk data. +gamma_l = 4.4 +pi_inf_l = 3.0e8 + +# Nitrogen (ideal gas) +rho_g = 1.165 # kg/m^3 at 1 atm, 20 C +gamma_g = 1.4 +pi_inf_g = 0.0 + +p_amb = 1.01325e5 # 1 atm + +# Droplet geometry and kinematics + +R = D / 2.0 # radius, 150 um (D = 300 um defined at top) + +# Weber number: We = rho_l * U_rel^2 * D / sigma, U_rel = 2U per Qian & Law. +# U_rel = 1.8823 m/s -> U = 0.9412 m/s per droplet +U_rel = math.sqrt(WE * sigma / (rho_l * D)) +U = 0.5 * U_rel + +# Initial axial gap between droplet surfaces (pre-drainage), and centers. +gap0 = 20e-6 # 20 um initial gap +x_center = R + 0.5 * gap0 # droplet centers at x = +/- x_center + +# Domain & grid + +Lx_half = 5.0 * R # = 750 um axial half-domain +Ly_max = 3.0 * R # = 450 um radial extent + + +def _post_stretch_extent(lo, hi, a, xa, xb, loops): + """Mimic MFC's hyperbolic-tangent stretching in pre_process/m_grid.f90 to + compute the post-stretch domain extent. The formula uses normalized coords: + x <- x / length; xa <- xa / length; xb <- xb / length; + x <- x/a * (a + log(cosh(a*(x-xa))) + log(cosh(a*(x-xb))) - 2*log(cosh(a*(xb-xa)/2))) + x <- x * length + This expands the domain beyond [lo, hi]; patches must cover the result or + the outer cells are left uninitialized and the sim NaNs at step 1. + """ + length = abs(hi - lo) + x0, x1 = lo / length, hi / length + xa_n, xb_n = xa / length, xb / length + for _ in range(loops): + const = 2.0 * math.log(math.cosh(a * (xb_n - xa_n) / 2.0)) + + def f(x): + return x / a * (a + math.log(math.cosh(a * (x - xa_n))) + math.log(math.cosh(a * (x - xb_n))) - const) + + x0, x1 = f(x0), f(x1) + return x0 * length, x1 * length + + +# Cell count to hit DX_FINE in the fine zone. +# MFC's tanh stretch (m_grid.f90 s_generate_serial_grid) starts with N uniform +# cells of size dx = domain/N, then remaps their boundaries. The remap +# derivative is ~1 inside [x_a, x_b] (cells stay at ~dx) and grows outside +# (cells expand). So to achieve DX_FINE in the fine zone: +# N = ceil(domain / DX_FINE) +# Bulk cell sizes after stretch are set by (a_*, loops_*) and are reported in +# the stderr diagnostic at the bottom of this file. +Nx = NX_OVERRIDE if NX_OVERRIDE is not None else int(math.ceil(2.0 * Lx_half / DX_FINE)) +Ny = NY_OVERRIDE if NY_OVERRIDE is not None else int(math.ceil(Ly_max / DX_FINE)) + +# Numerical scheme switches + +if SCHEME == "wenoz": + recon = { + "recon_type": 1, + "weno_order": 5, + "weno_eps": 1e-6, + "mapped_weno": "F", + "wenoz": "T", + "mp_weno": "T", + "teno": "F", + } +elif SCHEME == "muscl_thinc": + recon = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": 4, # Van Leer + "int_comp": "T", # THINC interface compression + "ic_eps": 1e-4, + "ic_beta": 1.6, + } +else: + raise ValueError(f"Unknown SCHEME: {SCHEME!r}") + +# Time stepping & grid stretching (depends on SMOKE_TEST) + +if SMOKE_TEST: + # Inching roadmap (edit dt / t_step_stop to advance stages): + # stage 0 (old): dt=1e-12, n=100, save=10 -> 100 ps total, 13 s/run + # stage 1 (now): dt=1e-11, n=10000, save=1000 -> 100 ns total, ~21 min @-n1 + # stage 2: dt=1e-11, n=100000, save=10000 -> 1 µs total, ~3.5 hr @-n1 + # stage 3: switch to cfl_adap_dt + t_stop=5e-6 (see production block) + # stage 4: refine DX_FINE to 20 nm -> 10 nm + # stage 5: enable stretch_y (currently NaNs at step ~10 — debug first) + # CFL on 40-nm grid with c_liquid=1316 m/s: dt_max(CFL=0.3)=9.1e-12 s. + # dt=1e-11 sits at CFL~0.33 — safely inside stability. + time_block = { + "dt": 1.0e-11, + "t_step_start": 0, + "t_step_stop": 10000, + "t_step_save": 1000, + } + stretch_block = { + "stretch_x": "T", + "a_x": 2.0, + "x_a": -FILM_HALF, + "x_b": +FILM_HALF, + "loops_x": 1, + "stretch_y": "F", + } +else: + # Production: CFL-adaptive dt, stretched grid. + # t_stop = 2e-4 s (~1.25 * tau_inertia = 1.59e-4 s), 40 frames. + time_block = { + "cfl_adap_dt": "T", + "cfl_target": 0.3, + "t_stop": 2.0e-4, + "t_save": 5.0e-6, + "n_start": 0, + } + # Production: uniform 300 nm cells everywhere (Nx/Ny already large enough + # to hit DX_FINE without tanh concentration). Disabling stretch keeps the + # domain at the intended [-Lx_half, +Lx_half] x [0, Ly_max]; with a=4 + # loops=2 the tanh remap was expanding to ±2892 um x 7 mm and wasting + # cells on empty bulk. Re-enable if you need non-uniform spacing later. + stretch_block = { + "stretch_x": "F", + "stretch_y": "F", + } + +# Compute post-stretch domain extent so patches cover the actual cell range. +# Add 10% margin on top, to account for cell centers vs boundaries. +if stretch_block["stretch_x"] == "T": + _xlo, _xhi = _post_stretch_extent(-Lx_half, +Lx_half, stretch_block["a_x"], stretch_block["x_a"], stretch_block["x_b"], stretch_block["loops_x"]) + patch1_length_x = 1.10 * (_xhi - _xlo) +else: + patch1_length_x = 2.0 * Lx_half + +if stretch_block["stretch_y"] == "T": + _ylo, _yhi = _post_stretch_extent(0.0, Ly_max, stretch_block["a_y"], stretch_block["y_a"], stretch_block["y_b"], stretch_block["loops_y"]) + patch1_length_y = 1.10 * (_yhi - _ylo) +else: + patch1_length_y = Ly_max + +# MFC input deck + +eps = 1e-9 # avoid exactly 0/1 volume fractions + +case = { + # grid + "x_domain%beg": -Lx_half, + "x_domain%end": +Lx_half, + "y_domain%beg": 0.0, + "y_domain%end": +Ly_max, + "m": Nx - 1, # MFC uses zero-based cell-count (m = Nx-1) + "n": Ny - 1, + "p": 0, # 2D + "cyl_coord": "T", # 2D axisymmetric: x is axial, y is radial + # boundary conditions + # Axial (x): outflow/extrapolation on both ends. + "bc_x%beg": -6, # TRIAGE: char non-reflecting (was -3) + "bc_x%end": -6, + # Radial (y): axisymmetric at y=0, outflow at y_max. + "bc_y%beg": -2, # axis + "bc_y%end": -6, # TRIAGE: char non-reflecting (was -3) + # time integration + "time_stepper": 3, # TVD-RK3 + # physical model + "model_eqns": 2, # 5-equation (Allaire) diffuse interface + "num_fluids": 2, + "num_patches": 3, + "mixture_err": "T", + "mpp_lim": "T", + "avg_state": 2, + "riemann_solver": 2, # HLLC + "wave_speeds": 1, + # surface tension + "surface_tension": "T", + "sigma": sigma, + # fluids (index 1 = tetradecane, 2 = nitrogen) + "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), # 1/(4.4-1) = 0.2941 + "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), # 4.4*3e8/3.4 ≈ 3.882e8 + "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), # 1/(1.4-1) = 2.5 + "fluid_pp(2)%pi_inf": 0.0, + # patch 1: nitrogen background filling the whole domain + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": patch1_length_y / 2.0, + "patch_icpp(1)%length_x": patch1_length_x, + "patch_icpp(1)%length_y": patch1_length_y, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": p_amb, + "patch_icpp(1)%alpha_rho(1)": eps * rho_l, + "patch_icpp(1)%alpha_rho(2)": (1.0 - eps) * rho_g, + "patch_icpp(1)%alpha(1)": eps, + "patch_icpp(1)%alpha(2)": 1.0 - eps, + "patch_icpp(1)%cf_val": 0.0, + # patch 2: left droplet, moving in +x + "patch_icpp(2)%alter_patch(1)": "T", + "patch_icpp(2)%smoothen": "T", + "patch_icpp(2)%smooth_patch_id": 1, + "patch_icpp(2)%smooth_coeff": 0.95, + "patch_icpp(2)%geometry": 2, # circle (axisymmetric disk becomes a sphere) + "patch_icpp(2)%x_centroid": -x_center, + "patch_icpp(2)%y_centroid": 0.0, + "patch_icpp(2)%radius": R, + "patch_icpp(2)%vel(1)": +U, + "patch_icpp(2)%vel(2)": 0.0, + "patch_icpp(2)%pres": p_amb, + "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_l, + "patch_icpp(2)%alpha_rho(2)": eps * rho_g, + "patch_icpp(2)%alpha(1)": 1.0 - eps, + "patch_icpp(2)%alpha(2)": eps, + "patch_icpp(2)%cf_val": 1.0, + # patch 3: right droplet, moving in -x + "patch_icpp(3)%alter_patch(1)": "T", + "patch_icpp(3)%smoothen": "T", + "patch_icpp(3)%smooth_patch_id": 1, + "patch_icpp(3)%smooth_coeff": 0.95, + "patch_icpp(3)%geometry": 2, + "patch_icpp(3)%x_centroid": +x_center, + "patch_icpp(3)%y_centroid": 0.0, + "patch_icpp(3)%radius": R, + "patch_icpp(3)%vel(1)": -U, + "patch_icpp(3)%vel(2)": 0.0, + "patch_icpp(3)%pres": p_amb, + "patch_icpp(3)%alpha_rho(1)": (1.0 - eps) * rho_l, + "patch_icpp(3)%alpha_rho(2)": eps * rho_g, + "patch_icpp(3)%alpha(1)": 1.0 - eps, + "patch_icpp(3)%alpha(2)": eps, + "patch_icpp(3)%cf_val": 1.0, + # post processing + "alpha_wrt(1)": "T", # liquid volume fraction (the field h_min lives in) + "alpha_wrt(2)": "T", # gas volume fraction + "rho_wrt": "T", + "pres_wrt": "T", + "mom_wrt(1)": "T", + "mom_wrt(2)": "T", + "cf_wrt": "T", # color function (requires surface_tension) + "parallel_io": "T", + "format": 1, # 1 = silo-HDF5; 2 = binary +} + +case.update(recon) +case.update(time_block) +case.update(stretch_block) + +# Diagnostic: report ACHIEVED cell sizes after MFC's tanh stretch (stderr so +# it doesn't pollute the JSON on stdout). +import sys as _sys + + +def _stretched_cb(m, beg, end, stretch, a, xa, xb, loops): + cb = [beg + (end - beg) * i / (m + 1) for i in range(-1, m + 1)] # m+2 bnds + if stretch: + length = abs(cb[-1] - cb[0]) + v = [x / length for x in cb] + xa_n, xb_n = xa / length, xb / length + for _ in range(loops): + c = 2.0 * math.log(math.cosh(a * (xb_n - xa_n) / 2.0)) + v = [x / a * (a + math.log(math.cosh(a * (x - xa_n))) + math.log(math.cosh(a * (x - xb_n))) - c) for x in v] + cb = [x * length for x in v] + return cb + + +def _report(axis, m, beg, end): + sk_stretch = case.get(f"stretch_{axis}", "F") == "T" + cb = _stretched_cb(m, beg, end, sk_stretch, case.get(f"a_{axis}", 1.0), case.get(f"{axis}_a", 0.0), case.get(f"{axis}_b", 0.0), case.get(f"loops_{axis}", 0)) + cells = [cb[i + 1] - cb[i] for i in range(m + 1)] + _sys.stderr.write(f"[case.py] {axis}: N={m + 1:6d} min={min(cells) * 1e9:>9.1f} nm max={max(cells) * 1e6:>8.2f} um span=[{cb[0] * 1e6:8.1f},{cb[-1] * 1e6:8.1f}] um\n") + + +_sys.stderr.write(f"[case.py] SMOKE_TEST={SMOKE_TEST} DX_FINE target={DX_FINE * 1e9:.0f} nm\n") +_report("x", Nx - 1, -Lx_half, +Lx_half) +_report("y", Ny - 1, 0.0, +Ly_max) +_sys.stderr.write(f"[case.py] total cells = {Nx * Ny:,}\n") +_sys.stderr.flush() + +print(json.dumps(case)) diff --git a/examples/3D_droplet_coalescence/case.py b/examples/3D_droplet_coalescence/case.py index 356cf01368..d01308b95f 100644 --- a/examples/3D_droplet_coalescence/case.py +++ b/examples/3D_droplet_coalescence/case.py @@ -1,188 +1,246 @@ #!/usr/bin/env python3 -# Head-on binary droplet coalescence — We = 2, Re ~ 236. -# -# Two water droplets (D = 306 um) in air collide along the x-axis. -# Ported from an OpenFOAM interFoam simulation (temp/ directory). -# -# Nondimensional parameters follow Meng (2016) convention: -# rho_l = 100, rho_g = 1, sigma = 0.72 -# gamma_l = 4.4, pi_inf_l = 100; gamma_g = 1.4, pi_inf_g = 0 -# D = 1, R = 0.5, p_gas = 1.0, p_liq = 1 + sigma/R = 2.44 -# Ur = sqrt(We * sigma / (rho_l * D)) = 0.12 -# Re = rho_l * Ur * D / mu_l = 236 -# -# 3 patches: gas background + two spherical droplets (no hardcoded IC needed). +# Phase 1 drainage validation case — acoustic_demulsification project +# Geometry: 2D axisymmetric, head-on collision of two equal droplets +# Fluids : tetradecane drops in nitrogen ambient +# Target : resolve gas-film drainage; compare h_min(t) to Klaseboer 2000 +# Stop : when h_min reaches ~100–200 nm (before numerical rupture) +# Roadmap reference: +Projects/acoustic_demulsification/simulation_roadmap.md +# (Phase 1 — Drainage validation) +# Run: +# ./mfc.sh run case.py -n # single node +# ./mfc.sh run case.py --gpu -n # GPU build +# Knobs at top of file: +# SCHEME — 'wenoz' or 'muscl_thinc' +# DX_FINE — cell size in drainage film (m); start 10 nm +# DX_BULK — cell size in bulk (m); start 1 micron +# FILM_HALF — axial half-width of the uniform fine zone (m) +# WE — target Weber number (default 30) -import math import json +import math + +# User-facing knobs + +SCHEME = "wenoz" # 'wenoz' | 'muscl_thinc' +DX_FINE = 10e-9 # 10 nm — film cell size target +DX_BULK = 1e-6 # 1 um — bulk cell size target +FILM_HALF = 2e-6 # ± half-width of the uniform fine zone on the axis (m) +WE = 30.0 # Weber number (coalescence regime) + +# For quick debugging, set coarser values, e.g. DX_FINE=40e-9, DX_BULK=4e-6. + +# Physical constants (tetradecane / nitrogen at ~25 C) + +# Tetradecane (liquid) +rho_l = 762.0 # kg/m^3 +mu_l = 2.1e-3 # Pa*s +sigma = 0.027 # N/m (≈ 27 mN/m) + +# Stiffened-gas fit for tetradecane: gamma = 4.4, p_inf = 3.0e8 Pa +# => c_l = sqrt(gamma*(p+p_inf)/rho) ≈ 1316 m/s at 1 atm — matches bulk data. +gamma_l = 4.4 +pi_inf_l = 3.0e8 + +# Nitrogen (ideal gas) +rho_g = 1.165 # kg/m^3 at 1 atm, 20 C +gamma_g = 1.4 +pi_inf_g = 0.0 + +p_amb = 1.01325e5 # 1 atm + +# Droplet geometry and kinematics + +D = 300e-6 # droplet diameter, 300 um +R = D / 2.0 # radius, 150 um + +# Weber number: We = rho_l * U_rel^2 * D / sigma, U_rel = 2U per Qian & Law. +# U_rel^2 = We*sigma/(rho_l*D) +# = 30 * 0.027 / (762 * 3e-4) +# = 0.810 / 0.2286 +# = 3.5433 m^2/s^2 +# U_rel = 1.8823 m/s -> U = 0.9412 m/s per droplet +U_rel = math.sqrt(WE * sigma / (rho_l * D)) +U = 0.5 * U_rel + +# Initial axial gap between droplet surfaces (pre-drainage), and centers. +gap0 = 20e-6 # 20 um initial gap — leaves inertial approach + drainage +x_center = R + 0.5 * gap0 # droplet centers at x = +/- x_center + +# Domain & grid + +# Axial extent — big enough that outgoing waves don't reflect back onto the film. +# Keep a few D each side of the droplets. +Lx_half = 5.0 * R # = 750 um axial half-domain +Ly_max = 3.0 * R # = 450 um radial extent + +# Uniform fine zone in x: [-FILM_HALF, +FILM_HALF] at DX_FINE. +# Stretched zone: DX_FINE -> DX_BULK, geometric growth handled by MFC stretch_x. +# Cell-count estimate (for sanity; MFC's hyperbolic stretch won't be exactly geometric): +# N_fine = 2*FILM_HALF / DX_FINE +# N_stretch (each side) ~ ln(DX_BULK/DX_FINE)/ln(r) with r = 1.1 -> ~48 +N_fine = int(2 * FILM_HALF / DX_FINE) # 400 for defaults +N_stretch = int(math.ceil(math.log(DX_BULK / DX_FINE) / math.log(1.1))) # ~48 each side +Nx = N_fine + 2 * N_stretch # 496 for defaults + +# Radial: fine near axis (film rim), coarsens outward. +# Same strategy on y (positive only, because axisymmetric). +N_fine_r = int(4 * FILM_HALF / DX_FINE) # 800 cells out to 8 um +N_stretch_r = int(math.ceil(math.log(DX_BULK / DX_FINE) / math.log(1.1))) # ~48 +Ny = N_fine_r + N_stretch_r # 848 for defaults + +# NOTE: These are generous baselines for a grid-convergence study. +# On first debug runs, set DX_FINE=40e-9, DX_BULK=4e-6 to cut Nx,Ny by ~4x. -# --- Nondimensional parameters (Meng 2016) --- -rho_l = 100.0 -rho_g = 1.0 -sigma = 0.72 -gamma_l, pi_inf_l = 4.4, 100.0 -gamma_g, pi_inf_g = 1.4, 0.0 -p_gas = 1.0 - -D = 1.0 -R = D / 2.0 -p_liq = p_gas + sigma / R # Laplace pressure jump - -# Collision parameters -We = 2.0 -Re = 236.0 -Ur = math.sqrt(We * sigma / (rho_l * D)) # ~0.12 relative velocity -Ud = Ur / 2.0 # each droplet moves at half the relative velocity - -# Viscosity (matching Re = 236) -mu_l = rho_l * Ur * D / Re # ~0.05085 -mu_g = mu_l / 48.1 # ~0.001057 (physical viscosity ratio water/air) - -# Domain: 9D x 6D x 6D centered at origin -x0, x1 = -4.5, 4.5 -y0, y1 = -3.0, 3.0 -z0, z1 = -3.0, 3.0 - -# Grid: 10 cells/D -Nx = 179 # 90 cells in x (9D * 10) -Ny = 119 # 60 cells in y (6D * 10) -Nz = 119 # 60 cells in z (6D * 10) - -eps = 1e-9 - -# Droplet placement: centers at (+-sep, 0, 0) -# Minimal gap (~0.3D) to avoid wake-induced deformation before contact -sep = 0.65 * D # center-to-center = 1.3D, surface gap = 0.3D - -# Time stepping -dx = (x1 - x0) / (Nx + 1) -c_l = math.sqrt(gamma_l * (p_gas + pi_inf_l) / rho_l) # ~2.11 -mydt = 0.1 * dx / c_l # CFL = 0.1 - -# Collision timescale = D / Ur ~ 8.33 -# Run for ~11 collision times (matching 5ms OpenFOAM end time) -t_end = 11.0 * D / Ur # ~92 -Nt = int(math.ceil(t_end / mydt)) -Ns = max(1, Nt // 100) # ~100 output frames - -# Configuration case dictionary -data = { - # Logistics - "run_time_info": "T", - # Computational Domain - "x_domain%beg": x0, - "x_domain%end": x1, - "y_domain%beg": y0, - "y_domain%end": y1, - "z_domain%beg": z0, - "z_domain%end": z1, - "m": Nx, - "n": Ny, - "p": Nz, - "cyl_coord": "F", - "dt": mydt, - "t_step_start": 0, - "t_step_stop": Nt, - "t_step_save": Ns, - # Simulation Algorithm - "model_eqns": 3, - "alt_soundspeed": "F", +# Time integration + +# Characteristic times: +# tau_inertia = D / U_rel = 3e-4 / 1.8823 ≈ 1.59e-4 s +# tau_cap = sqrt(rho_l D^3 / sigma) ≈ 8.73e-4 s +# Want to simulate until h_min ~ 100–200 nm. Experience says ~1–2 * tau_inertia. +t_stop = 2.0e-4 # 200 us +t_save = 5.0e-6 # 40 frames + +# Numerical scheme switches + +# WENO-Z 5 is the default from MFC 5.0 paper; MUSCL+THINC is the fallback if +# interface width exceeds ~5 cells or if we see early numerical rupture. +if SCHEME == "wenoz": + recon = { + "recon_type": 1, + "weno_order": 5, + "weno_eps": 1e-6, + "mapped_weno": "F", + "wenoz": "T", + "mp_weno": "T", + "teno": "F", + } +elif SCHEME == "muscl_thinc": + recon = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": 4, # Van Leer + "int_comp": "T", # THINC interface compression + "ic_eps": 1e-4, + "ic_beta": 1.6, + } +else: + raise ValueError(f"Unknown SCHEME: {SCHEME!r}") + +# MFC input deck + +eps = 1e-9 # avoid exactly 0/1 volume fractions + +case = { + # grid + "x_domain%beg": -Lx_half, + "x_domain%end": +Lx_half, + "y_domain%beg": 0.0, + "y_domain%end": +Ly_max, + "m": Nx - 1, # MFC uses zero-based cell-count (m = Nx-1) + "n": Ny - 1, + "p": 0, # 2D + "cyl_coord": "T", # 2D axisymmetric: x is axial, y is radial + # Grid stretching — fine uniform zone in [-FILM_HALF, +FILM_HALF], stretched outside. + "stretch_x": "T", + "a_x": 4.0, + "x_a": -FILM_HALF, + "x_b": +FILM_HALF, + "loops_x": 2, + "stretch_y": "T", + "a_y": 4.0, + "y_a": 0.0, + "y_b": 4.0 * FILM_HALF, # keep fine near the axis / film rim + "loops_y": 2, + # boundary conditions + # Axial (x): outflow/extrapolation on both ends. + "bc_x%beg": -3, + "bc_x%end": -3, + # Radial (y): axisymmetric at y=0, outflow at y_max. + "bc_y%beg": -2, # reflective / axis + "bc_y%end": -3, + # time integration + "time_stepper": 3, # TVD-RK3 + "cfl_adap_dt": "T", + "cfl_target": 0.3, # conservative; tighten/loosen after first runs + "t_stop": t_stop, + "t_save": t_save, + # physical model + "model_eqns": 2, # 5-equation (Allaire) diffuse interface + "num_fluids": 2, + "num_patches": 3, "mixture_err": "T", - "mpp_lim": "F", - "time_stepper": 3, - "weno_order": 3, + "mpp_lim": "T", "avg_state": 2, - "weno_eps": 1e-16, - "mapped_weno": "T", - "null_weights": "F", - "mp_weno": "F", - "weno_Re_flux": "T", - "riemann_solver": 2, + "riemann_solver": 2, # HLLC "wave_speeds": 1, - "viscous": "T", - "bc_x%beg": -6, - "bc_x%end": -6, - "bc_y%beg": -6, - "bc_y%end": -6, - "bc_z%beg": -6, - "bc_z%end": -6, - "num_patches": 3, - "num_fluids": 2, - "weno_avg": "T", + **recon, + # surface tension "surface_tension": "T", - # Database Structure Parameters - "format": 1, - "precision": 2, - "alpha_wrt(1)": "T", - "cf_wrt": "T", - "vel_wrt(1)": "T", - "vel_wrt(2)": "T", - "vel_wrt(3)": "T", - "parallel_io": "T", "sigma": sigma, - # Fluid Parameters (Liquid — fluid 1) - "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), - "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), - "fluid_pp(1)%Re(1)": 1.0 / mu_l, - # Fluid Parameters (Gas — fluid 2) - "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), + # fluids (index 1 = tetradecane, 2 = nitrogen) + "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), # 1/(4.4-1) = 0.2941 + "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), # 4.4*3e8/3.4 ≈ 3.882e8 + "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), # 1/(1.4-1) = 2.5 "fluid_pp(2)%pi_inf": 0.0, - "fluid_pp(2)%Re(1)": 1.0 / mu_g, - # ===== Patch 1: Gas background (fills entire domain) ===== - "patch_icpp(1)%geometry": 9, + # patch 1: nitrogen background filling the whole domain + "patch_icpp(1)%geometry": 3, # rectangle (2D) "patch_icpp(1)%x_centroid": 0.0, - "patch_icpp(1)%y_centroid": 0.0, - "patch_icpp(1)%z_centroid": 0.0, - "patch_icpp(1)%length_x": x1 - x0, - "patch_icpp(1)%length_y": y1 - y0, - "patch_icpp(1)%length_z": z1 - z0, + "patch_icpp(1)%y_centroid": Ly_max / 2.0, + "patch_icpp(1)%length_x": 2.5 * Lx_half, # cover w/ margin + "patch_icpp(1)%length_y": 2.0 * Ly_max, "patch_icpp(1)%vel(1)": 0.0, "patch_icpp(1)%vel(2)": 0.0, - "patch_icpp(1)%vel(3)": 0.0, - "patch_icpp(1)%pres": p_gas, + "patch_icpp(1)%pres": p_amb, "patch_icpp(1)%alpha_rho(1)": eps * rho_l, "patch_icpp(1)%alpha_rho(2)": (1.0 - eps) * rho_g, "patch_icpp(1)%alpha(1)": eps, "patch_icpp(1)%alpha(2)": 1.0 - eps, - "patch_icpp(1)%cf_val": 0, - # ===== Patch 2: Left droplet (center at -sep, moves in +x) ===== - "patch_icpp(2)%geometry": 8, - "patch_icpp(2)%x_centroid": -sep, - "patch_icpp(2)%y_centroid": 0.0, - "patch_icpp(2)%z_centroid": 0.0, - "patch_icpp(2)%radius": R, + "patch_icpp(1)%cf_val": 0.0, + # patch 2: left droplet, moving in +x "patch_icpp(2)%alter_patch(1)": "T", "patch_icpp(2)%smoothen": "T", "patch_icpp(2)%smooth_patch_id": 1, - "patch_icpp(2)%smooth_coeff": 0.95, - "patch_icpp(2)%vel(1)": Ud, + "patch_icpp(2)%smooth_coeff": 0.5, + "patch_icpp(2)%geometry": 2, # circle (axisymmetric disk becomes a sphere) + "patch_icpp(2)%x_centroid": -x_center, + "patch_icpp(2)%y_centroid": 0.0, + "patch_icpp(2)%radius": R, + "patch_icpp(2)%vel(1)": +U, "patch_icpp(2)%vel(2)": 0.0, - "patch_icpp(2)%vel(3)": 0.0, - "patch_icpp(2)%pres": p_liq, + "patch_icpp(2)%pres": p_amb, "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_l, "patch_icpp(2)%alpha_rho(2)": eps * rho_g, "patch_icpp(2)%alpha(1)": 1.0 - eps, "patch_icpp(2)%alpha(2)": eps, - "patch_icpp(2)%cf_val": 1, - # ===== Patch 3: Right droplet (center at +sep, moves in -x) ===== - "patch_icpp(3)%geometry": 8, - "patch_icpp(3)%x_centroid": sep, - "patch_icpp(3)%y_centroid": 0.0, - "patch_icpp(3)%z_centroid": 0.0, - "patch_icpp(3)%radius": R, + "patch_icpp(2)%cf_val": 1.0, + # patch 3: right droplet, moving in -x "patch_icpp(3)%alter_patch(1)": "T", "patch_icpp(3)%smoothen": "T", "patch_icpp(3)%smooth_patch_id": 1, - "patch_icpp(3)%smooth_coeff": 0.95, - "patch_icpp(3)%vel(1)": -Ud, + "patch_icpp(3)%smooth_coeff": 0.5, + "patch_icpp(3)%geometry": 2, + "patch_icpp(3)%x_centroid": +x_center, + "patch_icpp(3)%y_centroid": 0.0, + "patch_icpp(3)%radius": R, + "patch_icpp(3)%vel(1)": -U, "patch_icpp(3)%vel(2)": 0.0, - "patch_icpp(3)%vel(3)": 0.0, - "patch_icpp(3)%pres": p_liq, + "patch_icpp(3)%pres": p_amb, "patch_icpp(3)%alpha_rho(1)": (1.0 - eps) * rho_l, "patch_icpp(3)%alpha_rho(2)": eps * rho_g, "patch_icpp(3)%alpha(1)": 1.0 - eps, "patch_icpp(3)%alpha(2)": eps, - "patch_icpp(3)%cf_val": 1, + "patch_icpp(3)%cf_val": 1.0, } -print(json.dumps(data)) +# Derived-quantity banner (printed to stderr on --case-optimization runs) +# Printed here as a comment for humans; MFC only reads the JSON below. +# U_rel = 1.8823 m/s (We = 30, head-on: U_rel = 2U) +# U = 0.9412 m/s per drop +# h_c = (A R / sigma)^(1/3) = 38 nm (A = 1e-20 J) +# tau_i = D/U_rel ≈ 1.59e-4 s +# tau_s = sqrt(rho D^3/sigma) ≈ 8.73e-4 s +# Re = rho U_rel D / mu ≈ 205 +# Nx = N_fine + 2*N_stretch (defaults: 400 + 96 = 496) +# Ny = N_fine_r + N_stretch_r (defaults: 800 + 48 = 848) + +print(json.dumps(case)) diff --git a/examples/3D_droplet_coalescence/case_muscl.py b/examples/3D_droplet_coalescence/case_muscl.py new file mode 100644 index 0000000000..54e215f32d --- /dev/null +++ b/examples/3D_droplet_coalescence/case_muscl.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 + +import argparse +import json +import math + +# Parameter table: (We, Re, B) +CASES = { + "a": (0.2, 14.8, 0.20), + "b": (0.5, 23.6, 0.10), + "c": (8.6, 105.9, 0.08), + "d": (15.2, 139.8, 0.08), + "e": (19.4, 158.0, 0.05), + "f": (32.8, 210.8, 0.08), + "g": (37.2, 228.0, 0.01), + "h": (61.4, 296.5, 0.06), + "i": (61.3, 295.3, 0.11), + "j": (56.3, 288.9, 0.13), + "k": (70.8, 327.7, 0.25), + "l": (48.1, 270.1, 0.39), + "m": (60.1, 302.8, 0.55), + "n": (65.1, 320.3, 0.49), + "o": (60.8, 313.7, 0.68), + "p": (64.9, 312.8, 0.71), + "q": (48.8, 260.3, 0.72), + "r": (14.5, 149.1, 0.34), +} + +parser = argparse.ArgumentParser() +parser.add_argument("--mfc", type=json.loads, default="{}") +parser.add_argument("--case", type=str, default="a", choices=CASES.keys()) +args = parser.parse_args() + +We, Re, B = CASES[args.case] + +# Physical properties (tetradecane / nitrogen) — SI units +rho_l = 763.0 # kg/m^3 +rho_g = rho_l / 666 # kg/m^3 +sigma = 0.0266 # N/m +gamma_l = 2.35 +gamma_g = 1.4 + +D = 240e-6 # m +R = D / 2.0 + +# Pressure and EOS parameters scaled from non-dimensional formulation +# (low effective pressure keeps the speed of sound manageable) +p_ref = sigma / (0.72 * D) +p_gas = p_ref +pi_inf_l = 190.0 * p_ref +pi_inf_g = 0.0 +p_liq = p_gas + 2 * sigma / R + +# Collision parameters from dimensionless numbers +Ur = math.sqrt(We * sigma / (rho_l * D)) +Ud = Ur / 2.0 + +mu_l = rho_l * Ur * D / Re +mu_g = mu_l / 119 + +# Domain (m) +x0, x1 = -2.25 * D, 2.25 * D +y0, y1 = -1.5 * D, 1.5 * D +z0, z1 = -1.5 * D, 1.5 * D + +Nx = 199 +Ny = 139 +Nz = 139 + +eps = 1e-9 + +sep = 0.5 * D +b = B * D + +# Simulation time +t_end = 1e-3 # seconds +t_save = t_end / 1000 + +data = { + "run_time_info": "T", + "x_domain%beg": x0, + "x_domain%end": x1, + "y_domain%beg": y0, + "y_domain%end": y1, + "z_domain%beg": z0, + "z_domain%end": z1, + "m": Nx, + "n": Ny, + "p": Nz, + "cyl_coord": "F", + "cfl_adap_dt": "T", + "cfl_target": 0.1, + "n_start": 0, + "t_stop": t_end, + "t_save": t_save, + "model_eqns": 2, + "alt_soundspeed": "F", + "low_Mach": 0, + "mixture_err": "T", + "mpp_lim": "T", + "time_stepper": 3, + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": 4, + "int_comp": "T", + "ic_beta": 2.0, + "avg_state": 2, + "riemann_solver": 2, + "wave_speeds": 1, + "viscous": "T", + "bc_x%beg": -6, + "bc_x%end": -6, + "bc_y%beg": -6, + "bc_y%end": -6, + "bc_z%beg": -6, + "bc_z%end": -6, + "num_patches": 3, + "num_fluids": 2, + "surface_tension": "T", + "format": 1, + "precision": 2, + "alpha_wrt(1)": "T", + "cf_wrt": "T", + "vel_wrt(1)": "T", + "vel_wrt(2)": "T", + "vel_wrt(3)": "T", + "pres_wrt": "T", + "parallel_io": "T", + "sigma": sigma, + "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), + "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), + "fluid_pp(1)%Re(1)": 1.0 / mu_l, + "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), + "fluid_pp(2)%pi_inf": 0.0, + "fluid_pp(2)%Re(1)": 1.0 / mu_g, + # Patch 1: Gas background + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%z_centroid": 0.0, + "patch_icpp(1)%length_x": x1 - x0, + "patch_icpp(1)%length_y": y1 - y0, + "patch_icpp(1)%length_z": z1 - z0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": p_gas, + "patch_icpp(1)%alpha_rho(1)": eps * rho_l, + "patch_icpp(1)%alpha_rho(2)": (1.0 - eps) * rho_g, + "patch_icpp(1)%alpha(1)": eps, + "patch_icpp(1)%alpha(2)": 1.0 - eps, + "patch_icpp(1)%cf_val": 0, + # Patch 2: Left droplet + "patch_icpp(2)%geometry": 8, + "patch_icpp(2)%x_centroid": -sep, + "patch_icpp(2)%y_centroid": -b / 2.0, + "patch_icpp(2)%z_centroid": 0.0, + "patch_icpp(2)%radius": R, + "patch_icpp(2)%alter_patch(1)": "T", + "patch_icpp(2)%smoothen": "T", + "patch_icpp(2)%smooth_patch_id": 1, + "patch_icpp(2)%smooth_coeff": 0.95, + "patch_icpp(2)%vel(1)": Ud, + "patch_icpp(2)%vel(2)": 0.0, + "patch_icpp(2)%vel(3)": 0.0, + "patch_icpp(2)%pres": p_liq, + "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_l, + "patch_icpp(2)%alpha_rho(2)": eps * rho_g, + "patch_icpp(2)%alpha(1)": 1.0 - eps, + "patch_icpp(2)%alpha(2)": eps, + "patch_icpp(2)%cf_val": 1, + # Patch 3: Right droplet + "patch_icpp(3)%geometry": 8, + "patch_icpp(3)%x_centroid": sep, + "patch_icpp(3)%y_centroid": b / 2.0, + "patch_icpp(3)%z_centroid": 0.0, + "patch_icpp(3)%radius": R, + "patch_icpp(3)%alter_patch(1)": "T", + "patch_icpp(3)%smoothen": "T", + "patch_icpp(3)%smooth_patch_id": 1, + "patch_icpp(3)%smooth_coeff": 0.95, + "patch_icpp(3)%vel(1)": -Ud, + "patch_icpp(3)%vel(2)": 0.0, + "patch_icpp(3)%vel(3)": 0.0, + "patch_icpp(3)%pres": p_liq, + "patch_icpp(3)%alpha_rho(1)": (1.0 - eps) * rho_l, + "patch_icpp(3)%alpha_rho(2)": eps * rho_g, + "patch_icpp(3)%alpha(1)": 1.0 - eps, + "patch_icpp(3)%alpha(2)": eps, + "patch_icpp(3)%cf_val": 1, +} + +print(json.dumps(data)) From 7c206455c156deac24ea8048cc83673669e813b4 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Tue, 21 Apr 2026 21:11:58 -0400 Subject: [PATCH 04/10] format: apply ruff formatting to run.py --- toolchain/mfc/run/run.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/toolchain/mfc/run/run.py b/toolchain/mfc/run/run.py index f7a502a039..2b1be58c46 100644 --- a/toolchain/mfc/run/run.py +++ b/toolchain/mfc/run/run.py @@ -143,10 +143,9 @@ def __execute_job_script(qsystem: queues.QueueSystem): # For interactive runs on Unix, tee output to a .out file so users get # the same log file that batch runs produce via the scheduler. - if isinstance(qsystem, queues.InteractiveSystem) and os.name != 'nt': + if isinstance(qsystem, queues.InteractiveSystem) and os.name != "nt": out_filepath = os.path.abspath(os.path.join(cwd, f"{ARG('name')}.out")) - cmd = ["/bin/bash", "-c", - f"set -o pipefail; {shlex.join([str(x) for x in cmd])} 2>&1 | tee {shlex.quote(out_filepath)}"] + cmd = ["/bin/bash", "-c", f"set -o pipefail; {shlex.join([str(x) for x in cmd])} 2>&1 | tee {shlex.quote(out_filepath)}"] # Execute the job script with appropriate output handling # At verbosity >= 2, show print_cmd=True for system() calls From fe42aa45c9fee22517b6b4c09cbe9885553f4f04 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Tue, 21 Apr 2026 23:30:00 -0400 Subject: [PATCH 05/10] update --- .../2D_axisym_droplet_coalescence/case.py | 248 ++++-------- examples/3D_droplet_coalescence/case.py | 360 ++++++++---------- examples/3D_droplet_coalescence/case_muscl.py | 193 ---------- 3 files changed, 235 insertions(+), 566 deletions(-) delete mode 100644 examples/3D_droplet_coalescence/case_muscl.py diff --git a/examples/2D_axisym_droplet_coalescence/case.py b/examples/2D_axisym_droplet_coalescence/case.py index 5bb8c5a7f4..d01308b95f 100644 --- a/examples/2D_axisym_droplet_coalescence/case.py +++ b/examples/2D_axisym_droplet_coalescence/case.py @@ -10,40 +10,24 @@ # ./mfc.sh run case.py -n # single node # ./mfc.sh run case.py --gpu -n # GPU build # Knobs at top of file: -# SMOKE_TEST — True: fast end-to-end smoke test (no stretch, fixed dt, 100 steps) -# False: production (stretched grid, CFL-adaptive dt, t_stop = 200 us) # SCHEME — 'wenoz' or 'muscl_thinc' -# DX_FINE — cell size in drainage film (m) -# DX_BULK — cell size in bulk (m) +# DX_FINE — cell size in drainage film (m); start 10 nm +# DX_BULK — cell size in bulk (m); start 1 micron # FILM_HALF — axial half-width of the uniform fine zone (m) # WE — target Weber number (default 30) -# Known issues: -# - With SMOKE_TEST = False and stretching enabled, simulation hits NaN at -# corner cell (m, n) after ~10 steps. Cause under investigation. The -# smoke-test path (uniform grid) runs end-to-end cleanly. import json import math # User-facing knobs -SMOKE_TEST = True # True: coarse + uniform + 100 fixed-dt steps SCHEME = "wenoz" # 'wenoz' | 'muscl_thinc' +DX_FINE = 10e-9 # 10 nm — film cell size target +DX_BULK = 1e-6 # 1 um — bulk cell size target FILM_HALF = 2e-6 # ± half-width of the uniform fine zone on the axis (m) WE = 30.0 # Weber number (coalescence regime) -# Droplet diameter defined here so DX_FINE can reference it; also used below. -D = 300e-6 # droplet diameter, 300 um - -# DX_FINE = D/1000 = 300 nm for both modes. Sets the UNIFORM PRE-STRETCH cell -# size domain/N (see Nx/Ny derivation below). MFC's tanh stretch keeps cells -# inside [x_a, x_b] at ~this size and grows cells outward; effective bulk -# cell size is emergent, not a user knob. -DX_FINE = D / 1000.0 - -# Optional explicit overrides for (Nx, Ny). If set, DX_FINE becomes diagnostic. -NX_OVERRIDE = None -NY_OVERRIDE = None +# For quick debugging, set coarser values, e.g. DX_FINE=40e-9, DX_BULK=4e-6. # Physical constants (tetradecane / nitrogen at ~25 C) @@ -66,58 +50,60 @@ # Droplet geometry and kinematics -R = D / 2.0 # radius, 150 um (D = 300 um defined at top) +D = 300e-6 # droplet diameter, 300 um +R = D / 2.0 # radius, 150 um # Weber number: We = rho_l * U_rel^2 * D / sigma, U_rel = 2U per Qian & Law. +# U_rel^2 = We*sigma/(rho_l*D) +# = 30 * 0.027 / (762 * 3e-4) +# = 0.810 / 0.2286 +# = 3.5433 m^2/s^2 # U_rel = 1.8823 m/s -> U = 0.9412 m/s per droplet U_rel = math.sqrt(WE * sigma / (rho_l * D)) U = 0.5 * U_rel # Initial axial gap between droplet surfaces (pre-drainage), and centers. -gap0 = 20e-6 # 20 um initial gap +gap0 = 20e-6 # 20 um initial gap — leaves inertial approach + drainage x_center = R + 0.5 * gap0 # droplet centers at x = +/- x_center # Domain & grid +# Axial extent — big enough that outgoing waves don't reflect back onto the film. +# Keep a few D each side of the droplets. Lx_half = 5.0 * R # = 750 um axial half-domain Ly_max = 3.0 * R # = 450 um radial extent - -def _post_stretch_extent(lo, hi, a, xa, xb, loops): - """Mimic MFC's hyperbolic-tangent stretching in pre_process/m_grid.f90 to - compute the post-stretch domain extent. The formula uses normalized coords: - x <- x / length; xa <- xa / length; xb <- xb / length; - x <- x/a * (a + log(cosh(a*(x-xa))) + log(cosh(a*(x-xb))) - 2*log(cosh(a*(xb-xa)/2))) - x <- x * length - This expands the domain beyond [lo, hi]; patches must cover the result or - the outer cells are left uninitialized and the sim NaNs at step 1. - """ - length = abs(hi - lo) - x0, x1 = lo / length, hi / length - xa_n, xb_n = xa / length, xb / length - for _ in range(loops): - const = 2.0 * math.log(math.cosh(a * (xb_n - xa_n) / 2.0)) - - def f(x): - return x / a * (a + math.log(math.cosh(a * (x - xa_n))) + math.log(math.cosh(a * (x - xb_n))) - const) - - x0, x1 = f(x0), f(x1) - return x0 * length, x1 * length - - -# Cell count to hit DX_FINE in the fine zone. -# MFC's tanh stretch (m_grid.f90 s_generate_serial_grid) starts with N uniform -# cells of size dx = domain/N, then remaps their boundaries. The remap -# derivative is ~1 inside [x_a, x_b] (cells stay at ~dx) and grows outside -# (cells expand). So to achieve DX_FINE in the fine zone: -# N = ceil(domain / DX_FINE) -# Bulk cell sizes after stretch are set by (a_*, loops_*) and are reported in -# the stderr diagnostic at the bottom of this file. -Nx = NX_OVERRIDE if NX_OVERRIDE is not None else int(math.ceil(2.0 * Lx_half / DX_FINE)) -Ny = NY_OVERRIDE if NY_OVERRIDE is not None else int(math.ceil(Ly_max / DX_FINE)) +# Uniform fine zone in x: [-FILM_HALF, +FILM_HALF] at DX_FINE. +# Stretched zone: DX_FINE -> DX_BULK, geometric growth handled by MFC stretch_x. +# Cell-count estimate (for sanity; MFC's hyperbolic stretch won't be exactly geometric): +# N_fine = 2*FILM_HALF / DX_FINE +# N_stretch (each side) ~ ln(DX_BULK/DX_FINE)/ln(r) with r = 1.1 -> ~48 +N_fine = int(2 * FILM_HALF / DX_FINE) # 400 for defaults +N_stretch = int(math.ceil(math.log(DX_BULK / DX_FINE) / math.log(1.1))) # ~48 each side +Nx = N_fine + 2 * N_stretch # 496 for defaults + +# Radial: fine near axis (film rim), coarsens outward. +# Same strategy on y (positive only, because axisymmetric). +N_fine_r = int(4 * FILM_HALF / DX_FINE) # 800 cells out to 8 um +N_stretch_r = int(math.ceil(math.log(DX_BULK / DX_FINE) / math.log(1.1))) # ~48 +Ny = N_fine_r + N_stretch_r # 848 for defaults + +# NOTE: These are generous baselines for a grid-convergence study. +# On first debug runs, set DX_FINE=40e-9, DX_BULK=4e-6 to cut Nx,Ny by ~4x. + +# Time integration + +# Characteristic times: +# tau_inertia = D / U_rel = 3e-4 / 1.8823 ≈ 1.59e-4 s +# tau_cap = sqrt(rho_l D^3 / sigma) ≈ 8.73e-4 s +# Want to simulate until h_min ~ 100–200 nm. Experience says ~1–2 * tau_inertia. +t_stop = 2.0e-4 # 200 us +t_save = 5.0e-6 # 40 frames # Numerical scheme switches +# WENO-Z 5 is the default from MFC 5.0 paper; MUSCL+THINC is the fallback if +# interface width exceeds ~5 cells or if we see early numerical rupture. if SCHEME == "wenoz": recon = { "recon_type": 1, @@ -140,66 +126,6 @@ def f(x): else: raise ValueError(f"Unknown SCHEME: {SCHEME!r}") -# Time stepping & grid stretching (depends on SMOKE_TEST) - -if SMOKE_TEST: - # Inching roadmap (edit dt / t_step_stop to advance stages): - # stage 0 (old): dt=1e-12, n=100, save=10 -> 100 ps total, 13 s/run - # stage 1 (now): dt=1e-11, n=10000, save=1000 -> 100 ns total, ~21 min @-n1 - # stage 2: dt=1e-11, n=100000, save=10000 -> 1 µs total, ~3.5 hr @-n1 - # stage 3: switch to cfl_adap_dt + t_stop=5e-6 (see production block) - # stage 4: refine DX_FINE to 20 nm -> 10 nm - # stage 5: enable stretch_y (currently NaNs at step ~10 — debug first) - # CFL on 40-nm grid with c_liquid=1316 m/s: dt_max(CFL=0.3)=9.1e-12 s. - # dt=1e-11 sits at CFL~0.33 — safely inside stability. - time_block = { - "dt": 1.0e-11, - "t_step_start": 0, - "t_step_stop": 10000, - "t_step_save": 1000, - } - stretch_block = { - "stretch_x": "T", - "a_x": 2.0, - "x_a": -FILM_HALF, - "x_b": +FILM_HALF, - "loops_x": 1, - "stretch_y": "F", - } -else: - # Production: CFL-adaptive dt, stretched grid. - # t_stop = 2e-4 s (~1.25 * tau_inertia = 1.59e-4 s), 40 frames. - time_block = { - "cfl_adap_dt": "T", - "cfl_target": 0.3, - "t_stop": 2.0e-4, - "t_save": 5.0e-6, - "n_start": 0, - } - # Production: uniform 300 nm cells everywhere (Nx/Ny already large enough - # to hit DX_FINE without tanh concentration). Disabling stretch keeps the - # domain at the intended [-Lx_half, +Lx_half] x [0, Ly_max]; with a=4 - # loops=2 the tanh remap was expanding to ±2892 um x 7 mm and wasting - # cells on empty bulk. Re-enable if you need non-uniform spacing later. - stretch_block = { - "stretch_x": "F", - "stretch_y": "F", - } - -# Compute post-stretch domain extent so patches cover the actual cell range. -# Add 10% margin on top, to account for cell centers vs boundaries. -if stretch_block["stretch_x"] == "T": - _xlo, _xhi = _post_stretch_extent(-Lx_half, +Lx_half, stretch_block["a_x"], stretch_block["x_a"], stretch_block["x_b"], stretch_block["loops_x"]) - patch1_length_x = 1.10 * (_xhi - _xlo) -else: - patch1_length_x = 2.0 * Lx_half - -if stretch_block["stretch_y"] == "T": - _ylo, _yhi = _post_stretch_extent(0.0, Ly_max, stretch_block["a_y"], stretch_block["y_a"], stretch_block["y_b"], stretch_block["loops_y"]) - patch1_length_y = 1.10 * (_yhi - _ylo) -else: - patch1_length_y = Ly_max - # MFC input deck eps = 1e-9 # avoid exactly 0/1 volume fractions @@ -214,15 +140,30 @@ def f(x): "n": Ny - 1, "p": 0, # 2D "cyl_coord": "T", # 2D axisymmetric: x is axial, y is radial + # Grid stretching — fine uniform zone in [-FILM_HALF, +FILM_HALF], stretched outside. + "stretch_x": "T", + "a_x": 4.0, + "x_a": -FILM_HALF, + "x_b": +FILM_HALF, + "loops_x": 2, + "stretch_y": "T", + "a_y": 4.0, + "y_a": 0.0, + "y_b": 4.0 * FILM_HALF, # keep fine near the axis / film rim + "loops_y": 2, # boundary conditions # Axial (x): outflow/extrapolation on both ends. - "bc_x%beg": -6, # TRIAGE: char non-reflecting (was -3) - "bc_x%end": -6, + "bc_x%beg": -3, + "bc_x%end": -3, # Radial (y): axisymmetric at y=0, outflow at y_max. - "bc_y%beg": -2, # axis - "bc_y%end": -6, # TRIAGE: char non-reflecting (was -3) + "bc_y%beg": -2, # reflective / axis + "bc_y%end": -3, # time integration "time_stepper": 3, # TVD-RK3 + "cfl_adap_dt": "T", + "cfl_target": 0.3, # conservative; tighten/loosen after first runs + "t_stop": t_stop, + "t_save": t_save, # physical model "model_eqns": 2, # 5-equation (Allaire) diffuse interface "num_fluids": 2, @@ -232,6 +173,7 @@ def f(x): "avg_state": 2, "riemann_solver": 2, # HLLC "wave_speeds": 1, + **recon, # surface tension "surface_tension": "T", "sigma": sigma, @@ -241,11 +183,11 @@ def f(x): "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), # 1/(1.4-1) = 2.5 "fluid_pp(2)%pi_inf": 0.0, # patch 1: nitrogen background filling the whole domain - "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%geometry": 3, # rectangle (2D) "patch_icpp(1)%x_centroid": 0.0, - "patch_icpp(1)%y_centroid": patch1_length_y / 2.0, - "patch_icpp(1)%length_x": patch1_length_x, - "patch_icpp(1)%length_y": patch1_length_y, + "patch_icpp(1)%y_centroid": Ly_max / 2.0, + "patch_icpp(1)%length_x": 2.5 * Lx_half, # cover w/ margin + "patch_icpp(1)%length_y": 2.0 * Ly_max, "patch_icpp(1)%vel(1)": 0.0, "patch_icpp(1)%vel(2)": 0.0, "patch_icpp(1)%pres": p_amb, @@ -258,7 +200,7 @@ def f(x): "patch_icpp(2)%alter_patch(1)": "T", "patch_icpp(2)%smoothen": "T", "patch_icpp(2)%smooth_patch_id": 1, - "patch_icpp(2)%smooth_coeff": 0.95, + "patch_icpp(2)%smooth_coeff": 0.5, "patch_icpp(2)%geometry": 2, # circle (axisymmetric disk becomes a sphere) "patch_icpp(2)%x_centroid": -x_center, "patch_icpp(2)%y_centroid": 0.0, @@ -275,7 +217,7 @@ def f(x): "patch_icpp(3)%alter_patch(1)": "T", "patch_icpp(3)%smoothen": "T", "patch_icpp(3)%smooth_patch_id": 1, - "patch_icpp(3)%smooth_coeff": 0.95, + "patch_icpp(3)%smooth_coeff": 0.5, "patch_icpp(3)%geometry": 2, "patch_icpp(3)%x_centroid": +x_center, "patch_icpp(3)%y_centroid": 0.0, @@ -288,51 +230,17 @@ def f(x): "patch_icpp(3)%alpha(1)": 1.0 - eps, "patch_icpp(3)%alpha(2)": eps, "patch_icpp(3)%cf_val": 1.0, - # post processing - "alpha_wrt(1)": "T", # liquid volume fraction (the field h_min lives in) - "alpha_wrt(2)": "T", # gas volume fraction - "rho_wrt": "T", - "pres_wrt": "T", - "mom_wrt(1)": "T", - "mom_wrt(2)": "T", - "cf_wrt": "T", # color function (requires surface_tension) - "parallel_io": "T", - "format": 1, # 1 = silo-HDF5; 2 = binary } -case.update(recon) -case.update(time_block) -case.update(stretch_block) - -# Diagnostic: report ACHIEVED cell sizes after MFC's tanh stretch (stderr so -# it doesn't pollute the JSON on stdout). -import sys as _sys - - -def _stretched_cb(m, beg, end, stretch, a, xa, xb, loops): - cb = [beg + (end - beg) * i / (m + 1) for i in range(-1, m + 1)] # m+2 bnds - if stretch: - length = abs(cb[-1] - cb[0]) - v = [x / length for x in cb] - xa_n, xb_n = xa / length, xb / length - for _ in range(loops): - c = 2.0 * math.log(math.cosh(a * (xb_n - xa_n) / 2.0)) - v = [x / a * (a + math.log(math.cosh(a * (x - xa_n))) + math.log(math.cosh(a * (x - xb_n))) - c) for x in v] - cb = [x * length for x in v] - return cb - - -def _report(axis, m, beg, end): - sk_stretch = case.get(f"stretch_{axis}", "F") == "T" - cb = _stretched_cb(m, beg, end, sk_stretch, case.get(f"a_{axis}", 1.0), case.get(f"{axis}_a", 0.0), case.get(f"{axis}_b", 0.0), case.get(f"loops_{axis}", 0)) - cells = [cb[i + 1] - cb[i] for i in range(m + 1)] - _sys.stderr.write(f"[case.py] {axis}: N={m + 1:6d} min={min(cells) * 1e9:>9.1f} nm max={max(cells) * 1e6:>8.2f} um span=[{cb[0] * 1e6:8.1f},{cb[-1] * 1e6:8.1f}] um\n") - - -_sys.stderr.write(f"[case.py] SMOKE_TEST={SMOKE_TEST} DX_FINE target={DX_FINE * 1e9:.0f} nm\n") -_report("x", Nx - 1, -Lx_half, +Lx_half) -_report("y", Ny - 1, 0.0, +Ly_max) -_sys.stderr.write(f"[case.py] total cells = {Nx * Ny:,}\n") -_sys.stderr.flush() +# Derived-quantity banner (printed to stderr on --case-optimization runs) +# Printed here as a comment for humans; MFC only reads the JSON below. +# U_rel = 1.8823 m/s (We = 30, head-on: U_rel = 2U) +# U = 0.9412 m/s per drop +# h_c = (A R / sigma)^(1/3) = 38 nm (A = 1e-20 J) +# tau_i = D/U_rel ≈ 1.59e-4 s +# tau_s = sqrt(rho D^3/sigma) ≈ 8.73e-4 s +# Re = rho U_rel D / mu ≈ 205 +# Nx = N_fine + 2*N_stretch (defaults: 400 + 96 = 496) +# Ny = N_fine_r + N_stretch_r (defaults: 800 + 48 = 848) print(json.dumps(case)) diff --git a/examples/3D_droplet_coalescence/case.py b/examples/3D_droplet_coalescence/case.py index d01308b95f..c0cab686ba 100644 --- a/examples/3D_droplet_coalescence/case.py +++ b/examples/3D_droplet_coalescence/case.py @@ -1,246 +1,200 @@ #!/usr/bin/env python3 -# Phase 1 drainage validation case — acoustic_demulsification project -# Geometry: 2D axisymmetric, head-on collision of two equal droplets -# Fluids : tetradecane drops in nitrogen ambient -# Target : resolve gas-film drainage; compare h_min(t) to Klaseboer 2000 -# Stop : when h_min reaches ~100–200 nm (before numerical rupture) -# Roadmap reference: +Projects/acoustic_demulsification/simulation_roadmap.md -# (Phase 1 — Drainage validation) -# Run: -# ./mfc.sh run case.py -n # single node -# ./mfc.sh run case.py --gpu -n # GPU build -# Knobs at top of file: -# SCHEME — 'wenoz' or 'muscl_thinc' -# DX_FINE — cell size in drainage film (m); start 10 nm -# DX_BULK — cell size in bulk (m); start 1 micron -# FILM_HALF — axial half-width of the uniform fine zone (m) -# WE — target Weber number (default 30) +import argparse import json import math -# User-facing knobs - -SCHEME = "wenoz" # 'wenoz' | 'muscl_thinc' -DX_FINE = 10e-9 # 10 nm — film cell size target -DX_BULK = 1e-6 # 1 um — bulk cell size target -FILM_HALF = 2e-6 # ± half-width of the uniform fine zone on the axis (m) -WE = 30.0 # Weber number (coalescence regime) +# Parameter table: (We, Re, B) +CASES = { + "a": (0.2, 14.8, 0.20), + "b": (0.5, 23.6, 0.10), + "c": (8.6, 105.9, 0.08), + "d": (15.2, 139.8, 0.08), + "e": (19.4, 158.0, 0.05), + "f": (32.8, 210.8, 0.08), + "g": (37.2, 228.0, 0.01), + "h": (61.4, 296.5, 0.06), + "i": (61.3, 295.3, 0.11), + "j": (56.3, 288.9, 0.13), + "k": (70.8, 327.7, 0.25), + "l": (48.1, 270.1, 0.39), + "m": (60.1, 302.8, 0.55), + "n": (65.1, 320.3, 0.49), + "o": (60.8, 313.7, 0.68), + "p": (64.9, 312.8, 0.71), + "q": (48.8, 260.3, 0.72), + "r": (14.5, 149.1, 0.34), +} -# For quick debugging, set coarser values, e.g. DX_FINE=40e-9, DX_BULK=4e-6. +parser = argparse.ArgumentParser() +parser.add_argument("--mfc", type=json.loads, default="{}") +parser.add_argument("--case", type=str, default="a", choices=CASES.keys()) +args = parser.parse_args() -# Physical constants (tetradecane / nitrogen at ~25 C) +We, Re, B = CASES[args.case] -# Tetradecane (liquid) -rho_l = 762.0 # kg/m^3 -mu_l = 2.1e-3 # Pa*s -sigma = 0.027 # N/m (≈ 27 mN/m) +# Physical properties (tetradecane / nitrogen) — SI units +rho_l = 763.0 # kg/m^3 +rho_g = rho_l / 666 # kg/m^3 +sigma = 0.0266 # N/m +gamma_l = 2.35 +gamma_g = 1.4 -# Stiffened-gas fit for tetradecane: gamma = 4.4, p_inf = 3.0e8 Pa -# => c_l = sqrt(gamma*(p+p_inf)/rho) ≈ 1316 m/s at 1 atm — matches bulk data. -gamma_l = 4.4 -pi_inf_l = 3.0e8 +D = 240e-6 # m +R = D / 2.0 -# Nitrogen (ideal gas) -rho_g = 1.165 # kg/m^3 at 1 atm, 20 C -gamma_g = 1.4 +# Pressure and EOS parameters scaled from non-dimensional formulation +# (low effective pressure keeps the speed of sound manageable) +p_ref = sigma / (0.72 * D) +p_gas = p_ref +pi_inf_l = 190.0 * p_ref pi_inf_g = 0.0 - -p_amb = 1.01325e5 # 1 atm - -# Droplet geometry and kinematics - -D = 300e-6 # droplet diameter, 300 um -R = D / 2.0 # radius, 150 um - -# Weber number: We = rho_l * U_rel^2 * D / sigma, U_rel = 2U per Qian & Law. -# U_rel^2 = We*sigma/(rho_l*D) -# = 30 * 0.027 / (762 * 3e-4) -# = 0.810 / 0.2286 -# = 3.5433 m^2/s^2 -# U_rel = 1.8823 m/s -> U = 0.9412 m/s per droplet -U_rel = math.sqrt(WE * sigma / (rho_l * D)) -U = 0.5 * U_rel - -# Initial axial gap between droplet surfaces (pre-drainage), and centers. -gap0 = 20e-6 # 20 um initial gap — leaves inertial approach + drainage -x_center = R + 0.5 * gap0 # droplet centers at x = +/- x_center - -# Domain & grid - -# Axial extent — big enough that outgoing waves don't reflect back onto the film. -# Keep a few D each side of the droplets. -Lx_half = 5.0 * R # = 750 um axial half-domain -Ly_max = 3.0 * R # = 450 um radial extent - -# Uniform fine zone in x: [-FILM_HALF, +FILM_HALF] at DX_FINE. -# Stretched zone: DX_FINE -> DX_BULK, geometric growth handled by MFC stretch_x. -# Cell-count estimate (for sanity; MFC's hyperbolic stretch won't be exactly geometric): -# N_fine = 2*FILM_HALF / DX_FINE -# N_stretch (each side) ~ ln(DX_BULK/DX_FINE)/ln(r) with r = 1.1 -> ~48 -N_fine = int(2 * FILM_HALF / DX_FINE) # 400 for defaults -N_stretch = int(math.ceil(math.log(DX_BULK / DX_FINE) / math.log(1.1))) # ~48 each side -Nx = N_fine + 2 * N_stretch # 496 for defaults - -# Radial: fine near axis (film rim), coarsens outward. -# Same strategy on y (positive only, because axisymmetric). -N_fine_r = int(4 * FILM_HALF / DX_FINE) # 800 cells out to 8 um -N_stretch_r = int(math.ceil(math.log(DX_BULK / DX_FINE) / math.log(1.1))) # ~48 -Ny = N_fine_r + N_stretch_r # 848 for defaults - -# NOTE: These are generous baselines for a grid-convergence study. -# On first debug runs, set DX_FINE=40e-9, DX_BULK=4e-6 to cut Nx,Ny by ~4x. - -# Time integration - -# Characteristic times: -# tau_inertia = D / U_rel = 3e-4 / 1.8823 ≈ 1.59e-4 s -# tau_cap = sqrt(rho_l D^3 / sigma) ≈ 8.73e-4 s -# Want to simulate until h_min ~ 100–200 nm. Experience says ~1–2 * tau_inertia. -t_stop = 2.0e-4 # 200 us -t_save = 5.0e-6 # 40 frames - -# Numerical scheme switches - -# WENO-Z 5 is the default from MFC 5.0 paper; MUSCL+THINC is the fallback if -# interface width exceeds ~5 cells or if we see early numerical rupture. -if SCHEME == "wenoz": - recon = { - "recon_type": 1, - "weno_order": 5, - "weno_eps": 1e-6, - "mapped_weno": "F", - "wenoz": "T", - "mp_weno": "T", - "teno": "F", - } -elif SCHEME == "muscl_thinc": - recon = { - "recon_type": 2, - "muscl_order": 2, - "muscl_lim": 4, # Van Leer - "int_comp": "T", # THINC interface compression - "ic_eps": 1e-4, - "ic_beta": 1.6, - } -else: - raise ValueError(f"Unknown SCHEME: {SCHEME!r}") - -# MFC input deck - -eps = 1e-9 # avoid exactly 0/1 volume fractions - -case = { - # grid - "x_domain%beg": -Lx_half, - "x_domain%end": +Lx_half, - "y_domain%beg": 0.0, - "y_domain%end": +Ly_max, - "m": Nx - 1, # MFC uses zero-based cell-count (m = Nx-1) - "n": Ny - 1, - "p": 0, # 2D - "cyl_coord": "T", # 2D axisymmetric: x is axial, y is radial - # Grid stretching — fine uniform zone in [-FILM_HALF, +FILM_HALF], stretched outside. - "stretch_x": "T", - "a_x": 4.0, - "x_a": -FILM_HALF, - "x_b": +FILM_HALF, - "loops_x": 2, - "stretch_y": "T", - "a_y": 4.0, - "y_a": 0.0, - "y_b": 4.0 * FILM_HALF, # keep fine near the axis / film rim - "loops_y": 2, - # boundary conditions - # Axial (x): outflow/extrapolation on both ends. - "bc_x%beg": -3, - "bc_x%end": -3, - # Radial (y): axisymmetric at y=0, outflow at y_max. - "bc_y%beg": -2, # reflective / axis - "bc_y%end": -3, - # time integration - "time_stepper": 3, # TVD-RK3 - "cfl_adap_dt": "T", - "cfl_target": 0.3, # conservative; tighten/loosen after first runs - "t_stop": t_stop, - "t_save": t_save, - # physical model - "model_eqns": 2, # 5-equation (Allaire) diffuse interface - "num_fluids": 2, - "num_patches": 3, +p_liq = p_gas + 2 * sigma / R + +# Collision parameters from dimensionless numbers +Ur = math.sqrt(We * sigma / (rho_l * D)) +Ud = Ur / 2.0 + +mu_l = rho_l * Ur * D / Re +mu_g = mu_l / 119 + +# Domain (m) +x0, x1 = -2.25 * D, 2.25 * D +y0, y1 = -1.5 * D, 1.5 * D +z0, z1 = -1.5 * D, 1.5 * D + +Nx = 399 +Ny = 279 +Nz = 279 + +eps = 1e-9 + +sep = 0.55 * D +b = B * D + +# Simulation time +t_end = 1e-3 # seconds +t_save = t_end / 1000 + + +data = { + "run_time_info": "T", + "x_domain%beg": x0, + "x_domain%end": x1, + "y_domain%beg": y0, + "y_domain%end": y1, + "z_domain%beg": z0, + "z_domain%end": z1, + "m": Nx, + "n": Ny, + "p": Nz, + "cyl_coord": "F", + "n_start": 0, + # adaptive timestepping (disabled) + # "cfl_adap_dt": "T", + # "cfl_target": 0.1, + # "t_stop": t_end, + # "t_save": t_save, + # fixed timestepping + "dt": 1e-7, + "t_step_start": 0, + "t_step_stop": 1, + "t_step_save": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "low_Mach": 0, "mixture_err": "T", "mpp_lim": "T", + "time_stepper": 3, + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": 4, + "int_comp": "T", + "ic_beta": 1.6, "avg_state": 2, - "riemann_solver": 2, # HLLC + "riemann_solver": 2, "wave_speeds": 1, - **recon, - # surface tension + "viscous": "T", + "bc_x%beg": -6, + "bc_x%end": -6, + "bc_y%beg": -6, + "bc_y%end": -6, + "bc_z%beg": -6, + "bc_z%end": -6, + "num_patches": 3, + "num_fluids": 2, "surface_tension": "T", + "format": 1, + "precision": 2, + "alpha_wrt(1)": "T", + "cf_wrt": "T", + "vel_wrt(1)": "T", + "vel_wrt(2)": "T", + "vel_wrt(3)": "T", + "pres_wrt": "T", + "parallel_io": "T", "sigma": sigma, - # fluids (index 1 = tetradecane, 2 = nitrogen) - "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), # 1/(4.4-1) = 0.2941 - "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), # 4.4*3e8/3.4 ≈ 3.882e8 - "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), # 1/(1.4-1) = 2.5 + "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), + "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), + "fluid_pp(1)%Re(1)": 1.0 / mu_l, + "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), "fluid_pp(2)%pi_inf": 0.0, - # patch 1: nitrogen background filling the whole domain - "patch_icpp(1)%geometry": 3, # rectangle (2D) + "fluid_pp(2)%Re(1)": 1.0 / mu_g, + # Patch 1: Gas background + "patch_icpp(1)%geometry": 9, "patch_icpp(1)%x_centroid": 0.0, - "patch_icpp(1)%y_centroid": Ly_max / 2.0, - "patch_icpp(1)%length_x": 2.5 * Lx_half, # cover w/ margin - "patch_icpp(1)%length_y": 2.0 * Ly_max, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%z_centroid": 0.0, + "patch_icpp(1)%length_x": x1 - x0, + "patch_icpp(1)%length_y": y1 - y0, + "patch_icpp(1)%length_z": z1 - z0, "patch_icpp(1)%vel(1)": 0.0, "patch_icpp(1)%vel(2)": 0.0, - "patch_icpp(1)%pres": p_amb, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": p_gas, "patch_icpp(1)%alpha_rho(1)": eps * rho_l, "patch_icpp(1)%alpha_rho(2)": (1.0 - eps) * rho_g, "patch_icpp(1)%alpha(1)": eps, "patch_icpp(1)%alpha(2)": 1.0 - eps, - "patch_icpp(1)%cf_val": 0.0, - # patch 2: left droplet, moving in +x + "patch_icpp(1)%cf_val": 0, + # Patch 2: Left droplet + "patch_icpp(2)%geometry": 8, + "patch_icpp(2)%x_centroid": -sep, + "patch_icpp(2)%y_centroid": -b / 2.0, + "patch_icpp(2)%z_centroid": 0.0, + "patch_icpp(2)%radius": R, "patch_icpp(2)%alter_patch(1)": "T", "patch_icpp(2)%smoothen": "T", "patch_icpp(2)%smooth_patch_id": 1, - "patch_icpp(2)%smooth_coeff": 0.5, - "patch_icpp(2)%geometry": 2, # circle (axisymmetric disk becomes a sphere) - "patch_icpp(2)%x_centroid": -x_center, - "patch_icpp(2)%y_centroid": 0.0, - "patch_icpp(2)%radius": R, - "patch_icpp(2)%vel(1)": +U, + "patch_icpp(2)%smooth_coeff": 0.95, + "patch_icpp(2)%vel(1)": Ud, "patch_icpp(2)%vel(2)": 0.0, - "patch_icpp(2)%pres": p_amb, + "patch_icpp(2)%vel(3)": 0.0, + "patch_icpp(2)%pres": p_liq, "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_l, "patch_icpp(2)%alpha_rho(2)": eps * rho_g, "patch_icpp(2)%alpha(1)": 1.0 - eps, "patch_icpp(2)%alpha(2)": eps, - "patch_icpp(2)%cf_val": 1.0, - # patch 3: right droplet, moving in -x + "patch_icpp(2)%cf_val": 1, + # Patch 3: Right droplet + "patch_icpp(3)%geometry": 8, + "patch_icpp(3)%x_centroid": sep, + "patch_icpp(3)%y_centroid": b / 2.0, + "patch_icpp(3)%z_centroid": 0.0, + "patch_icpp(3)%radius": R, "patch_icpp(3)%alter_patch(1)": "T", "patch_icpp(3)%smoothen": "T", "patch_icpp(3)%smooth_patch_id": 1, - "patch_icpp(3)%smooth_coeff": 0.5, - "patch_icpp(3)%geometry": 2, - "patch_icpp(3)%x_centroid": +x_center, - "patch_icpp(3)%y_centroid": 0.0, - "patch_icpp(3)%radius": R, - "patch_icpp(3)%vel(1)": -U, + "patch_icpp(3)%smooth_coeff": 0.95, + "patch_icpp(3)%vel(1)": -Ud, "patch_icpp(3)%vel(2)": 0.0, - "patch_icpp(3)%pres": p_amb, + "patch_icpp(3)%vel(3)": 0.0, + "patch_icpp(3)%pres": p_liq, "patch_icpp(3)%alpha_rho(1)": (1.0 - eps) * rho_l, "patch_icpp(3)%alpha_rho(2)": eps * rho_g, "patch_icpp(3)%alpha(1)": 1.0 - eps, "patch_icpp(3)%alpha(2)": eps, - "patch_icpp(3)%cf_val": 1.0, + "patch_icpp(3)%cf_val": 1, } -# Derived-quantity banner (printed to stderr on --case-optimization runs) -# Printed here as a comment for humans; MFC only reads the JSON below. -# U_rel = 1.8823 m/s (We = 30, head-on: U_rel = 2U) -# U = 0.9412 m/s per drop -# h_c = (A R / sigma)^(1/3) = 38 nm (A = 1e-20 J) -# tau_i = D/U_rel ≈ 1.59e-4 s -# tau_s = sqrt(rho D^3/sigma) ≈ 8.73e-4 s -# Re = rho U_rel D / mu ≈ 205 -# Nx = N_fine + 2*N_stretch (defaults: 400 + 96 = 496) -# Ny = N_fine_r + N_stretch_r (defaults: 800 + 48 = 848) - -print(json.dumps(case)) +print(json.dumps(data)) diff --git a/examples/3D_droplet_coalescence/case_muscl.py b/examples/3D_droplet_coalescence/case_muscl.py deleted file mode 100644 index 54e215f32d..0000000000 --- a/examples/3D_droplet_coalescence/case_muscl.py +++ /dev/null @@ -1,193 +0,0 @@ -#!/usr/bin/env python3 - -import argparse -import json -import math - -# Parameter table: (We, Re, B) -CASES = { - "a": (0.2, 14.8, 0.20), - "b": (0.5, 23.6, 0.10), - "c": (8.6, 105.9, 0.08), - "d": (15.2, 139.8, 0.08), - "e": (19.4, 158.0, 0.05), - "f": (32.8, 210.8, 0.08), - "g": (37.2, 228.0, 0.01), - "h": (61.4, 296.5, 0.06), - "i": (61.3, 295.3, 0.11), - "j": (56.3, 288.9, 0.13), - "k": (70.8, 327.7, 0.25), - "l": (48.1, 270.1, 0.39), - "m": (60.1, 302.8, 0.55), - "n": (65.1, 320.3, 0.49), - "o": (60.8, 313.7, 0.68), - "p": (64.9, 312.8, 0.71), - "q": (48.8, 260.3, 0.72), - "r": (14.5, 149.1, 0.34), -} - -parser = argparse.ArgumentParser() -parser.add_argument("--mfc", type=json.loads, default="{}") -parser.add_argument("--case", type=str, default="a", choices=CASES.keys()) -args = parser.parse_args() - -We, Re, B = CASES[args.case] - -# Physical properties (tetradecane / nitrogen) — SI units -rho_l = 763.0 # kg/m^3 -rho_g = rho_l / 666 # kg/m^3 -sigma = 0.0266 # N/m -gamma_l = 2.35 -gamma_g = 1.4 - -D = 240e-6 # m -R = D / 2.0 - -# Pressure and EOS parameters scaled from non-dimensional formulation -# (low effective pressure keeps the speed of sound manageable) -p_ref = sigma / (0.72 * D) -p_gas = p_ref -pi_inf_l = 190.0 * p_ref -pi_inf_g = 0.0 -p_liq = p_gas + 2 * sigma / R - -# Collision parameters from dimensionless numbers -Ur = math.sqrt(We * sigma / (rho_l * D)) -Ud = Ur / 2.0 - -mu_l = rho_l * Ur * D / Re -mu_g = mu_l / 119 - -# Domain (m) -x0, x1 = -2.25 * D, 2.25 * D -y0, y1 = -1.5 * D, 1.5 * D -z0, z1 = -1.5 * D, 1.5 * D - -Nx = 199 -Ny = 139 -Nz = 139 - -eps = 1e-9 - -sep = 0.5 * D -b = B * D - -# Simulation time -t_end = 1e-3 # seconds -t_save = t_end / 1000 - -data = { - "run_time_info": "T", - "x_domain%beg": x0, - "x_domain%end": x1, - "y_domain%beg": y0, - "y_domain%end": y1, - "z_domain%beg": z0, - "z_domain%end": z1, - "m": Nx, - "n": Ny, - "p": Nz, - "cyl_coord": "F", - "cfl_adap_dt": "T", - "cfl_target": 0.1, - "n_start": 0, - "t_stop": t_end, - "t_save": t_save, - "model_eqns": 2, - "alt_soundspeed": "F", - "low_Mach": 0, - "mixture_err": "T", - "mpp_lim": "T", - "time_stepper": 3, - "recon_type": 2, - "muscl_order": 2, - "muscl_lim": 4, - "int_comp": "T", - "ic_beta": 2.0, - "avg_state": 2, - "riemann_solver": 2, - "wave_speeds": 1, - "viscous": "T", - "bc_x%beg": -6, - "bc_x%end": -6, - "bc_y%beg": -6, - "bc_y%end": -6, - "bc_z%beg": -6, - "bc_z%end": -6, - "num_patches": 3, - "num_fluids": 2, - "surface_tension": "T", - "format": 1, - "precision": 2, - "alpha_wrt(1)": "T", - "cf_wrt": "T", - "vel_wrt(1)": "T", - "vel_wrt(2)": "T", - "vel_wrt(3)": "T", - "pres_wrt": "T", - "parallel_io": "T", - "sigma": sigma, - "fluid_pp(1)%gamma": 1.0 / (gamma_l - 1.0), - "fluid_pp(1)%pi_inf": gamma_l * pi_inf_l / (gamma_l - 1.0), - "fluid_pp(1)%Re(1)": 1.0 / mu_l, - "fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0), - "fluid_pp(2)%pi_inf": 0.0, - "fluid_pp(2)%Re(1)": 1.0 / mu_g, - # Patch 1: Gas background - "patch_icpp(1)%geometry": 9, - "patch_icpp(1)%x_centroid": 0.0, - "patch_icpp(1)%y_centroid": 0.0, - "patch_icpp(1)%z_centroid": 0.0, - "patch_icpp(1)%length_x": x1 - x0, - "patch_icpp(1)%length_y": y1 - y0, - "patch_icpp(1)%length_z": z1 - z0, - "patch_icpp(1)%vel(1)": 0.0, - "patch_icpp(1)%vel(2)": 0.0, - "patch_icpp(1)%vel(3)": 0.0, - "patch_icpp(1)%pres": p_gas, - "patch_icpp(1)%alpha_rho(1)": eps * rho_l, - "patch_icpp(1)%alpha_rho(2)": (1.0 - eps) * rho_g, - "patch_icpp(1)%alpha(1)": eps, - "patch_icpp(1)%alpha(2)": 1.0 - eps, - "patch_icpp(1)%cf_val": 0, - # Patch 2: Left droplet - "patch_icpp(2)%geometry": 8, - "patch_icpp(2)%x_centroid": -sep, - "patch_icpp(2)%y_centroid": -b / 2.0, - "patch_icpp(2)%z_centroid": 0.0, - "patch_icpp(2)%radius": R, - "patch_icpp(2)%alter_patch(1)": "T", - "patch_icpp(2)%smoothen": "T", - "patch_icpp(2)%smooth_patch_id": 1, - "patch_icpp(2)%smooth_coeff": 0.95, - "patch_icpp(2)%vel(1)": Ud, - "patch_icpp(2)%vel(2)": 0.0, - "patch_icpp(2)%vel(3)": 0.0, - "patch_icpp(2)%pres": p_liq, - "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_l, - "patch_icpp(2)%alpha_rho(2)": eps * rho_g, - "patch_icpp(2)%alpha(1)": 1.0 - eps, - "patch_icpp(2)%alpha(2)": eps, - "patch_icpp(2)%cf_val": 1, - # Patch 3: Right droplet - "patch_icpp(3)%geometry": 8, - "patch_icpp(3)%x_centroid": sep, - "patch_icpp(3)%y_centroid": b / 2.0, - "patch_icpp(3)%z_centroid": 0.0, - "patch_icpp(3)%radius": R, - "patch_icpp(3)%alter_patch(1)": "T", - "patch_icpp(3)%smoothen": "T", - "patch_icpp(3)%smooth_patch_id": 1, - "patch_icpp(3)%smooth_coeff": 0.95, - "patch_icpp(3)%vel(1)": -Ud, - "patch_icpp(3)%vel(2)": 0.0, - "patch_icpp(3)%vel(3)": 0.0, - "patch_icpp(3)%pres": p_liq, - "patch_icpp(3)%alpha_rho(1)": (1.0 - eps) * rho_l, - "patch_icpp(3)%alpha_rho(2)": eps * rho_g, - "patch_icpp(3)%alpha(1)": 1.0 - eps, - "patch_icpp(3)%alpha(2)": eps, - "patch_icpp(3)%cf_val": 1, -} - -print(json.dumps(data)) From 92a7798d059f0c5e4450e0be97bfaebd299a3177 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Wed, 22 Apr 2026 23:20:03 -0400 Subject: [PATCH 06/10] Add droplet-coalescence analysis notebook Reads silo output via mfc.viz and reports cells per diameter and interface thickness (0.05<=alpha_1<=0.95) for a chosen step, with a |grad alpha| cross-check, a midplane sanity slice, and an interface-thickness time series. Works on both IC-check (t_step_stop=1) and full runs. Co-Authored-By: Claude Opus 4.7 --- .../3D_droplet_coalescence/analysis.ipynb | 318 ++++++++++++++++++ 1 file changed, 318 insertions(+) create mode 100644 examples/3D_droplet_coalescence/analysis.ipynb diff --git a/examples/3D_droplet_coalescence/analysis.ipynb b/examples/3D_droplet_coalescence/analysis.ipynb new file mode 100644 index 0000000000..d06a00c217 --- /dev/null +++ b/examples/3D_droplet_coalescence/analysis.ipynb @@ -0,0 +1,318 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Droplet-coalescence analysis\n", + "\n", + "Reads the silo output from this case directory and reports two resolution diagnostics:\n", + "\n", + "1. **Cells per diameter (CPD)** \u2014 grid resolution across the droplet, per axis.\n", + "2. **Interface thickness** \u2014 width of the diffuse interface ($0.05 \\le \\alpha_1 \\le 0.95$), in metres and in cells.\n", + "\n", + "Works on both an IC-check run (`t_step_stop=1, t_step_save=1`) and a full run. The time-series cell is a no-op when only one step exists.\n", + "\n", + "**Prerequisite:** `silo_hdf5/` must exist in this directory (run `./mfc.sh run case.py --case a` with `post_process` enabled, format=1).\n", + "\n", + "**Python env:** must have `h5py`, `numpy`, `matplotlib`. The MFC toolchain venv works:\n", + "\n", + "```\n", + "source ../../build/venv/bin/activate\n", + "jupyter lab analysis.ipynb\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "from pathlib import Path\n", + "\n", + "NOTEBOOK_DIR = Path.cwd().resolve()\n", + "MFC_ROOT = NOTEBOOK_DIR.parents[1]\n", + "sys.path.insert(0, str(MFC_ROOT / \"toolchain\"))\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from mfc.viz import assemble_silo, discover_timesteps\n", + "\n", + "D = 240e-6\n", + "R = D / 2.0\n", + "\n", + "# Mirror of the CASES table in case.py. Keep in sync if case.py changes.\n", + "CASES = {\n", + " \"a\": (0.2, 14.8, 0.20), \"b\": (0.5, 23.6, 0.10), \"c\": (8.6, 105.9, 0.08),\n", + " \"d\": (15.2, 139.8, 0.08), \"e\": (19.4, 158.0, 0.05), \"f\": (32.8, 210.8, 0.08),\n", + " \"g\": (37.2, 228.0, 0.01), \"h\": (61.4, 296.5, 0.06), \"i\": (61.3, 295.3, 0.11),\n", + " \"j\": (56.3, 288.9, 0.13), \"k\": (70.8, 327.7, 0.25), \"l\": (48.1, 270.1, 0.39),\n", + " \"m\": (60.1, 302.8, 0.55), \"n\": (65.1, 320.3, 0.49), \"o\": (60.8, 313.7, 0.68),\n", + " \"p\": (64.9, 312.8, 0.71), \"q\": (48.8, 260.3, 0.72), \"r\": (14.5, 149.1, 0.34),\n", + "}\n", + "\n", + "CASE_LABEL = \"a\"\n", + "_, _, B = CASES[CASE_LABEL]\n", + "sep = 0.55 * D\n", + "b_impact = B * D\n", + "\n", + "DROPLETS = {\n", + " \"left\": (-sep, -b_impact / 2.0, 0.0),\n", + " \"right\": (+sep, +b_impact / 2.0, 0.0),\n", + "}\n", + "\n", + "CASE_DIR = str(NOTEBOOK_DIR)\n", + "STEP = None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "steps = discover_timesteps(CASE_DIR, \"silo\")\n", + "if not steps:\n", + " raise RuntimeError(\n", + " f\"No silo timesteps found under {CASE_DIR!r}. \"\n", + " \"Run pre_process + post_process (format=1) first.\"\n", + " )\n", + "\n", + "step = STEP if STEP is not None else steps[0]\n", + "if step not in steps:\n", + " raise ValueError(f\"step={step} not in available steps {steps}\")\n", + "\n", + "print(f\"Available steps ({len(steps)}): {steps[:10]}{' ...' if len(steps) > 10 else ''}\")\n", + "print(f\"Using step: {step}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data = assemble_silo(CASE_DIR, step)\n", + "print(f\"ndim = {data.ndim}\")\n", + "print(f\"grid = {len(data.x_cc)} x {len(data.y_cc)} x {len(data.z_cc)}\")\n", + "print(f\"variables = {sorted(data.variables)}\")\n", + "\n", + "ALPHA_NAME = next((n for n in (\"alpha1\", \"alpha_1\") if n in data.variables), None)\n", + "if ALPHA_NAME is None:\n", + " raise KeyError(f\"No alpha1-like variable in {sorted(data.variables)}\")\n", + "alpha = data.variables[ALPHA_NAME]\n", + "print(f\"alpha: {ALPHA_NAME} shape={alpha.shape} range=[{alpha.min():.3g}, {alpha.max():.3g}]\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x, y, z = data.x_cc, data.y_cc, data.z_cc\n", + "\n", + "def _axis_stats(coord):\n", + " d = np.diff(coord)\n", + " return d.mean(), d.min(), d.max()\n", + "\n", + "stats = {\n", + " \"x\": (len(x), *_axis_stats(x)),\n", + " \"y\": (len(y), *_axis_stats(y)),\n", + " \"z\": (len(z), *_axis_stats(z)),\n", + "}\n", + "\n", + "print(f\"{'axis':<5}{'N':>6}{'\u0394 [um]':>12}{'CPD = D/\u0394':>14}{'nonuniformity':>18}\")\n", + "print(\"-\" * 55)\n", + "for name, (n, d, dmin, dmax) in stats.items():\n", + " u = (dmax - dmin) / max(dmin, 1e-300)\n", + " tag = \"uniform\" if u < 1e-6 else f\"{u:.1e}\"\n", + " print(f\"{name:<5}{n:>6d}{d*1e6:>12.4f}{D/d:>14.2f}{tag:>18}\")\n", + "\n", + "dx = stats[\"x\"][1]\n", + "dy = stats[\"y\"][1]\n", + "dz = stats[\"z\"][1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def interface_thickness_x_ray(alpha, x_cc, y_cc, z_cc, xc, yc, zc, R):\n", + " jy = int(np.argmin(np.abs(y_cc - yc)))\n", + " kz = int(np.argmin(np.abs(z_cc - zc)))\n", + "\n", + " if xc >= 0:\n", + " mask = (x_cc >= xc) & (x_cc <= xc + 2.0 * R)\n", + " else:\n", + " mask = (x_cc <= xc) & (x_cc >= xc - 2.0 * R)\n", + "\n", + " x_seg = x_cc[mask]\n", + " a_seg = alpha[mask, jy, kz]\n", + "\n", + " if a_seg.min() > 0.05 or a_seg.max() < 0.95:\n", + " raise RuntimeError(\n", + " f\"alpha range on ray = [{a_seg.min():.3g}, {a_seg.max():.3g}]; \"\n", + " \"interface not fully captured on this ray\"\n", + " )\n", + "\n", + " order = np.argsort(a_seg)\n", + " a_sorted, x_sorted = a_seg[order], x_seg[order]\n", + " x_05 = float(np.interp(0.05, a_sorted, x_sorted))\n", + " x_95 = float(np.interp(0.95, a_sorted, x_sorted))\n", + "\n", + " t_m = abs(x_95 - x_05)\n", + " dx_local = float(np.mean(np.diff(x_seg)))\n", + " return t_m, t_m / dx_local, x_seg, a_seg\n", + "\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 4))\n", + "results = {}\n", + "for label, (xc, yc, zc) in DROPLETS.items():\n", + " t_m, t_c, x_seg, a_seg = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", + " results[label] = (t_m, t_c)\n", + " ax.plot((x_seg - xc) / D, a_seg, \"o-\", ms=3, label=f\"{label}: {t_c:.2f} cells ({t_m*1e6:.2f} um)\")\n", + "\n", + "ax.axhline(0.05, color=\"k\", lw=0.5, ls=\":\")\n", + "ax.axhline(0.95, color=\"k\", lw=0.5, ls=\":\")\n", + "ax.set_xlabel(r\"$(x - x_c) / D$\")\n", + "ax.set_ylabel(r\"$\\alpha_1$\")\n", + "ax.set_title(f\"Interface profile along x (step {step})\")\n", + "ax.legend()\n", + "ax.grid(alpha=0.3)\n", + "plt.show()\n", + "\n", + "t_mean_cells = float(np.mean([v[1] for v in results.values()]))\n", + "print(f\"Interface thickness (mean over droplets): {t_mean_cells:.2f} cells\")\n", + "for k, (t_m, t_c) in results.items():\n", + " print(f\" {k:>5}: {t_c:.2f} cells = {t_m*1e6:.2f} um\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for label, (xc, yc, zc) in DROPLETS.items():\n", + " mx = (x >= xc - 1.5 * R) & (x <= xc + 1.5 * R)\n", + " my = (y >= yc - 1.5 * R) & (y <= yc + 1.5 * R)\n", + " mz = (z >= zc - 1.5 * R) & (z <= zc + 1.5 * R)\n", + " a_box = alpha[np.ix_(mx, my, mz)]\n", + " gx_box = np.gradient(a_box, x[mx], axis=0)\n", + " gy_box = np.gradient(a_box, y[my], axis=1)\n", + " gz_box = np.gradient(a_box, z[mz], axis=2)\n", + " gmag = np.sqrt(gx_box**2 + gy_box**2 + gz_box**2)\n", + " t_grad_m = 0.9 / gmag.max()\n", + " dx_local = float(np.mean(np.diff(x[mx])))\n", + " print(f\"{label:>5}: max|\u2207\u03b1|={gmag.max():.3g} 1/m \u2192 thickness \u2248 {t_grad_m*1e6:.2f} um ({t_grad_m/dx_local:.2f} cells)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "kz0 = int(np.argmin(np.abs(z - 0.0)))\n", + "slice_xy = alpha[:, :, kz0]\n", + "\n", + "fig, ax = plt.subplots(figsize=(9, 6))\n", + "im = ax.imshow(\n", + " slice_xy.T,\n", + " origin=\"lower\",\n", + " extent=[x[0] / D, x[-1] / D, y[0] / D, y[-1] / D],\n", + " cmap=\"viridis\",\n", + " aspect=\"equal\",\n", + " vmin=0.0,\n", + " vmax=1.0,\n", + ")\n", + "ax.contour(\n", + " x / D, y / D, slice_xy.T,\n", + " levels=[0.05, 0.5, 0.95],\n", + " colors=[\"w\", \"r\", \"w\"],\n", + " linewidths=[0.6, 1.0, 0.6],\n", + ")\n", + "for label, (xc, yc, _zc) in DROPLETS.items():\n", + " ax.plot(xc / D, yc / D, \"r+\", ms=12, mew=2)\n", + "ax.set_xlabel(\"x / D\")\n", + "ax.set_ylabel(\"y / D\")\n", + "ax.set_title(f\"alpha_1 on z=0 midplane (step {step})\")\n", + "fig.colorbar(im, ax=ax, label=r\"$\\alpha_1$\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if len(steps) < 2:\n", + " print(\"Only one step available \u2014 skipping time-series.\")\n", + "else:\n", + " series = {k: [] for k in DROPLETS}\n", + " for s in steps:\n", + " d_s = assemble_silo(CASE_DIR, s)\n", + " a_s = d_s.variables[ALPHA_NAME]\n", + " xs, ys, zs = d_s.x_cc, d_s.y_cc, d_s.z_cc\n", + " for label, (xc0, _yc0, _zc0) in DROPLETS.items():\n", + " mx = (xs >= xc0 - 1.5 * R) & (xs <= xc0 + 1.5 * R)\n", + " w = a_s[mx, :, :]\n", + " mass = float(w.sum())\n", + " if mass <= 0:\n", + " series[label].append(np.nan)\n", + " continue\n", + " xw = float(np.sum(w.sum(axis=(1, 2)) * xs[mx]) / mass)\n", + " yw = float(np.sum(w.sum(axis=(0, 2)) * ys) / mass)\n", + " zw = float(np.sum(w.sum(axis=(0, 1)) * zs) / mass)\n", + " try:\n", + " _, t_c, _, _ = interface_thickness_x_ray(a_s, xs, ys, zs, xw, yw, zw, R)\n", + " except RuntimeError:\n", + " t_c = np.nan\n", + " series[label].append(t_c)\n", + "\n", + " fig, ax = plt.subplots(figsize=(7, 4))\n", + " for label, vals in series.items():\n", + " ax.plot(steps, vals, \"o-\", label=label)\n", + " ax.set_xlabel(\"step\")\n", + " ax.set_ylabel(\"interface thickness [cells]\")\n", + " ax.set_title(\"Interface thickness over time\")\n", + " ax.legend()\n", + " ax.grid(alpha=0.3)\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"case: {CASE_LABEL} step: {step}\")\n", + "print(f\"grid: {stats['x'][0]} x {stats['y'][0]} x {stats['z'][0]}\")\n", + "print(f\"CPD : x={D/dx:.1f} y={D/dy:.1f} z={D/dz:.1f}\")\n", + "for k, (t_m, t_c) in results.items():\n", + " print(f\"interface thickness ({k:>5}): {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n", + "print(f\"mean interface thickness: {t_mean_cells:.2f} cells\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 6fc119421f8127fcc02f9609f6423c0c776651f7 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Thu, 23 Apr 2026 16:44:37 -0400 Subject: [PATCH 07/10] Rework droplet-coalescence analysis notebook to show the mesh first Lead with cell-resolved views before any derived scalars: (1) midplane pcolormesh with mesh lines, (2) close-up on the right-droplet outer face with per-cell alpha values annotated in the 0.05-0.95 interface band, (3) x-ray profile as step-style per-cell values with cell boundaries drawn. Thickness extraction + |grad alpha| cross-check moved after the plots. Co-Authored-By: Claude Opus 4.7 --- .../3D_droplet_coalescence/analysis.ipynb | 311 ++++++++++++------ 1 file changed, 217 insertions(+), 94 deletions(-) diff --git a/examples/3D_droplet_coalescence/analysis.ipynb b/examples/3D_droplet_coalescence/analysis.ipynb index d06a00c217..a113590ad9 100644 --- a/examples/3D_droplet_coalescence/analysis.ipynb +++ b/examples/3D_droplet_coalescence/analysis.ipynb @@ -6,12 +6,19 @@ "source": [ "# Droplet-coalescence analysis\n", "\n", - "Reads the silo output from this case directory and reports two resolution diagnostics:\n", + "Reads the silo output from this case directory and visualizes the mesh, the interface, and the actual per-cell $\\alpha_1$ values before reporting any derived scalars.\n", "\n", - "1. **Cells per diameter (CPD)** \u2014 grid resolution across the droplet, per axis.\n", - "2. **Interface thickness** \u2014 width of the diffuse interface ($0.05 \\le \\alpha_1 \\le 0.95$), in metres and in cells.\n", + "Order of cells:\n", "\n", - "Works on both an IC-check run (`t_step_stop=1, t_step_save=1`) and a full run. The time-series cell is a no-op when only one step exists.\n", + "1. Config, helpers\n", + "2. Discover timesteps, load the target step\n", + "3. Cells per diameter on each axis\n", + "4. Midplane $\\alpha_1$ field, cell-by-cell (`pcolormesh` with mesh lines)\n", + "5. Close-up at the outer face of the right droplet, with per-cell $\\alpha$ values annotated where $0.05 < \\alpha < 0.95$\n", + "6. $\\alpha_1$ sampled along an x-ray through each droplet centre, plotted as discrete per-cell values with cell boundaries marked\n", + "7. Interface thickness (from the same ray) and a $|\\nabla \\alpha|$ cross-check, printed after the plots\n", + "8. Time series (skipped when only one step exists)\n", + "9. Summary\n", "\n", "**Prerequisite:** `silo_hdf5/` must exist in this directory (run `./mfc.sh run case.py --case a` with `post_process` enabled, format=1).\n", "\n", @@ -64,7 +71,32 @@ "}\n", "\n", "CASE_DIR = str(NOTEBOOK_DIR)\n", - "STEP = None" + "STEP = None\n", + "\n", + "\n", + "def interface_thickness_x_ray(alpha, x_cc, y_cc, z_cc, xc, yc, zc, R):\n", + " jy = int(np.argmin(np.abs(y_cc - yc)))\n", + " kz = int(np.argmin(np.abs(z_cc - zc)))\n", + " if xc >= 0:\n", + " mask = (x_cc >= xc) & (x_cc <= xc + 2.0 * R)\n", + " else:\n", + " mask = (x_cc <= xc) & (x_cc >= xc - 2.0 * R)\n", + "\n", + " x_seg = x_cc[mask]\n", + " a_seg = alpha[mask, jy, kz]\n", + "\n", + " if a_seg.min() > 0.05 or a_seg.max() < 0.95:\n", + " raise RuntimeError(\n", + " f\"alpha range on ray = [{a_seg.min():.3g}, {a_seg.max():.3g}]; \"\n", + " \"interface not fully captured on this ray\"\n", + " )\n", + "\n", + " order = np.argsort(a_seg)\n", + " x_05 = float(np.interp(0.05, a_seg[order], x_seg[order]))\n", + " x_95 = float(np.interp(0.95, a_seg[order], x_seg[order]))\n", + " t_m = abs(x_95 - x_05)\n", + " dx_local = float(np.mean(np.diff(x_seg)))\n", + " return t_m, t_m / dx_local, x_seg, a_seg, x_05, x_95" ] }, { @@ -85,15 +117,8 @@ " raise ValueError(f\"step={step} not in available steps {steps}\")\n", "\n", "print(f\"Available steps ({len(steps)}): {steps[:10]}{' ...' if len(steps) > 10 else ''}\")\n", - "print(f\"Using step: {step}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + "print(f\"Using step: {step}\")\n", + "\n", "data = assemble_silo(CASE_DIR, step)\n", "print(f\"ndim = {data.ndim}\")\n", "print(f\"grid = {len(data.x_cc)} x {len(data.y_cc)} x {len(data.z_cc)}\")\n", @@ -103,6 +128,7 @@ "if ALPHA_NAME is None:\n", " raise KeyError(f\"No alpha1-like variable in {sorted(data.variables)}\")\n", "alpha = data.variables[ALPHA_NAME]\n", + "x, y, z = data.x_cc, data.y_cc, data.z_cc\n", "print(f\"alpha: {ALPHA_NAME} shape={alpha.shape} range=[{alpha.min():.3g}, {alpha.max():.3g}]\")" ] }, @@ -112,8 +138,6 @@ "metadata": {}, "outputs": [], "source": [ - "x, y, z = data.x_cc, data.y_cc, data.z_cc\n", - "\n", "def _axis_stats(coord):\n", " d = np.diff(coord)\n", " return d.mean(), d.min(), d.max()\n", @@ -136,60 +160,172 @@ "dz = stats[\"z\"][1]" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Midplane $\\alpha_1$ field on $z=0$\n", + "\n", + "`pcolormesh` with `shading=\"nearest\"` renders each grid cell as a discrete coloured quad centred on its cell-centre coordinate. At full 400x280 resolution the mesh lines look like noise; zoom in (toolbar or `ax.set_xlim / set_ylim`) to see them resolve into cells." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "def interface_thickness_x_ray(alpha, x_cc, y_cc, z_cc, xc, yc, zc, R):\n", - " jy = int(np.argmin(np.abs(y_cc - yc)))\n", - " kz = int(np.argmin(np.abs(z_cc - zc)))\n", - "\n", - " if xc >= 0:\n", - " mask = (x_cc >= xc) & (x_cc <= xc + 2.0 * R)\n", - " else:\n", - " mask = (x_cc <= xc) & (x_cc >= xc - 2.0 * R)\n", + "kz0 = int(np.argmin(np.abs(z - 0.0)))\n", + "slice_xy = alpha[:, :, kz0]\n", "\n", - " x_seg = x_cc[mask]\n", - " a_seg = alpha[mask, jy, kz]\n", + "fig, ax = plt.subplots(figsize=(10, 6))\n", + "pc = ax.pcolormesh(\n", + " x / D, y / D, slice_xy.T,\n", + " shading=\"nearest\",\n", + " cmap=\"viridis\",\n", + " vmin=0.0, vmax=1.0,\n", + " edgecolors=\"k\",\n", + " linewidths=0.05,\n", + ")\n", + "ax.contour(x / D, y / D, slice_xy.T, levels=[0.5], colors=\"r\", linewidths=0.8)\n", + "for label, (xc, yc, _zc) in DROPLETS.items():\n", + " ax.plot(xc / D, yc / D, \"r+\", ms=12, mew=2)\n", + "ax.set_xlabel(\"x / D\")\n", + "ax.set_ylabel(\"y / D\")\n", + "ax.set_aspect(\"equal\")\n", + "ax.set_title(f\"$\\\\alpha_1$ on z=0 midplane (step {step})\")\n", + "fig.colorbar(pc, ax=ax, label=r\"$\\alpha_1$\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interface close-up\n", "\n", - " if a_seg.min() > 0.05 or a_seg.max() < 0.95:\n", - " raise RuntimeError(\n", - " f\"alpha range on ray = [{a_seg.min():.3g}, {a_seg.max():.3g}]; \"\n", - " \"interface not fully captured on this ray\"\n", - " )\n", + "Zoom onto the outer face of the **right** droplet on the $z=0$ midplane. Cell edges are drawn, and the per-cell $\\alpha_1$ value is printed inside every cell where $0.05 < \\alpha_1 < 0.95$ (the interface band). Red line is the $\\alpha_1 = 0.5$ isocontour." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "xc_r, yc_r, _zc_r = DROPLETS[\"right\"]\n", + "x_lo, x_hi = xc_r + 0.2 * R, xc_r + 1.2 * R\n", + "y_lo, y_hi = yc_r - 0.5 * R, yc_r + 0.5 * R\n", "\n", - " order = np.argsort(a_seg)\n", - " a_sorted, x_sorted = a_seg[order], x_seg[order]\n", - " x_05 = float(np.interp(0.05, a_sorted, x_sorted))\n", - " x_95 = float(np.interp(0.95, a_sorted, x_sorted))\n", + "mx = (x >= x_lo) & (x <= x_hi)\n", + "my = (y >= y_lo) & (y <= y_hi)\n", + "sub_x = x[mx]\n", + "sub_y = y[my]\n", + "sub_alpha = slice_xy[np.ix_(mx, my)]\n", "\n", - " t_m = abs(x_95 - x_05)\n", - " dx_local = float(np.mean(np.diff(x_seg)))\n", - " return t_m, t_m / dx_local, x_seg, a_seg\n", + "fig, ax = plt.subplots(figsize=(10, 7))\n", + "pc = ax.pcolormesh(\n", + " sub_x / D, sub_y / D, sub_alpha.T,\n", + " shading=\"nearest\",\n", + " cmap=\"viridis\",\n", + " vmin=0.0, vmax=1.0,\n", + " edgecolors=\"k\",\n", + " linewidths=0.4,\n", + ")\n", + "ax.contour(sub_x / D, sub_y / D, sub_alpha.T, levels=[0.5], colors=\"r\", linewidths=1.5)\n", "\n", + "for i, xv in enumerate(sub_x):\n", + " for j, yv in enumerate(sub_y):\n", + " a_val = float(sub_alpha[i, j])\n", + " if 0.05 < a_val < 0.95:\n", + " txt_color = \"w\" if a_val < 0.5 else \"k\"\n", + " ax.text(xv / D, yv / D, f\"{a_val:.2f}\",\n", + " ha=\"center\", va=\"center\", fontsize=6, color=txt_color)\n", "\n", - "fig, ax = plt.subplots(figsize=(7, 4))\n", - "results = {}\n", - "for label, (xc, yc, zc) in DROPLETS.items():\n", - " t_m, t_c, x_seg, a_seg = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", - " results[label] = (t_m, t_c)\n", - " ax.plot((x_seg - xc) / D, a_seg, \"o-\", ms=3, label=f\"{label}: {t_c:.2f} cells ({t_m*1e6:.2f} um)\")\n", - "\n", - "ax.axhline(0.05, color=\"k\", lw=0.5, ls=\":\")\n", - "ax.axhline(0.95, color=\"k\", lw=0.5, ls=\":\")\n", - "ax.set_xlabel(r\"$(x - x_c) / D$\")\n", - "ax.set_ylabel(r\"$\\alpha_1$\")\n", - "ax.set_title(f\"Interface profile along x (step {step})\")\n", - "ax.legend()\n", - "ax.grid(alpha=0.3)\n", + "ax.set_xlabel(\"x / D\")\n", + "ax.set_ylabel(\"y / D\")\n", + "ax.set_aspect(\"equal\")\n", + "ax.set_title(f\"Right-droplet outer face, z=0 (step {step})\")\n", + "fig.colorbar(pc, ax=ax, label=r\"$\\alpha_1$\")\n", "plt.show()\n", "\n", - "t_mean_cells = float(np.mean([v[1] for v in results.values()]))\n", - "print(f\"Interface thickness (mean over droplets): {t_mean_cells:.2f} cells\")\n", - "for k, (t_m, t_c) in results.items():\n", - " print(f\" {k:>5}: {t_c:.2f} cells = {t_m*1e6:.2f} um\")" + "n_interface = int(((sub_alpha > 0.05) & (sub_alpha < 0.95)).sum())\n", + "print(f\"Interface cells in this close-up (0.05 < alpha < 0.95): {n_interface}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $\\alpha_1$ along an x-ray through each droplet centre\n", + "\n", + "One data point per cell. Step line uses `drawstyle=\"steps-mid\"` so the horizontal segment at each cell centre has length $\\Delta x$. Faint vertical lines are cell boundaries. Dotted horizontals at 0.05 and 0.95 show what the thickness extraction will interpolate between." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\n", + "for ax, (label, (xc, yc, zc)) in zip(axes, DROPLETS.items()):\n", + " try:\n", + " _, _, x_seg, a_seg, x_05, x_95 = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", + " have_crossings = True\n", + " except RuntimeError as err:\n", + " have_crossings = False\n", + " msg = str(err)\n", + " jy = int(np.argmin(np.abs(y - yc)))\n", + " kz = int(np.argmin(np.abs(z - zc)))\n", + " if xc >= 0:\n", + " mask = (x >= xc) & (x <= xc + 2.0 * R)\n", + " else:\n", + " mask = (x <= xc) & (x >= xc - 2.0 * R)\n", + " x_seg = x[mask]\n", + " a_seg = alpha[mask, jy, kz]\n", + "\n", + " xrel = (x_seg - xc) / D\n", + " ax.plot(xrel, a_seg, drawstyle=\"steps-mid\", lw=1.2, color=\"C0\")\n", + " ax.plot(xrel, a_seg, \"o\", color=\"C0\", ms=4)\n", + "\n", + " edges = np.concatenate((\n", + " [x_seg[0] - 0.5 * (x_seg[1] - x_seg[0])],\n", + " 0.5 * (x_seg[:-1] + x_seg[1:]),\n", + " [x_seg[-1] + 0.5 * (x_seg[-1] - x_seg[-2])],\n", + " ))\n", + " for e in edges:\n", + " ax.axvline((e - xc) / D, color=\"k\", lw=0.2, alpha=0.3)\n", + "\n", + " ax.axhline(0.05, color=\"k\", lw=0.5, ls=\":\")\n", + " ax.axhline(0.95, color=\"k\", lw=0.5, ls=\":\")\n", + "\n", + " if have_crossings:\n", + " ax.plot([(x_05 - xc) / D, (x_95 - xc) / D], [0.05, 0.95], \"x\", color=\"r\", ms=10, mew=2,\n", + " label=\"interp crossings\")\n", + " ax.legend(loc=\"best\", fontsize=8)\n", + "\n", + " ax.set_xlabel(r\"$(x - x_c) / D$\")\n", + " ax.set_title(label if have_crossings else f\"{label} \u2014 no 0.05/0.95 crossing on this ray\")\n", + " ax.grid(alpha=0.2)\n", + "\n", + "axes[0].set_ylabel(r\"$\\alpha_1$\")\n", + "fig.suptitle(f\"$\\\\alpha_1$ along x-ray through droplet centre (step {step})\")\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Derived: interface thickness\n", + "\n", + "Two independent estimates of interface thickness, derived from the field you just looked at:\n", + "- **Ray method** \u2014 $|x(\\alpha=0.95) - x(\\alpha=0.05)|$ interpolated along the x-ray through each droplet centre.\n", + "- **$|\\nabla \\alpha|$ method** \u2014 $0.9 / \\max |\\nabla \\alpha|$ over a 1.5R cube around each droplet. Assumes an approximately linear transition.\n", + "\n", + "If these two disagree by more than a factor of ~2, the interface is not well-represented by an axis-aligned ray (e.g. tilted interface, or the ray misses the interface)." ] }, { @@ -198,6 +334,25 @@ "metadata": {}, "outputs": [], "source": [ + "results = {}\n", + "for label, (xc, yc, zc) in DROPLETS.items():\n", + " try:\n", + " t_m, t_c, *_ = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", + " results[label] = (t_m, t_c)\n", + " except RuntimeError as err:\n", + " print(f\" {label:>5}: ray method failed \u2014 {err}\")\n", + " results[label] = (np.nan, np.nan)\n", + "\n", + "print(\"Ray method (0.05 <= alpha_1 <= 0.95):\")\n", + "for k, (t_m, t_c) in results.items():\n", + " if not np.isnan(t_c):\n", + " print(f\" {k:>5}: {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n", + "vals = [v[1] for v in results.values() if not np.isnan(v[1])]\n", + "t_mean_cells = float(np.mean(vals)) if vals else float(\"nan\")\n", + "print(f\" mean : {t_mean_cells:.2f} cells\")\n", + "\n", + "print()\n", + "print(\"|grad alpha| method (max over a 1.5R cube):\")\n", "for label, (xc, yc, zc) in DROPLETS.items():\n", " mx = (x >= xc - 1.5 * R) & (x <= xc + 1.5 * R)\n", " my = (y >= yc - 1.5 * R) & (y <= yc + 1.5 * R)\n", @@ -209,41 +364,7 @@ " gmag = np.sqrt(gx_box**2 + gy_box**2 + gz_box**2)\n", " t_grad_m = 0.9 / gmag.max()\n", " dx_local = float(np.mean(np.diff(x[mx])))\n", - " print(f\"{label:>5}: max|\u2207\u03b1|={gmag.max():.3g} 1/m \u2192 thickness \u2248 {t_grad_m*1e6:.2f} um ({t_grad_m/dx_local:.2f} cells)\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "kz0 = int(np.argmin(np.abs(z - 0.0)))\n", - "slice_xy = alpha[:, :, kz0]\n", - "\n", - "fig, ax = plt.subplots(figsize=(9, 6))\n", - "im = ax.imshow(\n", - " slice_xy.T,\n", - " origin=\"lower\",\n", - " extent=[x[0] / D, x[-1] / D, y[0] / D, y[-1] / D],\n", - " cmap=\"viridis\",\n", - " aspect=\"equal\",\n", - " vmin=0.0,\n", - " vmax=1.0,\n", - ")\n", - "ax.contour(\n", - " x / D, y / D, slice_xy.T,\n", - " levels=[0.05, 0.5, 0.95],\n", - " colors=[\"w\", \"r\", \"w\"],\n", - " linewidths=[0.6, 1.0, 0.6],\n", - ")\n", - "for label, (xc, yc, _zc) in DROPLETS.items():\n", - " ax.plot(xc / D, yc / D, \"r+\", ms=12, mew=2)\n", - "ax.set_xlabel(\"x / D\")\n", - "ax.set_ylabel(\"y / D\")\n", - "ax.set_title(f\"alpha_1 on z=0 midplane (step {step})\")\n", - "fig.colorbar(im, ax=ax, label=r\"$\\alpha_1$\")\n", - "plt.show()" + " print(f\" {label:>5}: max|grad alpha|={gmag.max():.3g} 1/m -> {t_grad_m*1e6:.2f} um ({t_grad_m/dx_local:.2f} cells)\")" ] }, { @@ -271,7 +392,7 @@ " yw = float(np.sum(w.sum(axis=(0, 2)) * ys) / mass)\n", " zw = float(np.sum(w.sum(axis=(0, 1)) * zs) / mass)\n", " try:\n", - " _, t_c, _, _ = interface_thickness_x_ray(a_s, xs, ys, zs, xw, yw, zw, R)\n", + " _, t_c, *_ = interface_thickness_x_ray(a_s, xs, ys, zs, xw, yw, zw, R)\n", " except RuntimeError:\n", " t_c = np.nan\n", " series[label].append(t_c)\n", @@ -297,8 +418,10 @@ "print(f\"grid: {stats['x'][0]} x {stats['y'][0]} x {stats['z'][0]}\")\n", "print(f\"CPD : x={D/dx:.1f} y={D/dy:.1f} z={D/dz:.1f}\")\n", "for k, (t_m, t_c) in results.items():\n", - " print(f\"interface thickness ({k:>5}): {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n", - "print(f\"mean interface thickness: {t_mean_cells:.2f} cells\")" + " if not np.isnan(t_c):\n", + " print(f\"interface thickness ({k:>5}): {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n", + "if not np.isnan(t_mean_cells):\n", + " print(f\"mean interface thickness: {t_mean_cells:.2f} cells\")" ] } ], From af9f42b34327b596a73fa56a6467fe7285736996 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Thu, 23 Apr 2026 18:09:47 -0400 Subject: [PATCH 08/10] Harden analysis notebook against empty slices and missing data - interface_thickness_x_ray: raise RuntimeError (not bare ValueError) when the x-ray mask covers no cells - x-ray plot, derived-thickness, and time-series cells catch (RuntimeError, ValueError) and guard against empty slices / zero gradients - Summary cell uses locals().get() so it renders even if earlier cells raised Co-Authored-By: Claude Opus 4.7 --- .../3D_droplet_coalescence/analysis.ipynb | 240 +++--------------- 1 file changed, 37 insertions(+), 203 deletions(-) diff --git a/examples/3D_droplet_coalescence/analysis.ipynb b/examples/3D_droplet_coalescence/analysis.ipynb index a113590ad9..c629b54f5f 100644 --- a/examples/3D_droplet_coalescence/analysis.ipynb +++ b/examples/3D_droplet_coalescence/analysis.ipynb @@ -2,6 +2,7 @@ "cells": [ { "cell_type": "markdown", + "id": "7fb27b941602401d91542211134fc71a", "metadata": {}, "source": [ "# Droplet-coalescence analysis\n", @@ -33,84 +34,21 @@ { "cell_type": "code", "execution_count": null, + "id": "acae54e37e7d407bbb7b55eff062a284", "metadata": {}, "outputs": [], - "source": [ - "import sys\n", - "from pathlib import Path\n", - "\n", - "NOTEBOOK_DIR = Path.cwd().resolve()\n", - "MFC_ROOT = NOTEBOOK_DIR.parents[1]\n", - "sys.path.insert(0, str(MFC_ROOT / \"toolchain\"))\n", - "\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from mfc.viz import assemble_silo, discover_timesteps\n", - "\n", - "D = 240e-6\n", - "R = D / 2.0\n", - "\n", - "# Mirror of the CASES table in case.py. Keep in sync if case.py changes.\n", - "CASES = {\n", - " \"a\": (0.2, 14.8, 0.20), \"b\": (0.5, 23.6, 0.10), \"c\": (8.6, 105.9, 0.08),\n", - " \"d\": (15.2, 139.8, 0.08), \"e\": (19.4, 158.0, 0.05), \"f\": (32.8, 210.8, 0.08),\n", - " \"g\": (37.2, 228.0, 0.01), \"h\": (61.4, 296.5, 0.06), \"i\": (61.3, 295.3, 0.11),\n", - " \"j\": (56.3, 288.9, 0.13), \"k\": (70.8, 327.7, 0.25), \"l\": (48.1, 270.1, 0.39),\n", - " \"m\": (60.1, 302.8, 0.55), \"n\": (65.1, 320.3, 0.49), \"o\": (60.8, 313.7, 0.68),\n", - " \"p\": (64.9, 312.8, 0.71), \"q\": (48.8, 260.3, 0.72), \"r\": (14.5, 149.1, 0.34),\n", - "}\n", - "\n", - "CASE_LABEL = \"a\"\n", - "_, _, B = CASES[CASE_LABEL]\n", - "sep = 0.55 * D\n", - "b_impact = B * D\n", - "\n", - "DROPLETS = {\n", - " \"left\": (-sep, -b_impact / 2.0, 0.0),\n", - " \"right\": (+sep, +b_impact / 2.0, 0.0),\n", - "}\n", - "\n", - "CASE_DIR = str(NOTEBOOK_DIR)\n", - "STEP = None\n", - "\n", - "\n", - "def interface_thickness_x_ray(alpha, x_cc, y_cc, z_cc, xc, yc, zc, R):\n", - " jy = int(np.argmin(np.abs(y_cc - yc)))\n", - " kz = int(np.argmin(np.abs(z_cc - zc)))\n", - " if xc >= 0:\n", - " mask = (x_cc >= xc) & (x_cc <= xc + 2.0 * R)\n", - " else:\n", - " mask = (x_cc <= xc) & (x_cc >= xc - 2.0 * R)\n", - "\n", - " x_seg = x_cc[mask]\n", - " a_seg = alpha[mask, jy, kz]\n", - "\n", - " if a_seg.min() > 0.05 or a_seg.max() < 0.95:\n", - " raise RuntimeError(\n", - " f\"alpha range on ray = [{a_seg.min():.3g}, {a_seg.max():.3g}]; \"\n", - " \"interface not fully captured on this ray\"\n", - " )\n", - "\n", - " order = np.argsort(a_seg)\n", - " x_05 = float(np.interp(0.05, a_seg[order], x_seg[order]))\n", - " x_95 = float(np.interp(0.95, a_seg[order], x_seg[order]))\n", - " t_m = abs(x_95 - x_05)\n", - " dx_local = float(np.mean(np.diff(x_seg)))\n", - " return t_m, t_m / dx_local, x_seg, a_seg, x_05, x_95" - ] + "source": "import sys\nfrom pathlib import Path\n\nNOTEBOOK_DIR = Path.cwd().resolve()\nMFC_ROOT = NOTEBOOK_DIR.parents[1]\nsys.path.insert(0, str(MFC_ROOT / \"toolchain\"))\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom mfc.viz import assemble_silo, discover_timesteps\n\nD = 240e-6\nR = D / 2.0\n\n# Mirror of the CASES table in case.py. Keep in sync if case.py changes.\nCASES = {\n \"a\": (0.2, 14.8, 0.20), \"b\": (0.5, 23.6, 0.10), \"c\": (8.6, 105.9, 0.08),\n \"d\": (15.2, 139.8, 0.08), \"e\": (19.4, 158.0, 0.05), \"f\": (32.8, 210.8, 0.08),\n \"g\": (37.2, 228.0, 0.01), \"h\": (61.4, 296.5, 0.06), \"i\": (61.3, 295.3, 0.11),\n \"j\": (56.3, 288.9, 0.13), \"k\": (70.8, 327.7, 0.25), \"l\": (48.1, 270.1, 0.39),\n \"m\": (60.1, 302.8, 0.55), \"n\": (65.1, 320.3, 0.49), \"o\": (60.8, 313.7, 0.68),\n \"p\": (64.9, 312.8, 0.71), \"q\": (48.8, 260.3, 0.72), \"r\": (14.5, 149.1, 0.34),\n}\n\nCASE_LABEL = \"a\"\n_, _, B = CASES[CASE_LABEL]\nsep = 0.55 * D\nb_impact = B * D\n\nDROPLETS = {\n \"left\": (-sep, -b_impact / 2.0, 0.0),\n \"right\": (+sep, +b_impact / 2.0, 0.0),\n}\n\nCASE_DIR = str(NOTEBOOK_DIR)\nSTEP = None\n\n\ndef interface_thickness_x_ray(alpha, x_cc, y_cc, z_cc, xc, yc, zc, R):\n jy = int(np.argmin(np.abs(y_cc - yc)))\n kz = int(np.argmin(np.abs(z_cc - zc)))\n if xc >= 0:\n mask = (x_cc >= xc) & (x_cc <= xc + 2.0 * R)\n else:\n mask = (x_cc <= xc) & (x_cc >= xc - 2.0 * R)\n\n x_seg = x_cc[mask]\n a_seg = alpha[mask, jy, kz]\n\n if a_seg.size == 0:\n raise RuntimeError(\n f\"x-ray at y={yc:.3g}, z={zc:.3g} through xc={xc:.3g} intersects no grid cells; \"\n \"check CASE_LABEL matches the run\"\n )\n\n if a_seg.min() > 0.05 or a_seg.max() < 0.95:\n raise RuntimeError(\n f\"alpha range on ray = [{a_seg.min():.3g}, {a_seg.max():.3g}]; \"\n \"interface not fully captured on this ray\"\n )\n\n order = np.argsort(a_seg)\n x_05 = float(np.interp(0.05, a_seg[order], x_seg[order]))\n x_95 = float(np.interp(0.95, a_seg[order], x_seg[order]))\n t_m = abs(x_95 - x_05)\n dx_local = float(np.mean(np.diff(x_seg)))\n return t_m, t_m / dx_local, x_seg, a_seg, x_05, x_95" }, { "cell_type": "code", "execution_count": null, + "id": "9a63283cbaf04dbcab1f6479b197f3a8", "metadata": {}, "outputs": [], "source": [ "steps = discover_timesteps(CASE_DIR, \"silo\")\n", "if not steps:\n", - " raise RuntimeError(\n", - " f\"No silo timesteps found under {CASE_DIR!r}. \"\n", - " \"Run pre_process + post_process (format=1) first.\"\n", - " )\n", + " raise RuntimeError(f\"No silo timesteps found under {CASE_DIR!r}. Run pre_process + post_process (format=1) first.\")\n", "\n", "step = STEP if STEP is not None else steps[0]\n", "if step not in steps:\n", @@ -135,6 +73,7 @@ { "cell_type": "code", "execution_count": null, + "id": "8dd0d8092fe74a7c96281538738b07e2", "metadata": {}, "outputs": [], "source": [ @@ -142,18 +81,19 @@ " d = np.diff(coord)\n", " return d.mean(), d.min(), d.max()\n", "\n", + "\n", "stats = {\n", " \"x\": (len(x), *_axis_stats(x)),\n", " \"y\": (len(y), *_axis_stats(y)),\n", " \"z\": (len(z), *_axis_stats(z)),\n", "}\n", "\n", - "print(f\"{'axis':<5}{'N':>6}{'\u0394 [um]':>12}{'CPD = D/\u0394':>14}{'nonuniformity':>18}\")\n", + "print(f\"{'axis':<5}{'N':>6}{'Δ [um]':>12}{'CPD = D/Δ':>14}{'nonuniformity':>18}\")\n", "print(\"-\" * 55)\n", "for name, (n, d, dmin, dmax) in stats.items():\n", " u = (dmax - dmin) / max(dmin, 1e-300)\n", " tag = \"uniform\" if u < 1e-6 else f\"{u:.1e}\"\n", - " print(f\"{name:<5}{n:>6d}{d*1e6:>12.4f}{D/d:>14.2f}{tag:>18}\")\n", + " print(f\"{name:<5}{n:>6d}{d * 1e6:>12.4f}{D / d:>14.2f}{tag:>18}\")\n", "\n", "dx = stats[\"x\"][1]\n", "dy = stats[\"y\"][1]\n", @@ -162,6 +102,7 @@ }, { "cell_type": "markdown", + "id": "72eea5119410473aa328ad9291626812", "metadata": {}, "source": [ "## Midplane $\\alpha_1$ field on $z=0$\n", @@ -172,6 +113,7 @@ { "cell_type": "code", "execution_count": null, + "id": "8edb47106e1a46a883d545849b8ab81b", "metadata": {}, "outputs": [], "source": [ @@ -180,10 +122,13 @@ "\n", "fig, ax = plt.subplots(figsize=(10, 6))\n", "pc = ax.pcolormesh(\n", - " x / D, y / D, slice_xy.T,\n", + " x / D,\n", + " y / D,\n", + " slice_xy.T,\n", " shading=\"nearest\",\n", " cmap=\"viridis\",\n", - " vmin=0.0, vmax=1.0,\n", + " vmin=0.0,\n", + " vmax=1.0,\n", " edgecolors=\"k\",\n", " linewidths=0.05,\n", ")\n", @@ -200,6 +145,7 @@ }, { "cell_type": "markdown", + "id": "10185d26023b46108eb7d9f57d49d2b3", "metadata": {}, "source": [ "## Interface close-up\n", @@ -210,6 +156,7 @@ { "cell_type": "code", "execution_count": null, + "id": "8763a12b2bbd4a93a75aff182afb95dc", "metadata": {}, "outputs": [], "source": [ @@ -225,10 +172,13 @@ "\n", "fig, ax = plt.subplots(figsize=(10, 7))\n", "pc = ax.pcolormesh(\n", - " sub_x / D, sub_y / D, sub_alpha.T,\n", + " sub_x / D,\n", + " sub_y / D,\n", + " sub_alpha.T,\n", " shading=\"nearest\",\n", " cmap=\"viridis\",\n", - " vmin=0.0, vmax=1.0,\n", + " vmin=0.0,\n", + " vmax=1.0,\n", " edgecolors=\"k\",\n", " linewidths=0.4,\n", ")\n", @@ -239,8 +189,7 @@ " a_val = float(sub_alpha[i, j])\n", " if 0.05 < a_val < 0.95:\n", " txt_color = \"w\" if a_val < 0.5 else \"k\"\n", - " ax.text(xv / D, yv / D, f\"{a_val:.2f}\",\n", - " ha=\"center\", va=\"center\", fontsize=6, color=txt_color)\n", + " ax.text(xv / D, yv / D, f\"{a_val:.2f}\", ha=\"center\", va=\"center\", fontsize=6, color=txt_color)\n", "\n", "ax.set_xlabel(\"x / D\")\n", "ax.set_ylabel(\"y / D\")\n", @@ -255,6 +204,7 @@ }, { "cell_type": "markdown", + "id": "7623eae2785240b9bd12b16a66d81610", "metadata": {}, "source": [ "## $\\alpha_1$ along an x-ray through each droplet centre\n", @@ -265,65 +215,21 @@ { "cell_type": "code", "execution_count": null, + "id": "7cdc8c89c7104fffa095e18ddfef8986", "metadata": {}, "outputs": [], - "source": [ - "fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\n", - "for ax, (label, (xc, yc, zc)) in zip(axes, DROPLETS.items()):\n", - " try:\n", - " _, _, x_seg, a_seg, x_05, x_95 = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", - " have_crossings = True\n", - " except RuntimeError as err:\n", - " have_crossings = False\n", - " msg = str(err)\n", - " jy = int(np.argmin(np.abs(y - yc)))\n", - " kz = int(np.argmin(np.abs(z - zc)))\n", - " if xc >= 0:\n", - " mask = (x >= xc) & (x <= xc + 2.0 * R)\n", - " else:\n", - " mask = (x <= xc) & (x >= xc - 2.0 * R)\n", - " x_seg = x[mask]\n", - " a_seg = alpha[mask, jy, kz]\n", - "\n", - " xrel = (x_seg - xc) / D\n", - " ax.plot(xrel, a_seg, drawstyle=\"steps-mid\", lw=1.2, color=\"C0\")\n", - " ax.plot(xrel, a_seg, \"o\", color=\"C0\", ms=4)\n", - "\n", - " edges = np.concatenate((\n", - " [x_seg[0] - 0.5 * (x_seg[1] - x_seg[0])],\n", - " 0.5 * (x_seg[:-1] + x_seg[1:]),\n", - " [x_seg[-1] + 0.5 * (x_seg[-1] - x_seg[-2])],\n", - " ))\n", - " for e in edges:\n", - " ax.axvline((e - xc) / D, color=\"k\", lw=0.2, alpha=0.3)\n", - "\n", - " ax.axhline(0.05, color=\"k\", lw=0.5, ls=\":\")\n", - " ax.axhline(0.95, color=\"k\", lw=0.5, ls=\":\")\n", - "\n", - " if have_crossings:\n", - " ax.plot([(x_05 - xc) / D, (x_95 - xc) / D], [0.05, 0.95], \"x\", color=\"r\", ms=10, mew=2,\n", - " label=\"interp crossings\")\n", - " ax.legend(loc=\"best\", fontsize=8)\n", - "\n", - " ax.set_xlabel(r\"$(x - x_c) / D$\")\n", - " ax.set_title(label if have_crossings else f\"{label} \u2014 no 0.05/0.95 crossing on this ray\")\n", - " ax.grid(alpha=0.2)\n", - "\n", - "axes[0].set_ylabel(r\"$\\alpha_1$\")\n", - "fig.suptitle(f\"$\\\\alpha_1$ along x-ray through droplet centre (step {step})\")\n", - "plt.tight_layout()\n", - "plt.show()" - ] + "source": "fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\nfor ax, (label, (xc, yc, zc)) in zip(axes, DROPLETS.items()):\n try:\n _, _, x_seg, a_seg, x_05, x_95 = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n have_crossings = True\n except (RuntimeError, ValueError) as err:\n have_crossings = False\n print(f\" {label:>5}: {err}\")\n jy = int(np.argmin(np.abs(y - yc)))\n kz = int(np.argmin(np.abs(z - zc)))\n if xc >= 0:\n mask = (x >= xc) & (x <= xc + 2.0 * R)\n else:\n mask = (x <= xc) & (x >= xc - 2.0 * R)\n x_seg = x[mask]\n a_seg = alpha[mask, jy, kz] if x_seg.size else np.array([])\n\n if x_seg.size == 0:\n ax.text(0.5, 0.5, \"no cells on ray\", ha=\"center\", va=\"center\", transform=ax.transAxes)\n ax.set_title(f\"{label} \\u2014 empty ray\")\n continue\n\n xrel = (x_seg - xc) / D\n ax.plot(xrel, a_seg, drawstyle=\"steps-mid\", lw=1.2, color=\"C0\")\n ax.plot(xrel, a_seg, \"o\", color=\"C0\", ms=4)\n\n if x_seg.size >= 2:\n edges = np.concatenate((\n [x_seg[0] - 0.5 * (x_seg[1] - x_seg[0])],\n 0.5 * (x_seg[:-1] + x_seg[1:]),\n [x_seg[-1] + 0.5 * (x_seg[-1] - x_seg[-2])],\n ))\n for e in edges:\n ax.axvline((e - xc) / D, color=\"k\", lw=0.2, alpha=0.3)\n\n ax.axhline(0.05, color=\"k\", lw=0.5, ls=\":\")\n ax.axhline(0.95, color=\"k\", lw=0.5, ls=\":\")\n\n if have_crossings:\n ax.plot([(x_05 - xc) / D, (x_95 - xc) / D], [0.05, 0.95], \"x\", color=\"r\", ms=10, mew=2,\n label=\"interp crossings\")\n ax.legend(loc=\"best\", fontsize=8)\n\n ax.set_xlabel(r\"$(x - x_c) / D$\")\n ax.set_title(label if have_crossings else f\"{label} \\u2014 no 0.05/0.95 crossing on this ray\")\n ax.grid(alpha=0.2)\n\naxes[0].set_ylabel(r\"$\\alpha_1$\")\nfig.suptitle(f\"$\\\\alpha_1$ along x-ray through droplet centre (step {step})\")\nplt.tight_layout()\nplt.show()" }, { "cell_type": "markdown", + "id": "b118ea5561624da68c537baed56e602f", "metadata": {}, "source": [ "## Derived: interface thickness\n", "\n", "Two independent estimates of interface thickness, derived from the field you just looked at:\n", - "- **Ray method** \u2014 $|x(\\alpha=0.95) - x(\\alpha=0.05)|$ interpolated along the x-ray through each droplet centre.\n", - "- **$|\\nabla \\alpha|$ method** \u2014 $0.9 / \\max |\\nabla \\alpha|$ over a 1.5R cube around each droplet. Assumes an approximately linear transition.\n", + "- **Ray method** — $|x(\\alpha=0.95) - x(\\alpha=0.05)|$ interpolated along the x-ray through each droplet centre.\n", + "- **$|\\nabla \\alpha|$ method** — $0.9 / \\max |\\nabla \\alpha|$ over a 1.5R cube around each droplet. Assumes an approximately linear transition.\n", "\n", "If these two disagree by more than a factor of ~2, the interface is not well-represented by an axis-aligned ray (e.g. tilted interface, or the ray misses the interface)." ] @@ -331,98 +237,26 @@ { "cell_type": "code", "execution_count": null, + "id": "938c804e27f84196a10c8828c723f798", "metadata": {}, "outputs": [], - "source": [ - "results = {}\n", - "for label, (xc, yc, zc) in DROPLETS.items():\n", - " try:\n", - " t_m, t_c, *_ = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", - " results[label] = (t_m, t_c)\n", - " except RuntimeError as err:\n", - " print(f\" {label:>5}: ray method failed \u2014 {err}\")\n", - " results[label] = (np.nan, np.nan)\n", - "\n", - "print(\"Ray method (0.05 <= alpha_1 <= 0.95):\")\n", - "for k, (t_m, t_c) in results.items():\n", - " if not np.isnan(t_c):\n", - " print(f\" {k:>5}: {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n", - "vals = [v[1] for v in results.values() if not np.isnan(v[1])]\n", - "t_mean_cells = float(np.mean(vals)) if vals else float(\"nan\")\n", - "print(f\" mean : {t_mean_cells:.2f} cells\")\n", - "\n", - "print()\n", - "print(\"|grad alpha| method (max over a 1.5R cube):\")\n", - "for label, (xc, yc, zc) in DROPLETS.items():\n", - " mx = (x >= xc - 1.5 * R) & (x <= xc + 1.5 * R)\n", - " my = (y >= yc - 1.5 * R) & (y <= yc + 1.5 * R)\n", - " mz = (z >= zc - 1.5 * R) & (z <= zc + 1.5 * R)\n", - " a_box = alpha[np.ix_(mx, my, mz)]\n", - " gx_box = np.gradient(a_box, x[mx], axis=0)\n", - " gy_box = np.gradient(a_box, y[my], axis=1)\n", - " gz_box = np.gradient(a_box, z[mz], axis=2)\n", - " gmag = np.sqrt(gx_box**2 + gy_box**2 + gz_box**2)\n", - " t_grad_m = 0.9 / gmag.max()\n", - " dx_local = float(np.mean(np.diff(x[mx])))\n", - " print(f\" {label:>5}: max|grad alpha|={gmag.max():.3g} 1/m -> {t_grad_m*1e6:.2f} um ({t_grad_m/dx_local:.2f} cells)\")" - ] + "source": "results = {}\nfor label, (xc, yc, zc) in DROPLETS.items():\n try:\n t_m, t_c, *_ = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n results[label] = (t_m, t_c)\n except (RuntimeError, ValueError) as err:\n print(f\" {label:>5}: ray method failed \\u2014 {err}\")\n results[label] = (np.nan, np.nan)\n\nprint(\"Ray method (0.05 <= alpha_1 <= 0.95):\")\nfor k, (t_m, t_c) in results.items():\n if not np.isnan(t_c):\n print(f\" {k:>5}: {t_c:.2f} cells = {t_m*1e6:.2f} um\")\nvals = [v[1] for v in results.values() if not np.isnan(v[1])]\nt_mean_cells = float(np.mean(vals)) if vals else float(\"nan\")\nif not np.isnan(t_mean_cells):\n print(f\" mean : {t_mean_cells:.2f} cells\")\n\nprint()\nprint(\"|grad alpha| method (max over a 1.5R cube):\")\nfor label, (xc, yc, zc) in DROPLETS.items():\n mx = (x >= xc - 1.5 * R) & (x <= xc + 1.5 * R)\n my = (y >= yc - 1.5 * R) & (y <= yc + 1.5 * R)\n mz = (z >= zc - 1.5 * R) & (z <= zc + 1.5 * R)\n a_box = alpha[np.ix_(mx, my, mz)]\n if a_box.size == 0:\n print(f\" {label:>5}: ROI has no cells\")\n continue\n gx_box = np.gradient(a_box, x[mx], axis=0)\n gy_box = np.gradient(a_box, y[my], axis=1)\n gz_box = np.gradient(a_box, z[mz], axis=2)\n gmag = np.sqrt(gx_box**2 + gy_box**2 + gz_box**2)\n gmax = float(gmag.max())\n if gmax <= 0:\n print(f\" {label:>5}: |grad alpha| is zero in ROI (no interface here)\")\n continue\n t_grad_m = 0.9 / gmax\n dx_local = float(np.mean(np.diff(x[mx])))\n print(f\" {label:>5}: max|grad alpha|={gmax:.3g} 1/m -> {t_grad_m*1e6:.2f} um ({t_grad_m/dx_local:.2f} cells)\")" }, { "cell_type": "code", "execution_count": null, + "id": "504fb2a444614c0babb325280ed9130a", "metadata": {}, "outputs": [], - "source": [ - "if len(steps) < 2:\n", - " print(\"Only one step available \u2014 skipping time-series.\")\n", - "else:\n", - " series = {k: [] for k in DROPLETS}\n", - " for s in steps:\n", - " d_s = assemble_silo(CASE_DIR, s)\n", - " a_s = d_s.variables[ALPHA_NAME]\n", - " xs, ys, zs = d_s.x_cc, d_s.y_cc, d_s.z_cc\n", - " for label, (xc0, _yc0, _zc0) in DROPLETS.items():\n", - " mx = (xs >= xc0 - 1.5 * R) & (xs <= xc0 + 1.5 * R)\n", - " w = a_s[mx, :, :]\n", - " mass = float(w.sum())\n", - " if mass <= 0:\n", - " series[label].append(np.nan)\n", - " continue\n", - " xw = float(np.sum(w.sum(axis=(1, 2)) * xs[mx]) / mass)\n", - " yw = float(np.sum(w.sum(axis=(0, 2)) * ys) / mass)\n", - " zw = float(np.sum(w.sum(axis=(0, 1)) * zs) / mass)\n", - " try:\n", - " _, t_c, *_ = interface_thickness_x_ray(a_s, xs, ys, zs, xw, yw, zw, R)\n", - " except RuntimeError:\n", - " t_c = np.nan\n", - " series[label].append(t_c)\n", - "\n", - " fig, ax = plt.subplots(figsize=(7, 4))\n", - " for label, vals in series.items():\n", - " ax.plot(steps, vals, \"o-\", label=label)\n", - " ax.set_xlabel(\"step\")\n", - " ax.set_ylabel(\"interface thickness [cells]\")\n", - " ax.set_title(\"Interface thickness over time\")\n", - " ax.legend()\n", - " ax.grid(alpha=0.3)\n", - " plt.show()" - ] + "source": "if len(steps) < 2:\n print(\"Only one step available \\u2014 skipping time-series.\")\nelse:\n series = {k: [] for k in DROPLETS}\n for s in steps:\n d_s = assemble_silo(CASE_DIR, s)\n a_s = d_s.variables[ALPHA_NAME]\n xs, ys, zs = d_s.x_cc, d_s.y_cc, d_s.z_cc\n for label, (xc0, _yc0, _zc0) in DROPLETS.items():\n mx = (xs >= xc0 - 1.5 * R) & (xs <= xc0 + 1.5 * R)\n if not mx.any():\n series[label].append(np.nan)\n continue\n w = a_s[mx, :, :]\n mass = float(w.sum())\n if mass <= 0:\n series[label].append(np.nan)\n continue\n xw = float(np.sum(w.sum(axis=(1, 2)) * xs[mx]) / mass)\n yw = float(np.sum(w.sum(axis=(0, 2)) * ys) / mass)\n zw = float(np.sum(w.sum(axis=(0, 1)) * zs) / mass)\n try:\n _, t_c, *_ = interface_thickness_x_ray(a_s, xs, ys, zs, xw, yw, zw, R)\n except (RuntimeError, ValueError):\n t_c = np.nan\n series[label].append(t_c)\n\n fig, ax = plt.subplots(figsize=(7, 4))\n for label, vals in series.items():\n ax.plot(steps, vals, \"o-\", label=label)\n ax.set_xlabel(\"step\")\n ax.set_ylabel(\"interface thickness [cells]\")\n ax.set_title(\"Interface thickness over time\")\n ax.legend()\n ax.grid(alpha=0.3)\n plt.show()" }, { "cell_type": "code", "execution_count": null, + "id": "59bbdb311c014d738909a11f9e486628", "metadata": {}, "outputs": [], - "source": [ - "print(f\"case: {CASE_LABEL} step: {step}\")\n", - "print(f\"grid: {stats['x'][0]} x {stats['y'][0]} x {stats['z'][0]}\")\n", - "print(f\"CPD : x={D/dx:.1f} y={D/dy:.1f} z={D/dz:.1f}\")\n", - "for k, (t_m, t_c) in results.items():\n", - " if not np.isnan(t_c):\n", - " print(f\"interface thickness ({k:>5}): {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n", - "if not np.isnan(t_mean_cells):\n", - " print(f\"mean interface thickness: {t_mean_cells:.2f} cells\")" - ] + "source": "print(f\"case: {CASE_LABEL} step: {step}\")\nprint(f\"grid: {stats['x'][0]} x {stats['y'][0]} x {stats['z'][0]}\")\nprint(f\"CPD : x={D/dx:.1f} y={D/dy:.1f} z={D/dz:.1f}\")\n_results = locals().get(\"results\", {})\nfor k, (t_m, t_c) in _results.items():\n if not np.isnan(t_c):\n print(f\"interface thickness ({k:>5}): {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n_t_mean = locals().get(\"t_mean_cells\", float(\"nan\"))\nif not np.isnan(_t_mean):\n print(f\"mean interface thickness: {_t_mean:.2f} cells\")" } ], "metadata": { @@ -438,4 +272,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file From e2d2aa015a079f6e75e7f3f1e41918e5afbcd73e Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Thu, 23 Apr 2026 18:10:24 -0400 Subject: [PATCH 09/10] Expand notebook source strings back to line lists after auto-format Co-Authored-By: Claude Opus 4.7 --- .../3D_droplet_coalescence/analysis.ipynb | 226 +++++++++++++++++- 1 file changed, 221 insertions(+), 5 deletions(-) diff --git a/examples/3D_droplet_coalescence/analysis.ipynb b/examples/3D_droplet_coalescence/analysis.ipynb index c629b54f5f..6d00cfde23 100644 --- a/examples/3D_droplet_coalescence/analysis.ipynb +++ b/examples/3D_droplet_coalescence/analysis.ipynb @@ -37,7 +37,82 @@ "id": "acae54e37e7d407bbb7b55eff062a284", "metadata": {}, "outputs": [], - "source": "import sys\nfrom pathlib import Path\n\nNOTEBOOK_DIR = Path.cwd().resolve()\nMFC_ROOT = NOTEBOOK_DIR.parents[1]\nsys.path.insert(0, str(MFC_ROOT / \"toolchain\"))\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom mfc.viz import assemble_silo, discover_timesteps\n\nD = 240e-6\nR = D / 2.0\n\n# Mirror of the CASES table in case.py. Keep in sync if case.py changes.\nCASES = {\n \"a\": (0.2, 14.8, 0.20), \"b\": (0.5, 23.6, 0.10), \"c\": (8.6, 105.9, 0.08),\n \"d\": (15.2, 139.8, 0.08), \"e\": (19.4, 158.0, 0.05), \"f\": (32.8, 210.8, 0.08),\n \"g\": (37.2, 228.0, 0.01), \"h\": (61.4, 296.5, 0.06), \"i\": (61.3, 295.3, 0.11),\n \"j\": (56.3, 288.9, 0.13), \"k\": (70.8, 327.7, 0.25), \"l\": (48.1, 270.1, 0.39),\n \"m\": (60.1, 302.8, 0.55), \"n\": (65.1, 320.3, 0.49), \"o\": (60.8, 313.7, 0.68),\n \"p\": (64.9, 312.8, 0.71), \"q\": (48.8, 260.3, 0.72), \"r\": (14.5, 149.1, 0.34),\n}\n\nCASE_LABEL = \"a\"\n_, _, B = CASES[CASE_LABEL]\nsep = 0.55 * D\nb_impact = B * D\n\nDROPLETS = {\n \"left\": (-sep, -b_impact / 2.0, 0.0),\n \"right\": (+sep, +b_impact / 2.0, 0.0),\n}\n\nCASE_DIR = str(NOTEBOOK_DIR)\nSTEP = None\n\n\ndef interface_thickness_x_ray(alpha, x_cc, y_cc, z_cc, xc, yc, zc, R):\n jy = int(np.argmin(np.abs(y_cc - yc)))\n kz = int(np.argmin(np.abs(z_cc - zc)))\n if xc >= 0:\n mask = (x_cc >= xc) & (x_cc <= xc + 2.0 * R)\n else:\n mask = (x_cc <= xc) & (x_cc >= xc - 2.0 * R)\n\n x_seg = x_cc[mask]\n a_seg = alpha[mask, jy, kz]\n\n if a_seg.size == 0:\n raise RuntimeError(\n f\"x-ray at y={yc:.3g}, z={zc:.3g} through xc={xc:.3g} intersects no grid cells; \"\n \"check CASE_LABEL matches the run\"\n )\n\n if a_seg.min() > 0.05 or a_seg.max() < 0.95:\n raise RuntimeError(\n f\"alpha range on ray = [{a_seg.min():.3g}, {a_seg.max():.3g}]; \"\n \"interface not fully captured on this ray\"\n )\n\n order = np.argsort(a_seg)\n x_05 = float(np.interp(0.05, a_seg[order], x_seg[order]))\n x_95 = float(np.interp(0.95, a_seg[order], x_seg[order]))\n t_m = abs(x_95 - x_05)\n dx_local = float(np.mean(np.diff(x_seg)))\n return t_m, t_m / dx_local, x_seg, a_seg, x_05, x_95" + "source": [ + "import sys\n", + "from pathlib import Path\n", + "\n", + "NOTEBOOK_DIR = Path.cwd().resolve()\n", + "MFC_ROOT = NOTEBOOK_DIR.parents[1]\n", + "sys.path.insert(0, str(MFC_ROOT / \"toolchain\"))\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from mfc.viz import assemble_silo, discover_timesteps\n", + "\n", + "D = 240e-6\n", + "R = D / 2.0\n", + "\n", + "# Mirror of the CASES table in case.py. Keep in sync if case.py changes.\n", + "CASES = {\n", + " \"a\": (0.2, 14.8, 0.20),\n", + " \"b\": (0.5, 23.6, 0.10),\n", + " \"c\": (8.6, 105.9, 0.08),\n", + " \"d\": (15.2, 139.8, 0.08),\n", + " \"e\": (19.4, 158.0, 0.05),\n", + " \"f\": (32.8, 210.8, 0.08),\n", + " \"g\": (37.2, 228.0, 0.01),\n", + " \"h\": (61.4, 296.5, 0.06),\n", + " \"i\": (61.3, 295.3, 0.11),\n", + " \"j\": (56.3, 288.9, 0.13),\n", + " \"k\": (70.8, 327.7, 0.25),\n", + " \"l\": (48.1, 270.1, 0.39),\n", + " \"m\": (60.1, 302.8, 0.55),\n", + " \"n\": (65.1, 320.3, 0.49),\n", + " \"o\": (60.8, 313.7, 0.68),\n", + " \"p\": (64.9, 312.8, 0.71),\n", + " \"q\": (48.8, 260.3, 0.72),\n", + " \"r\": (14.5, 149.1, 0.34),\n", + "}\n", + "\n", + "CASE_LABEL = \"a\"\n", + "_, _, B = CASES[CASE_LABEL]\n", + "sep = 0.55 * D\n", + "b_impact = B * D\n", + "\n", + "DROPLETS = {\n", + " \"left\": (-sep, -b_impact / 2.0, 0.0),\n", + " \"right\": (+sep, +b_impact / 2.0, 0.0),\n", + "}\n", + "\n", + "CASE_DIR = str(NOTEBOOK_DIR)\n", + "STEP = None\n", + "\n", + "\n", + "def interface_thickness_x_ray(alpha, x_cc, y_cc, z_cc, xc, yc, zc, R):\n", + " jy = int(np.argmin(np.abs(y_cc - yc)))\n", + " kz = int(np.argmin(np.abs(z_cc - zc)))\n", + " if xc >= 0:\n", + " mask = (x_cc >= xc) & (x_cc <= xc + 2.0 * R)\n", + " else:\n", + " mask = (x_cc <= xc) & (x_cc >= xc - 2.0 * R)\n", + "\n", + " x_seg = x_cc[mask]\n", + " a_seg = alpha[mask, jy, kz]\n", + "\n", + " if a_seg.size == 0:\n", + " raise RuntimeError(f\"x-ray at y={yc:.3g}, z={zc:.3g} through xc={xc:.3g} intersects no grid cells; check CASE_LABEL matches the run\")\n", + "\n", + " if a_seg.min() > 0.05 or a_seg.max() < 0.95:\n", + " raise RuntimeError(f\"alpha range on ray = [{a_seg.min():.3g}, {a_seg.max():.3g}]; interface not fully captured on this ray\")\n", + "\n", + " order = np.argsort(a_seg)\n", + " x_05 = float(np.interp(0.05, a_seg[order], x_seg[order]))\n", + " x_95 = float(np.interp(0.95, a_seg[order], x_seg[order]))\n", + " t_m = abs(x_95 - x_05)\n", + " dx_local = float(np.mean(np.diff(x_seg)))\n", + " return t_m, t_m / dx_local, x_seg, a_seg, x_05, x_95" + ] }, { "cell_type": "code", @@ -218,7 +293,60 @@ "id": "7cdc8c89c7104fffa095e18ddfef8986", "metadata": {}, "outputs": [], - "source": "fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\nfor ax, (label, (xc, yc, zc)) in zip(axes, DROPLETS.items()):\n try:\n _, _, x_seg, a_seg, x_05, x_95 = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n have_crossings = True\n except (RuntimeError, ValueError) as err:\n have_crossings = False\n print(f\" {label:>5}: {err}\")\n jy = int(np.argmin(np.abs(y - yc)))\n kz = int(np.argmin(np.abs(z - zc)))\n if xc >= 0:\n mask = (x >= xc) & (x <= xc + 2.0 * R)\n else:\n mask = (x <= xc) & (x >= xc - 2.0 * R)\n x_seg = x[mask]\n a_seg = alpha[mask, jy, kz] if x_seg.size else np.array([])\n\n if x_seg.size == 0:\n ax.text(0.5, 0.5, \"no cells on ray\", ha=\"center\", va=\"center\", transform=ax.transAxes)\n ax.set_title(f\"{label} \\u2014 empty ray\")\n continue\n\n xrel = (x_seg - xc) / D\n ax.plot(xrel, a_seg, drawstyle=\"steps-mid\", lw=1.2, color=\"C0\")\n ax.plot(xrel, a_seg, \"o\", color=\"C0\", ms=4)\n\n if x_seg.size >= 2:\n edges = np.concatenate((\n [x_seg[0] - 0.5 * (x_seg[1] - x_seg[0])],\n 0.5 * (x_seg[:-1] + x_seg[1:]),\n [x_seg[-1] + 0.5 * (x_seg[-1] - x_seg[-2])],\n ))\n for e in edges:\n ax.axvline((e - xc) / D, color=\"k\", lw=0.2, alpha=0.3)\n\n ax.axhline(0.05, color=\"k\", lw=0.5, ls=\":\")\n ax.axhline(0.95, color=\"k\", lw=0.5, ls=\":\")\n\n if have_crossings:\n ax.plot([(x_05 - xc) / D, (x_95 - xc) / D], [0.05, 0.95], \"x\", color=\"r\", ms=10, mew=2,\n label=\"interp crossings\")\n ax.legend(loc=\"best\", fontsize=8)\n\n ax.set_xlabel(r\"$(x - x_c) / D$\")\n ax.set_title(label if have_crossings else f\"{label} \\u2014 no 0.05/0.95 crossing on this ray\")\n ax.grid(alpha=0.2)\n\naxes[0].set_ylabel(r\"$\\alpha_1$\")\nfig.suptitle(f\"$\\\\alpha_1$ along x-ray through droplet centre (step {step})\")\nplt.tight_layout()\nplt.show()" + "source": [ + "fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\n", + "for ax, (label, (xc, yc, zc)) in zip(axes, DROPLETS.items()):\n", + " try:\n", + " _, _, x_seg, a_seg, x_05, x_95 = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", + " have_crossings = True\n", + " except (RuntimeError, ValueError) as err:\n", + " have_crossings = False\n", + " print(f\" {label:>5}: {err}\")\n", + " jy = int(np.argmin(np.abs(y - yc)))\n", + " kz = int(np.argmin(np.abs(z - zc)))\n", + " if xc >= 0:\n", + " mask = (x >= xc) & (x <= xc + 2.0 * R)\n", + " else:\n", + " mask = (x <= xc) & (x >= xc - 2.0 * R)\n", + " x_seg = x[mask]\n", + " a_seg = alpha[mask, jy, kz] if x_seg.size else np.array([])\n", + "\n", + " if x_seg.size == 0:\n", + " ax.text(0.5, 0.5, \"no cells on ray\", ha=\"center\", va=\"center\", transform=ax.transAxes)\n", + " ax.set_title(f\"{label} \\u2014 empty ray\")\n", + " continue\n", + "\n", + " xrel = (x_seg - xc) / D\n", + " ax.plot(xrel, a_seg, drawstyle=\"steps-mid\", lw=1.2, color=\"C0\")\n", + " ax.plot(xrel, a_seg, \"o\", color=\"C0\", ms=4)\n", + "\n", + " if x_seg.size >= 2:\n", + " edges = np.concatenate(\n", + " (\n", + " [x_seg[0] - 0.5 * (x_seg[1] - x_seg[0])],\n", + " 0.5 * (x_seg[:-1] + x_seg[1:]),\n", + " [x_seg[-1] + 0.5 * (x_seg[-1] - x_seg[-2])],\n", + " )\n", + " )\n", + " for e in edges:\n", + " ax.axvline((e - xc) / D, color=\"k\", lw=0.2, alpha=0.3)\n", + "\n", + " ax.axhline(0.05, color=\"k\", lw=0.5, ls=\":\")\n", + " ax.axhline(0.95, color=\"k\", lw=0.5, ls=\":\")\n", + "\n", + " if have_crossings:\n", + " ax.plot([(x_05 - xc) / D, (x_95 - xc) / D], [0.05, 0.95], \"x\", color=\"r\", ms=10, mew=2, label=\"interp crossings\")\n", + " ax.legend(loc=\"best\", fontsize=8)\n", + "\n", + " ax.set_xlabel(r\"$(x - x_c) / D$\")\n", + " ax.set_title(label if have_crossings else f\"{label} \\u2014 no 0.05/0.95 crossing on this ray\")\n", + " ax.grid(alpha=0.2)\n", + "\n", + "axes[0].set_ylabel(r\"$\\alpha_1$\")\n", + "fig.suptitle(f\"$\\\\alpha_1$ along x-ray through droplet centre (step {step})\")\n", + "plt.tight_layout()\n", + "plt.show()" + ] }, { "cell_type": "markdown", @@ -240,7 +368,47 @@ "id": "938c804e27f84196a10c8828c723f798", "metadata": {}, "outputs": [], - "source": "results = {}\nfor label, (xc, yc, zc) in DROPLETS.items():\n try:\n t_m, t_c, *_ = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n results[label] = (t_m, t_c)\n except (RuntimeError, ValueError) as err:\n print(f\" {label:>5}: ray method failed \\u2014 {err}\")\n results[label] = (np.nan, np.nan)\n\nprint(\"Ray method (0.05 <= alpha_1 <= 0.95):\")\nfor k, (t_m, t_c) in results.items():\n if not np.isnan(t_c):\n print(f\" {k:>5}: {t_c:.2f} cells = {t_m*1e6:.2f} um\")\nvals = [v[1] for v in results.values() if not np.isnan(v[1])]\nt_mean_cells = float(np.mean(vals)) if vals else float(\"nan\")\nif not np.isnan(t_mean_cells):\n print(f\" mean : {t_mean_cells:.2f} cells\")\n\nprint()\nprint(\"|grad alpha| method (max over a 1.5R cube):\")\nfor label, (xc, yc, zc) in DROPLETS.items():\n mx = (x >= xc - 1.5 * R) & (x <= xc + 1.5 * R)\n my = (y >= yc - 1.5 * R) & (y <= yc + 1.5 * R)\n mz = (z >= zc - 1.5 * R) & (z <= zc + 1.5 * R)\n a_box = alpha[np.ix_(mx, my, mz)]\n if a_box.size == 0:\n print(f\" {label:>5}: ROI has no cells\")\n continue\n gx_box = np.gradient(a_box, x[mx], axis=0)\n gy_box = np.gradient(a_box, y[my], axis=1)\n gz_box = np.gradient(a_box, z[mz], axis=2)\n gmag = np.sqrt(gx_box**2 + gy_box**2 + gz_box**2)\n gmax = float(gmag.max())\n if gmax <= 0:\n print(f\" {label:>5}: |grad alpha| is zero in ROI (no interface here)\")\n continue\n t_grad_m = 0.9 / gmax\n dx_local = float(np.mean(np.diff(x[mx])))\n print(f\" {label:>5}: max|grad alpha|={gmax:.3g} 1/m -> {t_grad_m*1e6:.2f} um ({t_grad_m/dx_local:.2f} cells)\")" + "source": [ + "results = {}\n", + "for label, (xc, yc, zc) in DROPLETS.items():\n", + " try:\n", + " t_m, t_c, *_ = interface_thickness_x_ray(alpha, x, y, z, xc, yc, zc, R)\n", + " results[label] = (t_m, t_c)\n", + " except (RuntimeError, ValueError) as err:\n", + " print(f\" {label:>5}: ray method failed \\u2014 {err}\")\n", + " results[label] = (np.nan, np.nan)\n", + "\n", + "print(\"Ray method (0.05 <= alpha_1 <= 0.95):\")\n", + "for k, (t_m, t_c) in results.items():\n", + " if not np.isnan(t_c):\n", + " print(f\" {k:>5}: {t_c:.2f} cells = {t_m * 1e6:.2f} um\")\n", + "vals = [v[1] for v in results.values() if not np.isnan(v[1])]\n", + "t_mean_cells = float(np.mean(vals)) if vals else float(\"nan\")\n", + "if not np.isnan(t_mean_cells):\n", + " print(f\" mean : {t_mean_cells:.2f} cells\")\n", + "\n", + "print()\n", + "print(\"|grad alpha| method (max over a 1.5R cube):\")\n", + "for label, (xc, yc, zc) in DROPLETS.items():\n", + " mx = (x >= xc - 1.5 * R) & (x <= xc + 1.5 * R)\n", + " my = (y >= yc - 1.5 * R) & (y <= yc + 1.5 * R)\n", + " mz = (z >= zc - 1.5 * R) & (z <= zc + 1.5 * R)\n", + " a_box = alpha[np.ix_(mx, my, mz)]\n", + " if a_box.size == 0:\n", + " print(f\" {label:>5}: ROI has no cells\")\n", + " continue\n", + " gx_box = np.gradient(a_box, x[mx], axis=0)\n", + " gy_box = np.gradient(a_box, y[my], axis=1)\n", + " gz_box = np.gradient(a_box, z[mz], axis=2)\n", + " gmag = np.sqrt(gx_box**2 + gy_box**2 + gz_box**2)\n", + " gmax = float(gmag.max())\n", + " if gmax <= 0:\n", + " print(f\" {label:>5}: |grad alpha| is zero in ROI (no interface here)\")\n", + " continue\n", + " t_grad_m = 0.9 / gmax\n", + " dx_local = float(np.mean(np.diff(x[mx])))\n", + " print(f\" {label:>5}: max|grad alpha|={gmax:.3g} 1/m -> {t_grad_m * 1e6:.2f} um ({t_grad_m / dx_local:.2f} cells)\")" + ] }, { "cell_type": "code", @@ -248,7 +416,44 @@ "id": "504fb2a444614c0babb325280ed9130a", "metadata": {}, "outputs": [], - "source": "if len(steps) < 2:\n print(\"Only one step available \\u2014 skipping time-series.\")\nelse:\n series = {k: [] for k in DROPLETS}\n for s in steps:\n d_s = assemble_silo(CASE_DIR, s)\n a_s = d_s.variables[ALPHA_NAME]\n xs, ys, zs = d_s.x_cc, d_s.y_cc, d_s.z_cc\n for label, (xc0, _yc0, _zc0) in DROPLETS.items():\n mx = (xs >= xc0 - 1.5 * R) & (xs <= xc0 + 1.5 * R)\n if not mx.any():\n series[label].append(np.nan)\n continue\n w = a_s[mx, :, :]\n mass = float(w.sum())\n if mass <= 0:\n series[label].append(np.nan)\n continue\n xw = float(np.sum(w.sum(axis=(1, 2)) * xs[mx]) / mass)\n yw = float(np.sum(w.sum(axis=(0, 2)) * ys) / mass)\n zw = float(np.sum(w.sum(axis=(0, 1)) * zs) / mass)\n try:\n _, t_c, *_ = interface_thickness_x_ray(a_s, xs, ys, zs, xw, yw, zw, R)\n except (RuntimeError, ValueError):\n t_c = np.nan\n series[label].append(t_c)\n\n fig, ax = plt.subplots(figsize=(7, 4))\n for label, vals in series.items():\n ax.plot(steps, vals, \"o-\", label=label)\n ax.set_xlabel(\"step\")\n ax.set_ylabel(\"interface thickness [cells]\")\n ax.set_title(\"Interface thickness over time\")\n ax.legend()\n ax.grid(alpha=0.3)\n plt.show()" + "source": [ + "if len(steps) < 2:\n", + " print(\"Only one step available \\u2014 skipping time-series.\")\n", + "else:\n", + " series = {k: [] for k in DROPLETS}\n", + " for s in steps:\n", + " d_s = assemble_silo(CASE_DIR, s)\n", + " a_s = d_s.variables[ALPHA_NAME]\n", + " xs, ys, zs = d_s.x_cc, d_s.y_cc, d_s.z_cc\n", + " for label, (xc0, _yc0, _zc0) in DROPLETS.items():\n", + " mx = (xs >= xc0 - 1.5 * R) & (xs <= xc0 + 1.5 * R)\n", + " if not mx.any():\n", + " series[label].append(np.nan)\n", + " continue\n", + " w = a_s[mx, :, :]\n", + " mass = float(w.sum())\n", + " if mass <= 0:\n", + " series[label].append(np.nan)\n", + " continue\n", + " xw = float(np.sum(w.sum(axis=(1, 2)) * xs[mx]) / mass)\n", + " yw = float(np.sum(w.sum(axis=(0, 2)) * ys) / mass)\n", + " zw = float(np.sum(w.sum(axis=(0, 1)) * zs) / mass)\n", + " try:\n", + " _, t_c, *_ = interface_thickness_x_ray(a_s, xs, ys, zs, xw, yw, zw, R)\n", + " except (RuntimeError, ValueError):\n", + " t_c = np.nan\n", + " series[label].append(t_c)\n", + "\n", + " fig, ax = plt.subplots(figsize=(7, 4))\n", + " for label, vals in series.items():\n", + " ax.plot(steps, vals, \"o-\", label=label)\n", + " ax.set_xlabel(\"step\")\n", + " ax.set_ylabel(\"interface thickness [cells]\")\n", + " ax.set_title(\"Interface thickness over time\")\n", + " ax.legend()\n", + " ax.grid(alpha=0.3)\n", + " plt.show()" + ] }, { "cell_type": "code", @@ -256,7 +461,18 @@ "id": "59bbdb311c014d738909a11f9e486628", "metadata": {}, "outputs": [], - "source": "print(f\"case: {CASE_LABEL} step: {step}\")\nprint(f\"grid: {stats['x'][0]} x {stats['y'][0]} x {stats['z'][0]}\")\nprint(f\"CPD : x={D/dx:.1f} y={D/dy:.1f} z={D/dz:.1f}\")\n_results = locals().get(\"results\", {})\nfor k, (t_m, t_c) in _results.items():\n if not np.isnan(t_c):\n print(f\"interface thickness ({k:>5}): {t_c:.2f} cells = {t_m*1e6:.2f} um\")\n_t_mean = locals().get(\"t_mean_cells\", float(\"nan\"))\nif not np.isnan(_t_mean):\n print(f\"mean interface thickness: {_t_mean:.2f} cells\")" + "source": [ + "print(f\"case: {CASE_LABEL} step: {step}\")\n", + "print(f\"grid: {stats['x'][0]} x {stats['y'][0]} x {stats['z'][0]}\")\n", + "print(f\"CPD : x={D / dx:.1f} y={D / dy:.1f} z={D / dz:.1f}\")\n", + "_results = locals().get(\"results\", {})\n", + "for k, (t_m, t_c) in _results.items():\n", + " if not np.isnan(t_c):\n", + " print(f\"interface thickness ({k:>5}): {t_c:.2f} cells = {t_m * 1e6:.2f} um\")\n", + "_t_mean = locals().get(\"t_mean_cells\", float(\"nan\"))\n", + "if not np.isnan(_t_mean):\n", + " print(f\"mean interface thickness: {_t_mean:.2f} cells\")" + ] } ], "metadata": { From 8f577d813de5214201b645582c1e33698d8b2b76 Mon Sep 17 00:00:00 2001 From: Davis Cole Date: Wed, 29 Apr 2026 17:31:33 -0400 Subject: [PATCH 10/10] Fix post_process silo coord corruption from assumed-shape re-mapping s_read_grid_data_direction declares dummy arguments dimension(-1:) / dimension(0:). When the actual x_cb / dx / x_cc arrays have lower bounds shifted by offset_x%beg or buff_size (set non-zero for format=1 parallel 3D ranks not at the left domain boundary), Fortran's assumed-shape re-mapping makes the dummy's first element correspond to the actual's first storage slot, not the intended interior index. The read then writes the file data into the wrong range of the actual array and leaves the last offset_x%beg interior cell-boundaries uninitialized, which s_populate_grid_variables_buffers later propagates into the right ghost zone and the silo output. Pass explicit slices x_cb(-1:m), dx(0:m), x_cc(0:m) (and y/z equivalents) so the dummy maps to the intended index range regardless of allocation bounds. Regression from PR #902 (June 2025); pre-refactor code used inline read (1) x_cb(-1:m) which avoided the re-mapping. Verified on a 128-rank 3D droplet-coalescence case: 9/128 ranks previously had NaN / extreme values in coord arrays, now 0/128. Most other ranks were silently shifted by 2 cells and undetectable on uniform grids. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/post_process/m_data_input.f90 | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index e7d1bd3e07..287a8a505c 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -211,13 +211,17 @@ impure subroutine s_read_serial_data_files(t_step) call s_assign_default_bc_type(bc_type) end if - call s_read_grid_data_direction(t_step_dir, 'x', x_cb, dx, x_cc, m) + ! Pass explicit slices so the dummy `dimension(-1:)` / `dimension(0:)` arguments map to the correct interior indices of the + ! actual arrays. Without slicing, when offset_x%beg or buff_size > 0 (i.e. format=1 parallel 3D ranks), Fortran's + ! assumed-shape re-mapping shifts the read by that many slots and leaves the last interior cells uninitialized - corrupting + ! downstream ghost-cell extrapolation. + call s_read_grid_data_direction(t_step_dir, 'x', x_cb(-1:m), dx(0:m), x_cc(0:m), m) if (n > 0) then - call s_read_grid_data_direction(t_step_dir, 'y', y_cb, dy, y_cc, n) + call s_read_grid_data_direction(t_step_dir, 'y', y_cb(-1:n), dy(0:n), y_cc(0:n), n) if (p > 0) then - call s_read_grid_data_direction(t_step_dir, 'z', z_cb, dz, z_cc, p) + call s_read_grid_data_direction(t_step_dir, 'z', z_cb(-1:p), dz(0:p), z_cc(0:p), p) end if end if