Moving EL Bubbles with MPI Decomposition

.gitignore

@@ -77,8 +77,7 @@ examples/*/workloads/
 examples/*/run-*/
 examples/*/logs/
 examples/**/*.f90
-!examples/3D_lag_bubbles_shbubcollapse/input/lag_bubbles.dat
-!examples/3D_lag_bubbles_bubblescreen/input/lag_bubbles.dat
+examples/**/*.csv
 workloads/
 
 benchmarks/*batch/*/

docs/documentation/case.md

@@ -276,6 +276,22 @@ Setup: Only requires specifying `files_dir` and filename pattern via `file_exten
 Implementation: All variables and file handling are managed in the `case.py` file of the simulation.
 Usage: Ideal for initializing simulations from lower-dimensional solutions, enabling users to add perturbations or modifications to the base extruded fields for flow instability studies.
 
+The following parameters support hardcoded initial conditions that read interface data from files:
+
+| Parameter        | Type    | Description                                               |
+| ---:             | :---:   | :---                                                      |
+| `interface_file` | String  | Path to interface geometry data file                      |
+| `normFac`        | Real    | Interface normalization factor                            |
+| `normMag`        | Real    | Interface normal magnitude                                |
+| `g0_ic`          | Real    | Gravitational acceleration for the interfacial IC pressure field |
+| `p0_ic`          | Real    | Reference pressure at the interface                       |
+
+These parameters are only read by the file-based hardcoded-IC patches (`hcid = 304` and `305` in `src/common/include/3dHardcodedIC.fpp`); they are ignored otherwise.
+
+- `interface_file` gives the path to a text file that supplies the interface-position field \f$h(i,j)\f$ used to place the material interface. The run aborts if the file is not found.
+- `normMag` and `normFac` rescale and offset the raw interface data, \f$h \leftarrow \texttt{normMag}\,h + \texttt{normFac}\f$. Each is applied only when set (defaults leave the data unchanged).
+- `p0_ic` and `g0_ic` set the initial (hydrostatic) pressure field about the interface, \f$p = p_{0} + \rho\, g_{0}\,\big(h - x\big)\f$, where \f$x\f$ is the coordinate normal to the interface.
+
 #### Parameter Descriptions
 
 - `num_patches` defines the total number of patches defined in the domain.
@@ -930,23 +946,33 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 
 #### 9.2 Volume-Averaged Bubble Model
 
-| Parameter             | Type    | Description                                               |
-| ---:                  | :---:   | :---                                                      |
-| `bubbles_lagrange`    | Logical | Lagrangian subgrid bubble model switch                    |
-| `nBubs_glb`           | Integer | Global number of bubbles                                  |
-| `solver_approach`     | Integer | 1: One-way coupling, 2: two-way coupling                  |
-| `cluster_type`        | Integer | Method to find p_inf                                      |
-| `pressure_corrector`  | Logical | Cell pressure correction term                             |
-| `smooth_type`         | Integer | Smoothing function. 1: Gaussian, 2:Delta 3x3              |
-| `heatTransfer_model`  | Logical | Activates the interface heat transfer model               |
-| `massTransfer_model`  | Logical | Activates the interface mass transfer model               |
-| `write_bubbles`       | Logical | Write files to track the bubble evolution each time step  |
-| `write_bubbles_stats` | Logical | Write the maximum and minimum radius of each bubble       |
-| `epsilonb`            | Real    | Standard deviation scaling for the gaussian function      |
-| `charwidth`           | Real    | Domain virtual depth (z direction, for 2D simulations)    |
-| `valmaxvoid`          | Real    | Maximum void fraction permitted                           |
-
-- `nBubs_glb` Total number of bubbles. Their initial conditions need to be specified in the ./input/lag_bubbles.dat file. See the example cases for additional information.
+| Parameter             | Type    | Description                                                    |
+| ---:                  | :---:   | :---                                                           |
+| `bubbles_lagrange`    | Logical | Lagrangian subgrid bubble model switch                         |
+| `nBubs_glb`           | Integer | Global number of bubbles                                       |
+| `solver_approach`     | Integer | 1: One-way coupling, 2: two-way coupling                       |
+| `cluster_type`        | Integer | Method to find p_inf                                           |
+| `pressure_corrector`  | Logical | Cell pressure correction term                                  |
+| `smooth_type`         | Integer | Smoothing function. 1: Gaussian, 2:Delta 3x3                   |
+| `heatTransfer_model`  | Logical | Activates the interface heat transfer model                    |
+| `massTransfer_model`  | Logical | Activates the interface mass transfer model                    |
+| `write_bubbles`       | Logical | Write files to track the bubble evolution each time step       |
+| `write_bubbles_stats` | Logical | Write the maximum and minimum radius of each bubble            |
+| `write_void_evol`     | Logical | Write stats about the void fraction evolution over time        |
+| `epsilonb`            | Real    | Standard deviation scaling for the gaussian function           |
+| `charwidth`           | Real    | Domain virtual depth (z direction, for 2D simulations)         |
+| `charNz`              | Integer | Number of cells in the virtual depth direction                 |
+| `valmaxvoid`          | Real    | Maximum void fraction permitted                                |
+| `vel_model`           | Integer | Model for translational motion (default 0, disabled)           |
+| `drag_model`          | Integer | Model for drag force (default 0, disabled)                     |
+| `gravity_force`       | Logical | Enable gravity force (default false)                           |
+| `pressure_force`      | Logical | Enable pressure force (default true)                           |
+| `input_path`          | String  | Path to bubble input file (default: `./input/lag_bubbles.dat`) |
+| `kahan_summation`     | Logical | Use Kahan compensated summation when accumulating the void fraction |
+
+- `nBubs_glb` Total number of bubbles. Their initial conditions are read from the file given by `input_path`.
+
+- `input_path` Path to the bubble input file (default `./input/lag_bubbles.dat`). Each row specifies the initial state of one bubble, with columns `xPosition/x0  yPosition/x0  zPosition/x0  xVel/c0  yVel/c0  zVel/c0  radius/x0  interfaceVelocity/c0`. See `examples/3D_lagrange_shbubcollapse/input/lag_bubbles.dat` for a checked-in example, or the other Lagrange example cases (e.g. `examples/2D_lagrange_bubblescreen/case.py`), which generate this file programmatically.
 
 - `solver_approach` Specifies the Euler-Lagrange coupling method: [1] enables a one-way coupling approach, where the bubbles do not influence the Eulerian field. [2] activates the two-way coupling approach based on \cite Maeda18, where the effect of the bubbles is added in the Eulerian field as source terms.
 
@@ -958,6 +984,14 @@ When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due
 
 - `massTransfer_model` Activates the mass transfer model at the bubble's interface based on (\cite Preston07).
 
+- `vel_model` activates translational motion of the bubbles (\cite Wilfong26): [1] tracer bubbles, which are advected with the local carrier velocity \f$\underline{u}\f$ so that \f$\dot{\underline{x}}_b = \underline{u}(\underline{x}_b)\f$; [2] Newton's second law, which integrates \f$m_b \ddot{\underline{x}}_b = \underline{F}_D + \underline{F}_p + \underline{F}_g\f$, where \f$m_b\f$ is the bubble mass and the right-hand side collects the drag, pressure, and gravity forces below. The carrier velocity at the bubble is interpolated with a Lagrange polynomial of order set by `fd_order`, which must be specified when `vel_model > 0`.
+
+    - `drag_model` selects the drag force \f$\underline{F}_D\f$ acting on the slip velocity \f$\underline{u}_{\rm rel} = \underline{u}_b - \underline{u}\f$, with bubble radius \f$a\f$ and Reynolds number \f$Re\f$: [0] no drag (default); [1] free-slip drag (clean-interface creeping-flow limit, \cite Hadamard1911, \cite Rybczynski1911), \f$\underline{F}_D = -4\pi a\,\underline{u}_{\rm rel}/Re\f$; [2] no-slip Stokes drag (rigid sphere, \cite Stokes1851), \f$\underline{F}_D = -6\pi a\,\underline{u}_{\rm rel}/Re\f$; [3] Levich drag (clean bubble at high \f$Re\f$, \cite Levich1962), \f$\underline{F}_D = -12\pi a\,\underline{u}_{\rm rel}/Re\f$. See \cite Magnaudet2000 for a review of these bubble-drag regimes.
+    - `pressure_force` (default true) enables the pressure-gradient force \f$\underline{F}_p = -V_b\,\nabla p\f$, where \f$V_b = \frac{4}{3}\pi a^3\f$ is the bubble volume.
+    - `gravity_force` (default false) enables the body force \f$\underline{F}_g = m_b\,\underline{g}\f$, with \f$\underline{g}\f$ the acceleration set by the body-force parameters.
+
+- `kahan_summation` uses Kahan compensated summation when smearing the bubble contributions onto the Eulerian void fraction, reducing the round-off sensitivity of the accumulation to the summation order. It is not compatible with `--mixed` precision builds.
+
 ### 10. Velocity Field Setup {#sec-velocity-field-setup}
 
 | Parameter              | Type    | Description |
@@ -1286,14 +1320,14 @@ Boundary is at polar angle \f$\theta = \mathrm{atan2}(y - y_{\mathrm{centroid}},
 
 ### Immersed Boundary Patch Types {#immersed-boundary-patch-types}
 
-| #    | Name               | Dim.   |
-| ---: | :----:             | :---   |
-| 2    | 2D Circle          | 2      |
-| 3    | 2D Rectangle       | 2      |
-| 4    | 2D Airfoil         | 2      |
-| 8    | 3D Sphere          | 3      |
+| #    | Name               | Dim.   | Notes                                             |
+| ---: | :----:             | :---:  | :---                                              |
+| 2    | 2D Circle          | 2      |                                                   |
+| 3    | 2D Rectangle       | 2      |                                                   |
+| 4    | 2D Airfoil         | 2      |                                                   |
+| 8    | 3D Sphere          | 3      |                                                   |
 | 10   | 3D Cylinder        | 3      | `length_x` sets the axial length of the cylinder. |
-| 11   | 3D Airfoil         | 3      |
+| 11   | 3D Airfoil         | 3      |                                                   |
 
 ### Acoustic Supports {#acoustic-supports}
 

docs/documentation/equations.md

@@ -522,6 +522,25 @@ with \f$\sigma = \varepsilon_b \max(\Delta x^{1/3}_\text{cell},\;R_\text{bubble}
 
 Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order adaptive Runge-Kutta time integration.
 
+**Translational motion (`vel_model > 0`):**
+
+Bubbles may also translate through the carrier flow. Let \f$\mathbf{x}_b\f$ be the bubble position, \f$\mathbf{u}_b\f$ its velocity, \f$a\f$ its radius, and \f$\mathbf{u}_l(\mathbf{x}_b)\f$ the carrier velocity interpolated to the bubble location (a Lagrange polynomial of order `fd_order`, which must be set when `vel_model > 0`).
+
+- **Tracer bubbles (`vel_model = 1`)** follow the local carrier velocity:
+  \f[\frac{d\mathbf{x}_b}{dt} = \mathbf{u}_l(\mathbf{x}_b).\f]
+
+- **Newton's second law (`vel_model = 2`)** integrates the bubble momentum:
+  \f[m_b\,\frac{d\mathbf{u}_b}{dt} = \mathbf{F}_D + \mathbf{F}_p + \mathbf{F}_g, \qquad \frac{d\mathbf{x}_b}{dt} = \mathbf{u}_b,\f]
+  with bubble mass \f$m_b\f$ and forces acting on the slip velocity \f$\mathbf{u}_\text{rel} = \mathbf{u}_b - \mathbf{u}_l\f$:
+
+| Force | `case.py` control | Non-dimensional form |
+|---|---|---|
+| Drag \f$\mathbf{F}_D\f$ | `drag_model` | \f$-\,c_D\,\pi\,a\,\mathbf{u}_\text{rel} / \text{Re}\f$, with \f$c_D = 4\f$ free-slip (\cite Hadamard1911; \cite Rybczynski1911), \f$c_D = 6\f$ no-slip Stokes (\cite Stokes1851), \f$c_D = 12\f$ Levich (\cite Levich1962) |
+| Pressure \f$\mathbf{F}_p\f$ | `pressure_force` | \f$-V_b\,\nabla p\f$, with bubble volume \f$V_b = \frac{4}{3}\pi a^3\f$ |
+| Gravity \f$\mathbf{F}_g\f$ | `gravity_force` | \f$m_b\,\mathbf{g}\f$, with \f$\mathbf{g}\f$ the body-force acceleration |
+
+Here \f$\text{Re}\f$ is the mixture Reynolds number — the same `fluid_pp%%Re(1)` that scales the viscous stress tensor (@ref sec-two-viscosities) — so the drag scales with the liquid viscosity. The three drag models increase in magnitude free-slip \f$<\f$ no-slip \f$<\f$ Levich; see \cite Magnaudet2000 for a review of these bubble-drag regimes.
+
 ---
 
 ## 7. Fluid-Structure Interaction

docs/references.bib

@@ -670,3 +670,47 @@ @article{Papanastasiou87
   year    = {1987},
   doi     = {10.1122/1.549926}
 }
+
+@article{Stokes1851,
+  author  = {G. G. Stokes},
+  title   = {On the effect of the internal friction of fluids on the motion of pendulums},
+  journal = {Transactions of the Cambridge Philosophical Society},
+  volume  = {9},
+  pages   = {8--106},
+  year    = {1851}
+}
+
+@article{Hadamard1911,
+  author  = {J. Hadamard},
+  title   = {Mouvement permanent lent d'une sph{\`e}re liquide et visqueuse dans un liquide visqueux},
+  journal = {Comptes Rendus de l'Acad{\'e}mie des Sciences},
+  volume  = {152},
+  pages   = {1735--1738},
+  year    = {1911}
+}
+
+@article{Rybczynski1911,
+  author  = {W. Rybczynski},
+  title   = {{\"U}ber die fortschreitende Bewegung einer fl{\"u}ssigen Kugel in einem z{\"a}hen Medium},
+  journal = {Bulletin International de l'Acad{\'e}mie des Sciences de Cracovie, S{\'e}rie A},
+  pages   = {40--46},
+  year    = {1911}
+}
+
+@book{Levich1962,
+  author    = {V. G. Levich},
+  title     = {Physicochemical Hydrodynamics},
+  publisher = {Prentice-Hall},
+  address   = {Englewood Cliffs, NJ},
+  year      = {1962}
+}
+
+@article{Magnaudet2000,
+  author  = {J. Magnaudet and I. Eames},
+  title   = {The motion of high-{R}eynolds-number bubbles in inhomogeneous flows},
+  journal = {Annual Review of Fluid Mechanics},
+  volume  = {32},
+  pages   = {659--708},
+  year    = {2000},
+  doi     = {10.1146/annurev.fluid.32.1.659}
+}

examples/2D_lagrange_bubblescreen/case.py

@@ -1,6 +1,8 @@
 #!/usr/bin/env python3
 import json
 import math
+import os
+import random
 
 # Bubble screen
 # Description: A planar acoustic wave interacts with a bubble cloud
@@ -56,6 +58,75 @@
 
 dt = 7.5e-9  # constant time-step - sec
 
+
+def generate_bubble_cloud():
+    """Generate bubble cloud with log-normal size distribution"""
+    # Bubble properties
+    void_fraction = 4e-5
+    mean_radius = 10e-6  # 10 μm in meters
+    shape_param = 0.3  # shape parameter for log-normal distribution
+
+    # Domain: 5mm x 5mm x 5mm cube centered at origin
+    box_size = 5.0e-3  # 5 mm in meters
+
+    # Convert to nondimensional units
+    mean_radius_nd = mean_radius / x0  # in units of x0
+    box_size_nd = box_size / x0
+
+    # Log-normal distribution parameters
+    # For log-normal: sigma is the shape parameter (std dev of log(r))
+    # mean = exp(mu + sigma^2/2)
+    # Solving: mean_radius = exp(mu + sigma^2/2)
+    sigma = shape_param
+    mu = math.log(mean_radius_nd) - sigma**2 / 2
+
+    # Set random seed for reproducibility
+    random.seed(42)
+
+    # Calculate box volume
+    box_volume = box_size_nd**3
+
+    # Estimate initial number of bubbles
+    # Average volume per bubble: (4/3)*pi*r_mean^3
+    mean_vol_per_bubble = (4.0 / 3.0) * math.pi * math.exp(mu + sigma**2 / 2.0) ** 3
+    n_bubbles_estimate = int(void_fraction * box_volume / mean_vol_per_bubble)
+
+    # Generate bubble radii using log-normal distribution
+    # Python's random.lognormvariate generates from lognormal(mu, sigma)
+    radii = [random.lognormvariate(mu, sigma) for _ in range(n_bubbles_estimate)]
+
+    # Adjust radii to match target void fraction
+    current_void = sum((4.0 / 3.0) * math.pi * r**3 for r in radii) / box_volume
+    if current_void > 0:
+        scale_factor = (void_fraction / current_void) ** (1.0 / 3.0)
+        radii = [r * scale_factor for r in radii]
+
+    n_bubbles = len(radii)
+
+    # Generate random positions in the cube
+    box_half = box_size_nd / 2.0
+    positions = [(random.uniform(-box_half, box_half), random.uniform(-box_half, box_half), random.uniform(-box_half, box_half)) for _ in range(n_bubbles)]
+
+    # Create output directory if needed
+    input_dir = os.path.join(os.path.dirname(__file__), "input")
+    os.makedirs(input_dir, exist_ok=True)
+
+    # Write bubble file
+    bubble_file = os.path.join(input_dir, "lag_bubbles.dat")
+    with open(bubble_file, "w") as f:
+        for i in range(n_bubbles):
+            # Format: x y z vx vy vz radius interface_velocity
+            # All velocities are zero at initialization
+            x, y, z = positions[i]
+            r = radii[i]
+            f.write(f"{x:.6e}\t{y:.6e}\t{z:.6e}\t0.0\t0.0\t0.0\t{r:.6e}\t0.0\n")
+
+    return n_bubbles
+
+
+# Generate bubble cloud
+nBubs = generate_bubble_cloud()
+
 # Configuring case dictionary
 print(
     json.dumps(
@@ -75,7 +146,7 @@
             "dt": dt * (c0 / x0),
             "t_step_start": 0,
             "t_step_stop": 3000,
-            "t_step_save": 500,
+            "t_step_save": 30,
             # Simulation Algorithm Parameters
             "model_eqns": "5eq",
             "time_stepper": "rk3",
@@ -130,7 +201,7 @@
             "bubble_model": "keller_miksis",  # Keller-Miksis model
             "thermal": 3,
             "polytropic": "F",
-            "lag_params%nBubs_glb": 1194,  # Number of bubbles
+            "lag_params%nBubs_glb": nBubs,  # Number of bubbles
             "lag_params%solver_approach": 2,
             "lag_params%cluster_type": 2,
             "lag_params%pressure_corrector": "T",
@@ -141,6 +212,8 @@
             "lag_params%valmaxvoid": 0.9,
             "lag_params%write_bubbles": "F",
             "lag_params%write_bubbles_stats": "F",
+            "lag_params%charwidth": 5e-3 / x0,
+            "lag_params%charNz": Ny,
             # Bubble parameters
             "bub_pp%R0ref": 1.0,
             "bub_pp%p0ref": 1.0,
@@ -169,6 +242,7 @@
             "fluid_pp(2)%gamma": 1.0 / (gam_g - 1.0),
             "fluid_pp(2)%pi_inf": 0.0e00,
             "fluid_pp(2)%Re(1)": 1.0 / (mu_g / (rho0 * c0 * x0)),
-        }
+        },
+        indent=4,
     )
 )