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..d01308b95f --- /dev/null +++ b/examples/2D_axisym_droplet_coalescence/case.py @@ -0,0 +1,246 @@ +#!/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 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. + +# 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": "T", + "avg_state": 2, + "riemann_solver": 2, # HLLC + "wave_speeds": 1, + **recon, + # 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, # rectangle (2D) + "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)%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.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)%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.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)%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, +} + +# 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/analysis.ipynb b/examples/3D_droplet_coalescence/analysis.ipynb new file mode 100644 index 0000000000..6d00cfde23 --- /dev/null +++ b/examples/3D_droplet_coalescence/analysis.ipynb @@ -0,0 +1,491 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7fb27b941602401d91542211134fc71a", + "metadata": {}, + "source": [ + "# Droplet-coalescence analysis\n", + "\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", + "Order of cells:\n", + "\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", + "**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, + "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 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", + "execution_count": null, + "id": "9a63283cbaf04dbcab1f6479b197f3a8", + "metadata": {}, + "outputs": [], + "source": [ + "steps = discover_timesteps(CASE_DIR, \"silo\")\n", + "if not steps:\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", + " 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}\")\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", + "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", + "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}]\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8dd0d8092fe74a7c96281538738b07e2", + "metadata": {}, + "outputs": [], + "source": [ + "def _axis_stats(coord):\n", + " 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}{'Δ [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", + "\n", + "dx = stats[\"x\"][1]\n", + "dy = stats[\"y\"][1]\n", + "dz = stats[\"z\"][1]" + ] + }, + { + "cell_type": "markdown", + "id": "72eea5119410473aa328ad9291626812", + "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, + "id": "8edb47106e1a46a883d545849b8ab81b", + "metadata": {}, + "outputs": [], + "source": [ + "kz0 = int(np.argmin(np.abs(z - 0.0)))\n", + "slice_xy = alpha[:, :, kz0]\n", + "\n", + "fig, ax = plt.subplots(figsize=(10, 6))\n", + "pc = ax.pcolormesh(\n", + " x / D,\n", + " y / D,\n", + " slice_xy.T,\n", + " shading=\"nearest\",\n", + " cmap=\"viridis\",\n", + " vmin=0.0,\n", + " 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", + "id": "10185d26023b46108eb7d9f57d49d2b3", + "metadata": {}, + "source": [ + "## Interface close-up\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, + "id": "8763a12b2bbd4a93a75aff182afb95dc", + "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", + "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", + "fig, ax = plt.subplots(figsize=(10, 7))\n", + "pc = ax.pcolormesh(\n", + " sub_x / D,\n", + " sub_y / D,\n", + " sub_alpha.T,\n", + " shading=\"nearest\",\n", + " cmap=\"viridis\",\n", + " vmin=0.0,\n", + " 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}\", ha=\"center\", va=\"center\", fontsize=6, color=txt_color)\n", + "\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", + "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", + "id": "7623eae2785240b9bd12b16a66d81610", + "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, + "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, 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", + "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** — $|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)." + ] + }, + { + "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, 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", + "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", + " 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", + "_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": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/examples/3D_droplet_coalescence/case.py b/examples/3D_droplet_coalescence/case.py new file mode 100644 index 0000000000..c0cab686ba --- /dev/null +++ b/examples/3D_droplet_coalescence/case.py @@ -0,0 +1,200 @@ +#!/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 = 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, + "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)) 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 diff --git a/toolchain/mfc/run/run.py b/toolchain/mfc/run/run.py index 82e886c064..0d3d759b40 100644 --- a/toolchain/mfc/run/run.py +++ b/toolchain/mfc/run/run.py @@ -133,6 +133,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") @@ -141,11 +142,17 @@ 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.")