OpenFOAM: dynamic meshes and moving walls

The Arbitrary Lagrangian–Eulerian formulation originates with C. W. Hirt, A. A. Amsden, J. L. Cook, “An arbitrary Lagrangian–Eulerian computing method for all flow speeds”, Journal of Computational Physics 14, 227–253, 1974, doi:10.1016/0021-9991(74)90051-5. In OpenFOAM 14 it appears in four places of the source, all verifiable:

• Mesh flux. In fvMesh.C the movePoints() method computes for every face the volume swept by its displacement and divides it by Δt: this is the meshPhi field, i.e. the face flux ∫ Uₘ·n dA.
• Relative flux. The incompressibleFluid module stores phi as the flux relative to the mesh: after the pressure correction it calls fvc::makeRelative(phi, U) (correctPressure.C). Convection therefore uses (U − Uₘ)·n, exactly the term of the section 2 equation.
• Courant number. basicFluidSolver.C computes “Courant Number” with that relative phi: it is the Courant number of the fluid motion relative to the faces, and maxCo limits the step only on this.
• Correction after motion. In moveMesh.C, if the mesh changed and correctPhi is on (default: on for dynamic meshes) or there was a topology change, the absolute flux is rebuilt from Uf (phi = Sf & Uf), made divergence-free on the new geometry by fv::correctPhi (which solves the pcorr equation, hence the "pcorr.*" entry of fvSolution) and finally made relative.

This is why, at the start of a dynamic-mesh case, the log reports “Constructing face momentum Uf”: the face velocity Uf is an extra state field, written to the time directories and needed to restart the computation.

The class names are those of the src/fvMeshMovers and src/fvMeshTopoChangers folders of the version 14 repository. Sliding interfaces between meshes in relative motion (NCC, Guide 15) are a fourth, complementary tool: the topology does not change, the coupling between faces does.

With constant ρ and U = 0 (fluid at rest) the section 2 equation reduces to a purely geometric relation between the volume change of each cell and the velocity of its faces:

d/dt ∫_V dV = ∮_∂V Uₘ·n dA          (continuous)
(V^(n+1) − V^n)/Δt = Σ_f meshPhi_f   (discrete, implicit Euler)

If the discrete version is not satisfied exactly, a uniform field at rest generates spurious fluxes just because the mesh moves. The principle is due to P. D. Thomas, C. K. Lombard, “Geometric conservation law and its application to flow computations on moving grids”, AIAA Journal 17, 1030–1037, 1979, doi:10.2514/3.61273; for finite volumes it was formulated by I. Demirdžić, M. Perić, “Space conservation law in finite volume calculations of fluid flow”, International Journal for Numerical Methods in Fluids 8, 1037–1050, 1988, doi:10.1002/fld.1650080906. Its link with the stability of ALE schemes is analysed by C. Farhat, P. Geuzaine, C. Grandmont, Journal of Computational Physics 174, 669–694, 2001, doi:10.1006/jcph.2001.6932.

OpenFOAM satisfies it by construction for point motion: meshPhi is not computed from a velocity but from the volume actually swept by each face between the old and new positions, divided by Δt. The sum over a cell therefore returns exactly its volume change. A topology change is a different matter: there fields are transferred from one mesh to another (section 5) and conservation depends on the mapping, not on the GCL.

Reading the tutorial’s blockMeshDict and boundary conditions (section 3), the domain is a cylinder of radius 2.5 mm and length 7.5 mm (between x = −7.5 mm and the fixed wall at x = 0), represented by a 5° wedge. The “cone” is a rigid truncated cone translating along x inside it. The outer wall stays still: the farFieldMoving patch is a portion of the outer surface whose points slide along x with the cone, sweeping no volume. The only open patch is left (totalPressure with p₀ = 0).

An exact result follows: the fluid volume is constant at every instant, so the flow rate through left is zero at every instant. The fluid squeezed between the cone base and the fixed wall flows back along the cone side and fills the region opening in front of the tip. A flow rate measured on left that differs from zero beyond the level of the pressure-solver tolerance is therefore a conservation error, and it is the most severe check the case offers, both during the motion and around the six remapping events.

A convenient comparison scale is the volume swept per second by the flat cone base on the wedge: Q_ref = (θ/2)R²U = (0.08727/2)·(2·10⁻³)²·1 ≈ 1.745·10⁻⁷ m³/s. The following script reads the output of the Qleft function (next practical tip) and reports the maximum of |Q|/Q_ref, the net volume exchanged and the behaviour around each remapping; the tolerance must be decided before running it.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Volume balance of the movingCone case from the Qleft flow-rate monitor (Guide 16).

In movingCone the cone translates rigidly inside a fixed outer boundary and 'left' is the only open
patch: the fluid volume is constant, so the flow rate through 'left' must be zero at every step.
This script reads postProcessing/Qleft/<startTime>/surfaceFieldValue.dat (written by
'#includeFunc patchFlowRate(name=Qleft, patch=left)'), and reports:
  * max |Q| and its time, relative to Q_ref = (wedge/2) R^2 U, the volume swept per second by the
    flat cone base (R = 2 mm, U = 1 m/s, 5 deg wedge by default);
  * the cumulative volume exchanged, integral of Q dt (trapezoidal), relative to the wedge volume
    swept by the cone over one stroke, Q_ref * T/2;
  * the largest |Q| within +-window of each mapping event of the tutorial schedule.
The verdict uses a tolerance chosen in advance (--tol, on max |Q|/Q_ref). Read-only.
"""
import argparse
import math
from pathlib import Path


def read_dat(path):
    t, q = [], []
    for line in path.read_text(encoding='utf-8').splitlines():
        s = line.strip()
        if not s or s.startswith('#'):
            continue
        p = s.split()
        t.append(float(p[0])); q.append(float(p[1]))
    if not t:
        raise SystemExit(f'{path}: no data')
    return t, q


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('dat', type=Path)
    ap.add_argument('--R', type=float, default=2e-3, help='cone base radius [m]')
    ap.add_argument('--U', type=float, default=1.0, help='wall speed [m/s]')
    ap.add_argument('--wedge-deg', type=float, default=5.0)
    ap.add_argument('--f', type=float, default=150.0, help='motion frequency [Hz]')
    ap.add_argument('--events', type=float, nargs='*', default=[0.0015, 0.003, 0.0036667, 0.0051667, 0.0081667, 0.0096667],
                    help='mapping instants [s] (default: tutorial schedule, see movingcone_schedule.py)')
    ap.add_argument('--window', type=float, default=2e-5, help='half-width of the window around each event [s]')
    ap.add_argument('--tol', type=float, required=True, help='tolerance on max |Q|/Q_ref chosen in advance')
    a = ap.parse_args()
    t, q = read_dat(a.dat)
    qref = 0.5 * math.radians(a.wedge_deg) * a.R ** 2 * a.U
    i = max(range(len(q)), key=lambda k: abs(q[k]))
    cum = sum(0.5 * (q[k] + q[k - 1]) * (t[k] - t[k - 1]) for k in range(1, len(t)))
    vstroke = qref / (2 * a.f)
    print(f'{len(t)} samples, t = {t[0]:.6g} .. {t[-1]:.6g} s')
    print(f'Q_ref = (wedge/2) R^2 U = {qref:.4g} m3/s')
    print(f'max |Q| = {abs(q[i]):.4g} m3/s at t = {t[i]:.6g} s  ->  max |Q|/Q_ref = {abs(q[i])/qref:.3g}')
    print(f'net volume through left = {cum:.4g} m3  ->  relative to one stroke Q_ref*T/2 = {cum/vstroke:.3g}')
    for e in a.events:
        w = [abs(q[k]) for k in range(len(t)) if abs(t[k] - e) <= a.window]
        print(f'  event t = {e*1e3:8.4f} ms: ' + (f'max |Q|/Q_ref in +-{a.window*1e6:g} us = {max(w)/qref:.3g} ({len(w)} samples)' if w else 'no samples in window'))
    print('RESULT: ' + ('within' if abs(q[i]) / qref <= a.tol else 'ABOVE') + f' tolerance {a.tol:g} on max |Q|/Q_ref')


if __name__ == '__main__':
    main()

python3 volume_balance.py "$HOME/OpenFOAM/guide16_cases/baseline/postProcessing/Qleft/0/surfaceFieldValue.dat" --tol 1e-3

The tutorial records only a cutting plane. The checks of sections 2, 10 and 15 need the flow rate through left, mean pressures on cone and wall, the force on the cone, probes in the gap and a mesh check at every step. The following file adds them while keeping the original cutPlane:

/*--------------------------------*- C++ -*----------------------------------*\
  Guide 16 - monitoring functions for movingCone (OpenFOAM Foundation 14)
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// The tutorial's own function: p and U on the plane z = 0, VTK at write times
#includeFunc cutPlane

// Flow rate through the only open patch [m3/s on the 5 deg wedge].
// The cone translates rigidly inside a fixed outer boundary, so the fluid
// volume is constant: Qleft must stay at the level set by the pressure
// solver tolerance at every step (see volume_balance.py).
#includeFunc patchFlowRate(name=Qleft, patch=left)

// Area-averaged kinematic pressure on the cone and on the fixed wall [m2/s2]
#includeFunc patchAverage(name=pCone, patch=movingWall, fields=(p))
#includeFunc patchAverage(name=pWall, patch=fixedWall, fields=(p))

// Force on the moving cone (kinematic: rhoInf = 1 gives m4/s2 per wedge;
// multiply by rho and by 360/5 = 72 for the full axisymmetric body)
#includeFunc forcesIncompressible(name=coneForces, patches=(movingWall), rhoInf=1, CofR=(0 0 0), pitchAxis=(0 0 1))

// Probes [m]: in the squeeze gap (x = -0.05 mm stays fluid: the cone base
// never passes x = -0.167 mm) and in front of the cone near the open patch
#includeFunc probes(name=gapProbes, points=((-5e-05 0.0005 0) (-5e-05 0.0015 0) (-0.0074 0.001 0)), fields=(p U))

// Mesh quality at every step in which the mesh moved or was mapped (log only).
// To stop and write the fields at the first failed check use instead:
//   #includeFunc checkMesh(stopAt=writeNow)
#includeFunc checkMesh

// ************************************************************************* //

The gap probes are at x = −0.05 mm: the cone base never passes x = −0.167 mm (section 4), so they always stay in the fluid. On the 5° wedge, forces and flow rates refer to 1/72 of the axisymmetric body: multiply by 72 for the real value and, for dynamic quantities, by the density. The checkMesh function (checkMesh.H) runs the full mesh check in every step in which the mesh moved or was remapped; with stopAt writeNow it stops the run and writes the fields at the first failed check.

The script installs the file in a new copy and enables the printing of the mesh Courant number, which in OpenFOAM 14 appears in the log only with checkMeshCourantNo yes in the PIMPLE block (default: no; basicFluidSolver::meshCourantNo()). Run it after prepare_cases.py and before run_case.sh:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Add the Guide 16 monitors to a NEW, never-run copy of the movingCone tutorial.

* system/functions is replaced by functions_guide16 (the tutorial file is kept as
  system/functions.tutorial, never deleted);
* 'checkMeshCourantNo yes;' is added to the PIMPLE dictionary of system/fvSolution, so that
  the solver prints 'Mesh Courant Number mean: ... max: ...' after every mesh motion
  (OpenFOAM 14 prints it only when this switch is on; default off).
The script refuses a case that already has logs, time directories or generated meshes.
"""
import argparse
import re
import shutil
import sys
from pathlib import Path


def already_run(case):
    if list(case.glob('log.*')) or (case / 'run.full.log').exists() or list(case.glob('processor*')):
        return True
    if (case / 'constant' / 'polyMesh').exists() or (case / 'constant' / 'meshes').exists():
        return True
    return any(d.is_dir() and re.fullmatch(r'\d+(\.\d+)?(e-?\d+)?', d.name) and d.name != '0' for d in case.iterdir())


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('case', type=Path)
    ap.add_argument('--functions', type=Path, default=Path(__file__).with_name('functions_guide16'))
    a = ap.parse_args()
    case = a.case.expanduser().resolve()
    for p in ('system/functions', 'system/fvSolution', 'system/controlDict'):
        if not (case / p).is_file():
            sys.exit(f'ERROR: {case/p} missing: not a movingCone case')
    if already_run(case):
        sys.exit('ERROR: the case has logs, results or generated meshes: prepare a NEW copy')
    if not a.functions.is_file():
        sys.exit(f'ERROR: {a.functions} not found')
    keep = case / 'system' / 'functions.tutorial'
    if keep.exists():
        sys.exit(f'ERROR: {keep} exists: monitors already added?')
    fvs = case / 'system' / 'fvSolution'
    t = fvs.read_text(encoding='utf-8')
    if 'checkMeshCourantNo' in t:
        sys.exit('ERROR: checkMeshCourantNo already present in fvSolution')
    m = re.search(r'^PIMPLE\s*\{', t, re.M)
    if not m:
        sys.exit('ERROR: PIMPLE dictionary not found in fvSolution')
    t = t[:m.end()] + '\n    checkMeshCourantNo  yes;' + t[m.end():]
    shutil.copy2(case / 'system' / 'functions', keep)
    (case / 'system' / 'functions').write_text(a.functions.read_text(encoding='utf-8'), encoding='utf-8', newline='\n')
    fvs.write_text(t, encoding='utf-8', newline='\n')
    with (case / 'GUIDE16_CASE.txt').open('a', encoding='utf-8') as f:
        f.write('Monitors: system/functions <- functions_guide16 (tutorial copy in functions.tutorial); '
                'PIMPLE checkMeshCourantNo yes\n')
    print(f'Monitors added to {case}')


if __name__ == '__main__':
    main()

$ python3 enable_monitors.py guide16_cases/courant_half
Monitors added to guide16_cases/courant_half
$ python3 enable_monitors.py guide16_cases/baseline     # already run
ERROR: the case has logs, results or generated meshes: prepare a NEW copy
$ sed -n "/^PIMPLE/,/^}/p" guide16_cases/courant_half/system/fvSolution
PIMPLE
{
    checkMeshCourantNo  yes;
    correctPhi          yes;
    nOuterCorrectors    2;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

The vertices are in millimetres (units [mm]). The domain is a 5° wedge (front and back patches of type wedge, empty axis): an axisymmetric case. The moving body (movingWall) has a flat front face of radius 0.75 mm, a conical side up to radius 2 mm and a flat base of radius 2 mm; the fixed wall fixedWall is the plane x = 0. The three meshes differ only in the body position and in the number of cells along x of the blocks in front of and behind the cone:

The displacement between meshes, 1.5 and 3.0 mm, coincides with that of the imposed motion at times 0.0015 and 0.003 s (section 4): the remapping meshes are the same geometry with the cone already moved and the cells redistributed, not different geometries. The point of the operation is in the “gap cells” column: as the gap closes from 3.5 to 0.5 mm, the cells go from 20 to 5, preventing the motion from squashing them to degenerate thicknesses.

The following script derives the table from the files and adds order-of-magnitude estimates useful for reading the log (gap Reynolds number, squeeze pressure, Courant number of the wall motion):

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Geometry, cell counts and order-of-magnitude estimates for the OpenFOAM 14 movingCone case.

Reads system/blockMeshDict and system/meshes/*/blockMeshDict of a case (or of the tutorial),
reports for each mesh: cone position, gap between the cone base and the fixed wall, cells per block
and in the gap, wedge angle and total cells. Then, with the motion of 0/pointMotionUx (square wave,
amplitude 1 m/s, 150 Hz, read from the defaults below or from the command line), it estimates:
peak displacement, minimum gap, gap Reynolds number, the classical squeeze-film (Stefan) kinematic
pressure p/rho ~ 3 nu V R^2 / h^3 for the flat cone base, and the cell Courant numbers of the wall motion.
These are ESTIMATES to plan the run and to judge the log, not CFD results. Read-only: writes nothing.
"""
import argparse
import math
import os
import re
from pathlib import Path

NUM = r'[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?'


def parse_block_mesh(path):
    t = path.read_text(encoding='utf-8')
    t = re.sub(r'//.*', '', t)
    scale = 1.0
    m = re.search(r'^\s*units\s+\[(\w+)\]\s*;', t, re.M)
    if m:
        scale = {'mm': 1e-3, 'cm': 1e-2, 'm': 1.0, 'um': 1e-6}[m.group(1)]
    m = re.search(r'^\s*convertToMeters\s+(' + NUM + r')\s*;', t, re.M)
    if m:
        scale = float(m.group(1))
    vb = re.search(r'vertices\s*\((.*?)\n\s*\)\s*;', t, re.S).group(1)
    verts = [tuple(float(v) * scale for v in p.split()) for p in re.findall(r'\(([^()]*)\)', vb)]
    blocks = [(tuple(int(i) for i in h.split()), tuple(int(n) for n in c.split()))
              for h, c in re.findall(r'hex\s*\(([^)]*)\)\s*\(([^)]*)\)', t)]
    return verts, blocks


def describe(verts, blocks):
    xs = lambda b: sorted({verts[i][0] for i in b[0]})
    out = []
    for k, b in enumerate(blocks):
        x = xs(b)
        out.append({'block': k, 'x0': x[0], 'x1': x[-1], 'cells': b[1][0] * b[1][1] * b[1][2], 'n': b[1]})
    return out


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    default = Path(os.environ['FOAM_TUTORIALS']) / 'incompressibleFluid' / 'movingCone' if os.environ.get('FOAM_TUTORIALS') else None
    ap.add_argument('case', type=Path, nargs='?', default=default, help='case or tutorial directory')
    ap.add_argument('--U', type=float, default=1.0, help='wall speed amplitude [m/s] (pointMotionUx amplitude)')
    ap.add_argument('--f', type=float, default=150.0, help='square-wave frequency [Hz]')
    ap.add_argument('--nu', type=float, default=1e-5, help='kinematic viscosity [m2/s] (constant/physicalProperties)')
    ap.add_argument('--dt', type=float, default=5e-6, help='time step for the Courant estimates [s]')
    a = ap.parse_args()
    if a.case is None:
        ap.error('give the case directory or set FOAM_TUTORIALS')
    files = [('0 (system)', a.case / 'system' / 'blockMeshDict')]
    md = a.case / 'system' / 'meshes'
    if md.is_dir():
        files += [(d.name, d / 'blockMeshDict') for d in sorted(md.iterdir(), key=lambda p: float(p.name)) if (d / 'blockMeshDict').is_file()]

    print('Mesh            cone base x [mm]  gap to wall [mm]  gap cells  front cells  total cells  wedge [deg]')
    res = {}
    for name, f in files:
        v, b = parse_block_mesh(f)
        d = describe(v, b)
        xmax_wall = max(p[0] for p in v)               # fixed wall at x = 0
        base = d[1]['x0']                              # block 1 = cone base .. wall (gap block)
        gap = xmax_wall - base
        ymax = max(p[1] for p in v); zmax = max(p[2] for p in v)
        wedge = 2 * math.degrees(math.atan2(zmax, ymax))
        total = sum(x['cells'] for x in d)
        res[name] = gap
        print(f'{name:14s} {base*1e3:16.3f} {gap*1e3:17.3f} {d[1]["n"][0]:10d} {d[0]["n"][0]:12d} {total:12d} {wedge:12.3f}')

    T = 1 / a.f
    s_peak = a.U * T / 2
    g0 = res['0 (system)']
    gmin = g0 - s_peak
    R = 2e-3            # radius of the flat cone base (vertices (-3.5 2 .) and (0 2 .) in mm)
    print()
    print(f'Square wave: period {T*1e3:.4f} ms, peak displacement U*T/2 = {s_peak*1e3:.4f} mm at t = {T/2*1e3:.4f} ms')
    print(f'Minimum gap cone base - wall = {g0*1e3:.3f} - {s_peak*1e3:.4f} = {gmin*1e3:.4f} mm (on mesh 0.003: 5 cells -> {gmin/5*1e6:.1f} um mean thickness)')
    for h in (g0, 2.0e-3, 0.5e-3, gmin):
        Re = a.U * h / a.nu
        pk = 3 * a.nu * a.U * R**2 / h**3
        print(f'  gap {h*1e3:7.4f} mm: Re_h = U h / nu = {Re:7.1f};  squeeze estimate p/rho ~ 3 nu U R^2/h^3 = {pk:9.3g} m2/s2;'
              f'  radial exit speed ~ U R/(2h) = {a.U*R/(2*h):6.2f} m/s')
    print()
    print('Wall-motion Courant number U dt / dx for the time step given:')
    for h, n, lab in ((g0, 20, 'initial gap, 20 cells'), (0.5e-3, 5, 'mesh 0.003 at mapping, 5 cells'), (gmin, 5, 'peak, 5 cells')):
        dx = h / n
        print(f'  {lab:32s} dx_mean = {dx*1e6:7.1f} um  ->  U dt/dx = {a.U*a.dt/dx:6.3f}')
    print()
    print('Notes: the squeeze formula is the lubrication (Stefan) result for a flat disc approaching a wall with a free edge;')
    print('here Re_h is not small and the fluid leaves through a 0.5 mm annulus, so it gives only the order of magnitude.')


if __name__ == '__main__':
    main()

$ python3 g16_estimates.py "$FOAM_TUTORIALS/incompressibleFluid/movingCone"
Mesh            cone base x [mm]  gap to wall [mm]  gap cells  front cells  total cells  wedge [deg]
0 (system)               -3.500             3.500         20            8         1690        4.995
0.0015                   -2.000             2.000         10           16         1580        4.995
0.003                    -0.500             0.500          5           25         1675        4.995

Square wave: period 6.6667 ms, peak displacement U*T/2 = 3.3333 mm at t = 3.3333 ms
Minimum gap cone base - wall = 3.500 - 3.3333 = 0.1667 mm (on mesh 0.003: 5 cells -> 33.3 um mean thickness)
  gap  3.5000 mm: Re_h = U h / nu =   350.0;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =    0.0028 m2/s2;  radial exit speed ~ U R/(2h) =   0.29 m/s
  gap  2.0000 mm: Re_h = U h / nu =   200.0;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =     0.015 m2/s2;  radial exit speed ~ U R/(2h) =   0.50 m/s
  gap  0.5000 mm: Re_h = U h / nu =    50.0;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =      0.96 m2/s2;  radial exit speed ~ U R/(2h) =   2.00 m/s
  gap  0.1667 mm: Re_h = U h / nu =    16.7;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =      25.9 m2/s2;  radial exit speed ~ U R/(2h) =   6.00 m/s

Wall-motion Courant number U dt / dx for the time step given:
  initial gap, 20 cells            dx_mean =   175.0 um  ->  U dt/dx =  0.029
  mesh 0.003 at mapping, 5 cells   dx_mean =   100.0 um  ->  U dt/dx =  0.050
  peak, 5 cells                    dx_mean =    33.3 um  ->  U dt/dx =  0.150

Notes: the squeeze formula is the lubrication (Stefan) result for a flat disc approaching a wall with a free edge;
here Re_h is not small and the fluid leaves through a 0.5 mm annulus, so it gives only the order of magnitude.

The squeeze formula is the classical lubrication (Stefan) result for a flat disc approaching a wall: p ≈ 3μU(R² − r²)/h³ at the centre. Here the gap Reynolds number is not small and the fluid leaves through a 0.5 mm annulus between cone and outer wall, so the value must be read only as an order of magnitude: it says that the kinematic pressure in the gap grows by about four orders of magnitude between the start and the point of closest approach, which is where to look first at fields and probes.

File contents (without headers), to be compared with the copy on the server before any change. The full source is in the version 14 repository.

0/U

dimensions      [velocity];

internalField   uniform (0 0 0);

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    movingWall
    {
        type            movingWallVelocity;
        value           $internalField;
    }

    farFieldMoving
    {
        type            noSlip;
    }

    fixedWall
    {
        type            noSlip;
    }

    left
    {
        type            pressureInletOutletVelocity;
    }

    farField
    {
        type            noSlip;
    }
}

0/p

dimensions      [kinematicPressure];

internalField   uniform 0;

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    movingWall
    {
        type            zeroGradient;
    }

    farFieldMoving
    {
        type            zeroGradient;
    }

    fixedWall
    {
        type            zeroGradient;
    }

    left
    {
        type            totalPressure;
        p0              $internalField;
    }

    farField
    {
        type            zeroGradient;
    }
}

0/pointMotionUx

dimensions      [velocity];

internalField   uniform 0;

movement
{
    type            uniformFixedValue;

    uniformValue
    {
        type            square;

        amplitude       1;
        frequency       150;
        start           0;
        level           0;
    }
}


boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    movingWall
    {
        $movement;
    }

    farFieldMoving
    {
        $movement;
    }

    fixedWall
    {
        type            uniformFixedValue;
        uniformValue    constant 0;
    }

    left
    {
        type            uniformFixedValue;
        uniformValue    constant 0;
    }

    farField
    {
        type            slip;
    }
}

constant/physicalProperties, constant/momentumTransport

viscosityModel  constant;

nu      1e-05;
simulationType  laminar;

system/fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss linearUpwind grad(U);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;

    laplacian(diffusivity,cellMotionU) Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

system/fvSolution

solvers
{
    p
    {
        solver          GAMG;
        smoother        DIC;
        cacheAgglomeration no;

        tolerance       0;
        relTol          0.01;
    }

    pFinal
    {
        $p;
        tolerance       1e-06;
        relTol          0;
    }

    "pcorr.*"
    {
        $p
        tolerance       0.02;
        relTol          0;
    }

    U
    {
        solver          PBiCGStab;
        preconditioner  DILU;

        tolerance       1e-6;
        relTol          0.1;
    }

    UFinal
    {
        $U;
        relTol          0;
    }

    cellMotionUx
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    correctPhi          yes;
    nOuterCorrectors    2;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    equations
    {
        ".*"            1;
    }
}

system/functions + system/cutPlane

#includeFunc cutPlane

type            surfaces;
libs            ("libsampling.so");

writeControl    writeTime;

surfaceFormat   vtk;

fields          (p U);

interpolationScheme cellPoint;

surfaces
{
    zNormal
    {
        type            cutPlane;
        point           (0 0 0);
        normal          (0 0 1);
        interpolate     true;
    }
}

Allrun-parallel

#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

runApplication blockMesh
runApplication decomposePar -cellProc

for mapTime in 0.0015 0.003
do
    runApplication -s $mapTime blockMesh -mesh $mapTime
    ln -s meshes/$mapTime/decomposeParDict system/decomposeParDict
    runApplication -s $mapTime decomposePar -mesh $mapTime
    rm -f system/decomposeParDict
done

runParallel foamRun
runApplication reconstructPar -cellProc

#------------------------------------------------------------------------------

Three details that the guide text summarises: on farFieldMoving the fluid velocity is noSlip (the outer wall is still, the mesh points slide over it); on left the pressureInletOutletVelocity/totalPressure pair lets fluid enter or leave at zero total pressure; the kinematic viscosity is ν = 1·10⁻⁵ m²/s in laminar flow. Allrun-parallel also decomposes the two remapping meshes separately, each with its own decomposeParDict (scotch method, 4 subdomains).

In the version 14 source (SquareI.H) the square function is a·s(f(t − t₀)) + l, with s = +1 for the period fraction markSpace/(1 + markSpace) and −1 for the rest; markSpace is 1 if not given. With the tutorial values (a = 1, f = 150, t₀ = 0, l = 0):

• point velocity +1 m/s for 3.333 ms, then −1 m/s for 3.333 ms: period 6.667 ms;
• triangular displacement, maximum 3.333 mm at t = 3.333 ms, zero again at t = 6.667 ms;
• gap between cone base and wall: from 3.5 mm down to 0.167 mm at the point of closest approach;
• in 9.9 ms (endTime) the cone completes about one and a half periods: it is moving forward again, at 3 mm, at the end of the run.

The following script reproduces outside OpenFOAM the square function and the meshToMesh logic (section 5) and computes every remapping event, with the displacement at that instant, cross-checking them with an independent scan:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Motion and mesh-to-mesh mapping schedule of the OpenFOAM 14 movingCone tutorial.

Reproduces, outside OpenFOAM, two pieces of logic read in the OpenFOAM 14 source:
  * Function1 'square' (SquareI.H): value = amplitude*(+1 if fraction < markSpace/(1+markSpace) else -1) + level,
    fraction = frac(frequency*(t - start));
  * fvMeshTopoChangers::meshToMesh (meshToMesh_fvMeshTopoChanger.C): with 'cycle' the mesh sequence is
    walked forward in even cycles and backward in odd cycles; a mesh is loaded when the mesh time,
    rounded to 'timeDelta', equals one of 'times'.
It prints the wall velocity, the integrated cone displacement and every mapping event up to endTime,
with the displacement at that instant, so that each event can be checked against the x-offset of the
corresponding system/meshes/<time>/blockMeshDict. Defaults are the tutorial values; nothing is written
to disk unless --csv is given (and an existing file is never overwritten).
"""
import argparse
import csv
import math
from pathlib import Path


def square(t, amplitude, frequency, start=0.0, level=0.0, mark_space=1.0):
    mark_fraction = mark_space / (1.0 + mark_space)
    phi = frequency * (t - start)
    fraction = phi - math.floor(phi)
    return amplitude * (1.0 if fraction < mark_fraction else -1.0) + level


def square_integral(t, amplitude, frequency, start=0.0, level=0.0, mark_space=1.0):
    """Exact integral of square() from start to t (displacement for a velocity signal)."""
    mark_fraction = mark_space / (1.0 + mark_space)
    phi = frequency * (t - start)
    n = math.floor(phi)
    fraction = phi - n
    # each full period contributes (mark - space) * amplitude / frequency
    per_period = (mark_fraction - (1.0 - mark_fraction)) / frequency
    part = (fraction if fraction < mark_fraction else mark_fraction - (fraction - mark_fraction)) / frequency
    return amplitude * (n * per_period + part) + level * (t - start)


def mesh_time(t, begin, cycle):
    if cycle <= 0:
        return t
    k = int((t - begin) / cycle)
    r = math.fmod(t - begin, cycle)
    return begin + r if k % 2 == 0 else begin + cycle - r


def events(times, begin, cycle, time_delta, end_time):
    """User times at which the mesh time equals a listed mesh time.

    With 'cycle' the mesh time runs forward in even cycles (meshTime = begin + r) and backward in odd
    cycles (meshTime = begin + cycle - r), r = fmod(t - begin, cycle); the solver lands on these instants
    because meshToMeshAdjustTimeStep limits the time step with timeToNextMesh().
    """
    out = []
    if cycle <= 0:
        return [(m, m, True) for m in sorted(times) if m <= end_time + time_delta / 2]
    k = 0
    while begin + k * cycle <= end_time:
        fwd = (k % 2 == 0)
        seq = sorted(times) if fwd else sorted(times, reverse=True)
        for m in seq:
            r = (m - begin) if fwd else (cycle - (m - begin))
            if r < -time_delta / 2 or r > cycle + time_delta / 2:
                continue
            t = begin + k * cycle + r
            if t <= end_time + time_delta / 2 and (not out or t > out[-1][0] + time_delta / 2):
                out.append((t, m, fwd))
        k += 1
    return out


def check_scan(times, begin, cycle, time_delta, end_time, dt=1e-7):
    """Independent check: scan meshTime() with a small step and count distinct hits of the listed times."""
    targets = {int((m + time_delta / 2.0) / time_delta) for m in times}
    hits, last = 0, None
    for i in range(1, int(round(end_time / dt)) + 1):
        idx = int((mesh_time(i * dt, begin, cycle) + time_delta / 2.0) / time_delta)
        if idx in targets and idx != last:
            hits += 1
        last = idx if idx in targets else None
    return hits


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('--amplitude', type=float, default=1.0, help='square amplitude [m/s]')
    ap.add_argument('--frequency', type=float, default=150.0, help='square frequency [Hz]')
    ap.add_argument('--times', type=float, nargs='+', default=[0.0015, 0.003], help='meshToMesh times [s]')
    ap.add_argument('--cycle', type=float, default=1.0 / 300.0, help='meshToMesh cycle [s]')
    ap.add_argument('--begin', type=float, default=0.0)
    ap.add_argument('--time-delta', type=float, default=1e-6)
    ap.add_argument('--end-time', type=float, default=0.0099)
    ap.add_argument('--csv', type=Path, help='optional CSV of t, U_wall, x_wall (new file only)')
    a = ap.parse_args()

    T = 1.0 / a.frequency
    print(f'Square wave: period {T*1e3:.4f} ms, +{a.amplitude} m/s for {T/2*1e3:.4f} ms, then -{a.amplitude} m/s')
    print(f'Peak displacement amplitude*T/2 = {a.amplitude*T/2*1e3:.4f} mm at t = {T/2*1e3:.4f} ms')
    print(f'meshToMesh cycle {a.cycle*1e3:.4f} ms (= T/2: {abs(a.cycle - T/2) < 1e-12}); forward in even cycles, backward in odd cycles')
    print()
    print(f'{"event":>5} {"t [ms]":>9} {"mesh":>8} {"pass":>8} {"x_wall [mm]":>12} {"U_wall [m/s]":>12}')
    ev = events(a.times, a.begin, a.cycle, a.time_delta, a.end_time)
    for k, (t, m, fwd) in enumerate(ev, 1):
        x = square_integral(t, a.amplitude, a.frequency)
        u = square(t, a.amplitude, a.frequency)
        print(f'{k:5d} {t*1e3:9.4f} {m:8g} {"forward" if fwd else "backward":>8} {x*1e3:12.4f} {u:12.1f}')
    print(f'Independent scan of meshTime(): {check_scan(a.times, a.begin, a.cycle, a.time_delta, a.end_time)} events (table: {len(ev)})')
    print()
    for m in a.times:
        print(f'mesh {m:g}: displacement expected at that mesh time = {square_integral(m, a.amplitude, a.frequency)*1e3:.4f} mm')
    if a.csv:
        if a.csv.exists():
            raise SystemExit(f'Refusing to overwrite {a.csv}')
        with a.csv.open('x', newline='', encoding='utf-8') as f:
            w = csv.writer(f)
            w.writerow(('time_s', 'U_wall_m_s', 'x_wall_m'))
            n = 400
            for i in range(n + 1):
                t = a.end_time * i / n
                w.writerow((f'{t:.8g}', square(t, a.amplitude, a.frequency), f'{square_integral(t, a.amplitude, a.frequency):.8g}'))
        print(f'Wrote {a.csv}')


if __name__ == '__main__':
    main()

$ python3 movingcone_schedule.py
Square wave: period 6.6667 ms, +1.0 m/s for 3.3333 ms, then -1.0 m/s
Peak displacement amplitude*T/2 = 3.3333 mm at t = 3.3333 ms
meshToMesh cycle 3.3333 ms (= T/2: True); forward in even cycles, backward in odd cycles

event    t [ms]     mesh     pass  x_wall [mm] U_wall [m/s]
    1    1.5000   0.0015  forward       1.5000          1.0
    2    3.0000    0.003  forward       3.0000          1.0
    3    3.6667    0.003 backward       3.0000         -1.0
    4    5.1667   0.0015 backward       1.5000         -1.0
    5    8.1667   0.0015  forward       1.5000          1.0
    6    9.6667    0.003  forward       3.0000          1.0
Independent scan of meshTime(): 6 events (table: 6)

mesh 0.0015: displacement expected at that mesh time = 1.5000 mm
mesh 0.003: displacement expected at that mesh time = 3.0000 mm

The displacement at each event is exactly that of the mesh loaded (1.5 or 3.0 mm, section 3): this is the condition that makes remapping possible without geometry “jumps”. The script accepts different amplitude, frequency, times and cycle: use it to redesign the sequence before changing any single parameter.

The four quantities are linked: the topoChanger cycle must equal half the motion period (1/300 s = (1/150)/2), and the displacements encoded in the blockMeshDict files of system/meshes/<t> must equal U·t for each time t. For a different amplitude, in a new copy:

• choose the new remapping times so that U·t gives positions with a still reasonable gap and that they fall within half a period;
• create a system/meshes/<t> folder for each time, shifting the x coordinates of the cone vertices by U·t and redistributing the cells between the blocks in front and behind;
• update times and cycle in dynamicMeshDict and the for mapTime in ... list of Allrun;
• check with movingcone_schedule.py --amplitude ... --times ... --cycle ... that each event falls on the displacement of the corresponding mesh.

As an alternative to cycle (sequence walked forward and backward, suited to a symmetric motion) the topoChanger accepts repeat (sequence always repeated forward, suited to a cyclic motion such as an engine’s); specifying both is a fatal error.

From the source (velocityComponentLaplacian_fvMotionSolver.C): at every step the solver solves on the cell centres the Laplace equation for a single component of the point velocity, cellMotionUx, with boundary conditions taken from pointMotionUx; it interpolates the result to the points and moves the x coordinate of every point by Δt·u_x. The y and z coordinates never change.

∇·(γ ∇u_x) = 0            (cells, PCG/DIC solver, tolerance 1e-8)
γ_f = n_f · diag(1, 200, 0) · n_f = n_x² + 200 n_y²   (<a href="https://github.com/OpenFOAM/OpenFOAM-14/blob/master/src/fvMeshMovers/fvMotionSolvers/motionDiffusivity/directional/directionalDiffusivity.C">directionalDiffusivity.C</a>)
x_p ← x_p + Δt u_x,p

The diffusivity depends on the face orientation: it is 200 times larger on faces normal to y than on those normal to x. Motion therefore propagates “vertically” almost without attenuation (the cell columns aligned with the cone move together with it) while compression and stretching concentrate along x, in the blocks in front and behind. This is the choice that keeps cells regular in the gap. The alternatives available in version 14 are the motionDiffusivity classes: uniform, directional, motionDirectional, inverseDistance, inverseFaceDistance, inversePointDistance, inverseVolume, file and the quadratic and exponential modifiers.

Laplacian motion solvers with variable diffusivity are described by H. Jasak, Ž. Tuković, “Automatic mesh motion for the unstructured finite volume method”, Transactions of FAMENA 30(2), 1–20, 2006 (CROSBI record), and dynamic mesh handling in OpenFOAM by H. Jasak, “Dynamic mesh handling in OpenFOAM”, AIAA 2009-341, doi:10.2514/6.2009-341.

Reading of meshToMesh_fvMeshTopoChanger.C, with the practical consequences:

• Where the meshes are. The topoChanger reads constant/meshes/<time>/polyMesh. The files in system/meshes/<time>/ are only the dictionaries: blockMesh -mesh 0.0015 reads system/meshes/0.0015/blockMeshDict and writes the mesh to constant/meshes/0.0015 (setMeshPath.H: path "meshes"/<name>). Without the two blockMesh -mesh commands the run stops at the first event.
• When it fires. With cycle the “mesh time” is t in even cycles and cycle − (t mod cycle) in odd cycles: the sequence is walked forward and then backward. A mesh is loaded when the mesh time, rounded to timeDelta = 1 µs, coincides with one of the times. In the 9.9 ms of the tutorial there are six events: 1.5 and 3.0 ms forward, 3.667 and 5.167 ms backward, 8.167 and 9.667 ms forward again (section 4).
• Time step. The topoChanger adds the meshToMeshAdjustTimeStep function object by itself, which limits the step to the time left to the next event: the run lands exactly on the remapping instants, which can be recognised in the log by suddenly shorter steps.
• Field transfer. Volume fields are transferred with fvMeshToFvMesh and the intersection method (geometric intersection of old and new cells, volume-weighted average); surface and point fields are set to NaN and rebuilt: Uf is reinterpolated from U, and the flux is recomputed and made divergence-free by correctPhi (section 1). The step after the event therefore starts from a field consistent with the new mesh, but the history of old-time values (clearOldTimes) is reset.
• Log message. At every event “Mapping to mesh time 0.0015” (or 0.003) appears, printed before the “Time =” line of the next step.

The volume-weighted average over intersections conserves the volume integrals of the transferred fields (except for cell portions not covered): it is the principle of conservative projection described, for general meshes, by P. E. Farrell, J. R. Maddison, “Conservative interpolation between volume meshes by local Galerkin projection”, Computer Methods in Applied Mechanics and Engineering 200, 89–100, 2011, doi:10.1016/j.cma.2010.07.015. The mean value is conserved but gradients are smoothed: this is why the guide asks to quantify field jumps at the critical instants.

From movingWallVelocityFvPatchVectorField.C: if the mesh is moving, the face value is

Up = (C_f(t^(n+1)) − C_f(t^n)) / Δt          face-centre velocity
Un = meshPhi_f / |S_f|                      normal velocity consistent with the swept volume
U_f = Up + n (Un − n·Up)

The tangential component is that of the wall motion (no-slip), the normal component is the one that makes the relative flux through the face zero consistently with the swept volume, i.e. with the GCL. For a stationary wall the condition reduces to zero velocity. On farFieldMoving the tutorial uses noSlip because there the physical wall (the outer surface) is still and only the mesh points slide over it: movingWallVelocity would give the fluid the translation velocity of the mesh, i.e. a sliding outer wall, which is physically wrong.

• p is kinematic pressure [m²/s²] relative to the p₀ = 0 imposed by totalPressure on left: the physical gauge pressure is ρ·p.
• The tutorial fixes only ν = 1·10⁻⁵ m²/s: the density must be chosen according to the fluid to be represented, and the choice declared in the results.
• The coneForces function of functions_guide16 uses rhoInf=1: forces come out in m⁴/s², for the 5° wedge. Real force on the axisymmetric body = ρ · 72 · F_x (the y and z components cancel by symmetry).
• The sign of the force on the cone reverses with the motion: during the approach the gap pressure brakes it, during the retreat it pulls it back. It is a good indicator to recognise in plots the direction changes at 3.333 and 6.667 ms.

The official Allrun file:

#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

runApplication blockMesh

for mapTime in 0.0015 0.003
do
    runApplication -s $mapTime blockMesh -mesh $mapTime
done

runApplication foamRun

#------------------------------------------------------------------------------

runApplication -s $mapTime writes the logs log.blockMesh.0.0015 and log.blockMesh.0.003 next to log.blockMesh and log.foamRun. In the solver log, in the order they appear at every step: “Courant Number” (computed in preSolve()), the possible “Mapping to mesh time”, “deltaT =”, “Time =”, then the solutions of cellMotionUx, pcorr (if correctPhi), U and p, the possible “Mesh Courant Number” and “ExecutionTime”. The complete sequence of events:

Between 5.167 and 8.167 ms the run stays on mesh 0.0015, which the motion solver carries backward to zero displacement (6.667 ms) and forward again: the initial mesh is never reloaded. This is the stretch in which the quality of the cells in front of the cone must be checked most carefully, because that block is the most stretched.

• correctPhi yes: recomputes a divergence-free flux after every mesh motion (section 1); with a topology change it happens anyway.
• "pcorr.*" (GAMG, tolerance 0.02): the correction equation used by correctPhi; it is deliberately loose because it only cleans the flux before the PIMPLE loop.
• cellMotionUx (PCG/DIC, tolerance 1e-8, relTol 0): the mesh motion equation is solved tightly, because an error here becomes a point-position error that accumulates.
• laplacian(diffusivity,cellMotionU) Gauss linear uncorrected: the motion Laplacian has no non-orthogonal correction.
• nOuterCorrectors 2: two PIMPLE iterations per step; the mesh moves at the first one (moveMeshOuterCorrectors default no).
• adjustTimeStep yes, maxCo 0.5: the step grows by at most deltaTFactor = 1.2 per step (default, solver.C) and shortens to hit the remapping events.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Prepare independent OpenFOAM-14 movingCone teaching cases safely."""
import argparse
import os
from pathlib import Path
import re
import shutil
import sys

VARIANTS = {
    "baseline": {},
    "courant_half": {"deltaT": "2.5e-06", "maxCo": "0.25"},
    "dense_output": {"writeInterval": "5e-05"},
}
REQUIRED = (
    "Allrun", "0/pointMotionUx", "0/U", "0/p", "constant/dynamicMeshDict",
    "system/controlDict", "system/fvSolution", "system/fvSchemes",
    "system/blockMeshDict", "system/meshes/0.0015/blockMeshDict",
    "system/meshes/0.003/blockMeshDict",
)


def update_dict(text, key, value):
    pat = re.compile(r"(?m)^(\s*" + re.escape(key) + r"\s+)([^;\n]+)(;)")
    changed, n = pat.subn(lambda m: m.group(1) + value + m.group(3), text)
    if n != 1:
        raise ValueError(f"Expected exactly one {key} entry; found {n}")
    return changed


def prepare(source, target):
    source = source.expanduser().resolve()
    target = target.expanduser().resolve()
    if not source.is_dir():
        raise ValueError(f"Source tutorial directory does not exist: {source}")
    missing = [p for p in REQUIRED if not (source / p).is_file()]
    if missing:
        raise ValueError("Reference tutorial is incomplete: " + ", ".join(missing))
    if source == target or source in target.parents:
        raise ValueError("Output directory cannot be inside the source tutorial")
    if target.exists() and not target.is_dir():
        raise ValueError(f"Target is not a directory: {target}")
    destinations = {name: target/name for name in VARIANTS}
    existing = [str(p) for p in destinations.values() if p.exists() or p.is_symlink()]
    if existing:
        raise FileExistsError("Refusing to overwrite existing case(s): " + ", ".join(existing))
    target.mkdir(parents=True, exist_ok=True)
    for name, changes in VARIANTS.items():
        dest = destinations[name]
        # copytree refuses pre-existing destinations and never replaces files.
        shutil.copytree(source, dest, symlinks=False)
        conf = dest / "system/controlDict"
        raw = conf.read_text(encoding="utf-8")
        for key, val in changes.items():
            raw = update_dict(raw, key, val)
        conf.write_text(raw, encoding="utf-8", newline="\n")
        (dest / "GUIDE16_CASE.txt").write_text(
            "OpenFOAM Foundation 14 movingCone tutorial copy\n"
            f"Variant: {name}\nOverrides: {changes or 'none'}\n"
            "The original tutorial remains untouched.\n", encoding="utf-8")
        print(f"Created: {dest}")
    print("Run cases individually using run_case.sh, after reviewing their input files.")


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    default_src = (Path(os.environ["FOAM_TUTORIALS"]) / "incompressibleFluid" / "movingCone"
                   if os.environ.get("FOAM_TUTORIALS") else None)
    ap.add_argument("--source", type=Path, default=default_src,
                    help="Path to official movingCone tutorial (default: $FOAM_TUTORIALS/...)")
    ap.add_argument("--dest", type=Path, required=True,
                    help="Parent folder for three NEW case directories")
    args = ap.parse_args()
    if args.source is None:
        ap.error("Set FOAM_TUTORIALS or provide --source")
    try:
        prepare(args.source, args.dest)
    except (OSError, ValueError) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        sys.exit(1)


if __name__ == "__main__":
    main()

Test on a copy of the version 14 tutorial (the three folders created, the refusal to overwrite and the only changes in controlDict):

$ python3 prepare_cases.py --source movingCone_tutorial --dest guide16_cases
Created: guide16_cases/baseline
Created: guide16_cases/courant_half
Created: guide16_cases/dense_output
Run cases individually using run_case.sh, after reviewing their input files.
$ python3 prepare_cases.py --source movingCone_tutorial --dest guide16_cases   # second time
ERROR: Refusing to overwrite existing case(s): guide16_cases/baseline, guide16_cases/courant_half, guide16_cases/dense_output
$ diff movingCone_tutorial/system/controlDict guide16_cases/courant_half/system/controlDict
27c27
< deltaT          5e-06;
---
> deltaT          2.5e-06;
49c49
< maxCo           0.5;
---
> maxCo           0.25;
$ cat guide16_cases/courant_half/GUIDE16_CASE.txt
OpenFOAM Foundation 14 movingCone tutorial copy
Variant: courant_half
Overrides: {'deltaT': '2.5e-06', 'maxCo': '0.25'}
The original tutorial remains untouched.

The update_dict function requires exactly one occurrence of the entry: a controlDict with a repeated deltaT (for instance inside a function object) makes the preparation fail instead of changing the wrong entry. The script also copies Allrun-parallel and the decomposeParDict files, but does not change them.

examples/experiment_manifest.csv

case_id,initial_deltaT_s,maxCo,writeInterval_s,motion_model,topology_map_times_s
baseline,5e-06,0.5,1e-4,velocityComponentLaplacian,"0.0015; 0.003"
courant_half,2.5e-06,0.25,1e-4,velocityComponentLaplacian,"0.0015; 0.003"
dense_output,5e-06,0.5,5e-5,velocityComponentLaplacian,"0.0015; 0.003"

examples/measurement_sheet.csv

case_id,source_run,mesh_valid,mesh_co_max,co_max,conservation_check,elapsed_s,notes
baseline,,,,,,,
courant_half,,,,,,,
dense_output,,,,,,,

The first describes the three variants; the second is the sheet to fill in. In the conservation_check column it is convenient to report the maximum |Q|/Q_ref of volume_balance.py (section 2), in the mesh_co_max column the value extracted with checkMeshCourantNo yes enabled, otherwise the cell stays empty.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Execute a prepared OpenFOAM movingCone case, refusing any re-run.
set -euo pipefail
if [[ $# -ne 1 ]]; then
  echo "Usage: $0 /absolute/path/to/new/case" >&2
  exit 2
fi
case_dir=$(realpath -e -- "$1")
[[ -f "$case_dir/GUIDE16_CASE.txt" && -f "$case_dir/Allrun" ]] || {
  echo "ERROR: Not a Guide 16 prepared case" >&2; exit 1;
}
command -v blockMesh >/dev/null || { echo 'ERROR: OpenFOAM environment not loaded' >&2; exit 1; }
command -v foamRun >/dev/null || { echo 'ERROR: foamRun not available' >&2; exit 1; }
# Refuse to resume/overwrite data. Never delete an earlier mesh or result.
shopt -s nullglob
outputs=("$case_dir"/log.* "$case_dir"/processor* "$case_dir"/run.full.log)
for d in "$case_dir"/*; do
  [[ -d "$d" && "${d##*/}" =~ ^[0-9]+([.][0-9]+)?$ && "${d##*/}" != 0 ]] && outputs+=("$d")
done
(( ${#outputs[@]} == 0 )) || { echo 'ERROR: Existing logs/results: use a NEW prepared case.' >&2; exit 1; }
cd -- "$case_dir"
# Allrun invokes blockMesh, two remapping meshes, then foamRun.
bash ./Allrun 2>&1 | tee run.full.log
[[ -f log.foamRun ]] || { echo 'ERROR: log.foamRun missing' >&2; exit 1; }
if grep -Eq 'FOAM FATAL|Floating point exception|Segmentation fault' log.foamRun; then
  echo 'ERROR: Solver log contains a fatal failure' >&2; exit 1
fi
if ! grep -Eq '^[[:space:]]*End[[:space:]]*$' log.foamRun; then
  echo 'WARNING: No final End marker found. Inspect solver log before accepting results.' >&2
  exit 1
fi
echo "Run finished; inspect quality and mass conservation before accepting results."

With shopt -s nullglob a glob with no match disappears, but a literal path does not: "$case_dir"/run.full.log always enters the outputs array, which is therefore never empty. The script refuses even a freshly prepared case and the section 9 command does not start. The corrected version counts run.full.log only if it exists and, in addition, refuses a case whose meshes were already generated (blockMesh would overwrite them):

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Execute a prepared OpenFOAM movingCone case, refusing any re-run.
# Corrected version of the Guide 16 package script (giux-lab.it): the original
# counted run.full.log even when absent and therefore refused every case.
# It also refuses a case whose meshes were already generated (blockMesh would
# overwrite constant/polyMesh and constant/meshes).
set -euo pipefail
if [[ $# -ne 1 ]]; then
  echo "Usage: $0 /absolute/path/to/new/case" >&2
  exit 2
fi
case_dir=$(realpath -e -- "$1")
[[ -f "$case_dir/GUIDE16_CASE.txt" && -f "$case_dir/Allrun" ]] || {
  echo "ERROR: Not a Guide 16 prepared case" >&2; exit 1;
}
command -v blockMesh >/dev/null || { echo 'ERROR: OpenFOAM environment not loaded' >&2; exit 1; }
command -v foamRun >/dev/null || { echo 'ERROR: foamRun not available' >&2; exit 1; }
# Refuse to resume/overwrite data. Never delete an earlier mesh or result.
shopt -s nullglob
# Only paths that really exist: with nullglob an unmatched glob expands to nothing,
# but a literal path would always be counted, so run.full.log is tested explicitly.
outputs=("$case_dir"/log.* "$case_dir"/processor*)
[[ -e "$case_dir/run.full.log" ]] && outputs+=("$case_dir/run.full.log")
for d in "$case_dir"/constant/meshes "$case_dir"/constant/polyMesh; do
  [[ -e "$d" ]] && outputs+=("$d")
done
for d in "$case_dir"/*; do
  [[ -d "$d" && "${d##*/}" =~ ^[0-9]+([.][0-9]+)?$ && "${d##*/}" != 0 ]] && outputs+=("$d")
done
(( ${#outputs[@]} == 0 )) || { echo 'ERROR: Existing logs/results: use a NEW prepared case.' >&2; exit 1; }
cd -- "$case_dir"
# Allrun invokes blockMesh, two remapping meshes, then foamRun.
bash ./Allrun 2>&1 | tee run.full.log
[[ -f log.foamRun ]] || { echo 'ERROR: log.foamRun missing' >&2; exit 1; }
if grep -Eq 'FOAM FATAL|Floating point exception|Segmentation fault' log.foamRun; then
  echo 'ERROR: Solver log contains a fatal failure' >&2; exit 1
fi
if ! grep -Eq '^[[:space:]]*End[[:space:]]*$' log.foamRun; then
  echo 'WARNING: No final End marker found. Inspect solver log before accepting results.' >&2
  exit 1
fi
echo "Run finished; inspect quality and mass conservation before accepting results."

Test with the OpenFOAM commands replaced by stand-ins (same sequence of files and logs as the tutorial):

$ bash run_case.sh "$PWD/guide16_cases/baseline"      # original script, fresh case
ERROR: Existing logs/results: use a NEW prepared case.
exit code 1
$ bash run_case_fixed.sh "$PWD/guide16_cases/baseline"
Run finished; inspect quality and mass conservation before accepting results.
exit code 0
$ bash run_case_fixed.sh "$PWD/guide16_cases/baseline"   # second time
ERROR: Existing logs/results: use a NEW prepared case.
exit code 1

In the downloadable archive the corrected version is extra/run_case_fixed.sh: use it instead of run_case.sh in the same commands of exercises A, B and C.

cd "$HOME/OpenFOAM/guide16_cases/baseline"
grep -n "Mapping to mesh time" log.foamRun            # six events expected
grep -c "^Time = " log.foamRun                       # number of steps
grep "^deltaT = " log.foamRun | sort -g -k3 | head -3   # shortest steps
grep "Courant Number mean" log.foamRun | sort -g -k6 | tail -3
grep -n "Solving for cellMotionUx" log.foamRun | tail -2
grep -Ein "FOAM FATAL|nan|negative|failed" log.foamRun | head
grep -n "Mesh OK\|Failed" log.foamRun | tail     # with the checkMesh function enabled

The maximum Courant number is in the sixth column of the “Courant Number mean: … max: …” line. If there are fewer than six remapping events, the run stopped before endTime or the mesh times are not the expected ones.

With an adaptive step the two cases do not write at the same instants. The script reads the probe files (format of the OpenFOAM 14 probes function object, probes.C), checks that the locations coincide, interpolates case A onto the times of case B without extrapolating and applies the guide formula, with mandatory q_scale and tolerance:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Compare one probed quantity between two OpenFOAM runs at the same physical times (Guide 16, Exercise B).

Reads two OpenFOAM probe files (postProcessing/<name>/<startTime>/<field>, written by the
'probes' function object: '# Probe i (x y z)' header, then 'time v0 v1 ...', vectors as '(a b c)'),
checks that the probe locations coincide, interpolates run A linearly onto the times of run B
(only inside the common time range, never extrapolating), and computes for every probe

    eps_q(t) = |q_A(t) - q_B(t)| / max(|q_B(t)|, q_scale)

as in the guide, with q = the scalar field or a vector component / magnitude. q_scale and the
tolerance must be chosen BEFORE looking at the result and are required arguments.
Prints max and RMS of eps per probe and the time of the maximum; optional CSV (new file only).
"""
import argparse
import bisect
import csv
import math
import re
from pathlib import Path

NUM = r'[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?'


def read_probes(path, comp):
    locs, times, vals = [], [], []
    for line in path.read_text(encoding='utf-8').splitlines():
        s = line.strip()
        if not s:
            continue
        if s.startswith('#'):
            m = re.match(r'#\s*Probe\s+(\d+)\s*\(([^)]*)\)', s)
            if m:
                locs.append(tuple(float(x) for x in m.group(2).split()))
            continue
        t = float(s.split()[0])
        rest = s[len(s.split()[0]):]
        vecs = re.findall(r'\(([^)]*)\)', rest)
        if vecs:
            row = []
            for v in vecs:
                c = [float(x) for x in v.split()]
                row.append(math.sqrt(sum(x * x for x in c)) if comp == 'mag' else c['xyz'.index(comp)])
        else:
            if comp not in ('scalar', 'mag'):
                raise SystemExit(f'{path}: scalar field, use --component scalar')
            row = [abs(float(x)) if comp == 'mag' else float(x) for x in re.findall(NUM, rest)]
        times.append(t)
        vals.append(row)
    if not times:
        raise SystemExit(f'{path}: no data rows')
    return locs, times, vals


def interp(times, vals, t, k):
    i = bisect.bisect_left(times, t)
    if i < len(times) and abs(times[i] - t) <= 1e-12 * max(1.0, abs(t)):
        return vals[i][k]
    if i == 0 or i == len(times):
        return None
    t0, t1 = times[i - 1], times[i]
    w = (t - t0) / (t1 - t0)
    return (1 - w) * vals[i - 1][k] + w * vals[i][k]


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('fileA', type=Path, help='probe file of run A (e.g. courant_half)')
    ap.add_argument('fileB', type=Path, help='probe file of run B, the reference (e.g. baseline)')
    ap.add_argument('--component', default='scalar', choices=('scalar', 'x', 'y', 'z', 'mag'))
    ap.add_argument('--q-scale', type=float, required=True, help='q_scale of the guide formula, same units as q')
    ap.add_argument('--tol', type=float, required=True, help='acceptance tolerance on max eps')
    ap.add_argument('--csv', type=Path, help='optional output CSV (new file only)')
    a = ap.parse_args()
    la, ta, va = read_probes(a.fileA, a.component)
    lb, tb, vb = read_probes(a.fileB, a.component)
    if len(la) != len(lb) or any(max(abs(p - q) for p, q in zip(x, y)) > 1e-12 for x, y in zip(la, lb)):
        raise SystemExit('probe locations differ between the two files: not comparable')
    rows, ok = [], True
    for k, loc in enumerate(lb):
        eps = []
        for t, row in zip(tb, vb):
            qa = interp(ta, va, t, k)
            if qa is None:
                continue
            qb = row[k]
            e = abs(qa - qb) / max(abs(qb), a.q_scale)
            eps.append((t, e))
            rows.append((k, t, qa, qb, e))
        if not eps:
            raise SystemExit('no common time range')
        emax = max(eps, key=lambda x: x[1])
        rms = math.sqrt(sum(e * e for _, e in eps) / len(eps))
        verdict = 'within' if emax[1] <= a.tol else 'ABOVE'
        ok &= emax[1] <= a.tol
        print(f'probe {k} {loc}: {len(eps)} times, max eps {emax[1]:.4g} at t = {emax[0]:.6g} s, RMS {rms:.4g} -> {verdict} tolerance {a.tol:g}')
    print('RESULT: ' + ('all probes within tolerance' if ok else 'at least one probe above tolerance'))
    if a.csv:
        if a.csv.exists():
            raise SystemExit(f'refusing to overwrite {a.csv}')
        with a.csv.open('x', newline='', encoding='utf-8') as f:
            w = csv.writer(f)
            w.writerow(('probe', 'time_s', 'qA_interp', 'qB', 'eps'))
            w.writerows(rows)
        print(f'wrote {a.csv}')


if __name__ == '__main__':
    main()

B=$HOME/OpenFOAM/guide16_cases
python3 compare_probes.py "$B/courant_half/postProcessing/gapProbes/0/p" "$B/baseline/postProcessing/gapProbes/0/p" \
    --q-scale 0.5 --tol 0.05
python3 compare_probes.py "$B/courant_half/postProcessing/gapProbes/0/U" "$B/baseline/postProcessing/gapProbes/0/U" \
    --component x --q-scale 0.1 --tol 0.05 --csv eps_Ux.csv

The probes exist only if enable_monitors.py was run on both copies before the computation. The tutorial time scheme is implicit Euler, first order: halving the step roughly halves the time truncation error, so the difference between two cases estimates the error of the finer one. Richardson extrapolation in time needs three steps with a constant ratio (Guide 4); with maxCo the ratio of the actual steps is not guaranteed constant, and must be read from the log.

• besides U and p: Uf (face velocity, needed for restart), pointMotionUx and cellMotionUx, and the polyMesh folder with at least points, because the geometry changes at every step;
• after a remapping the polyMesh written is complete (faces, cells, patches), because the topology changed;
• with writeInterval 1e-4 and endTime 0.0099 there are 99 directories plus 0; with 5e-05 (dense_output case) there are 198 plus 0: disk space roughly doubles, computing time hardly changes;
• the remapping instants 3.667, 5.167, 8.167 and 9.667 ms are not multiples of the write interval: to see the mesh immediately before and after these events you must write at every step in that window, in a dedicated copy.

cd "$HOME/OpenFOAM/guide16_cases/dense_output"
foamListTimes | wc -l
du -sh . postProcessing
ls 0.0015/polyMesh 0.0016/polyMesh

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Extract optional time-step diagnostics from OpenFOAM solver logs to CSV."""
import argparse
import csv
from pathlib import Path
import re

FLOAT = r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?"
TIME = re.compile(r"^Time\s*=\s*(" + FLOAT + r")\b")
COURANT = re.compile(r"^Courant Number mean:\s*" + FLOAT + r"\s+max:\s*(" + FLOAT + r")")
MESH = re.compile(r"^Mesh Courant Number mean:\s*" + FLOAT + r"\s+max:\s*(" + FLOAT + r")")
EXEC = re.compile(r"^ExecutionTime\s*=\s*(" + FLOAT + r")\s+s")


def extract(lines):
    records = []
    current = None
    for line in lines:
        line = line.strip()
        m = TIME.search(line)
        if m:
            current = {"time_s": float(m.group(1)), "max_Co": "", "max_mesh_Co": "", "execution_s": ""}
            records.append(current)
        if current is None:
            continue
        for pattern, key in [(COURANT, "max_Co"), (MESH, "max_mesh_Co"), (EXEC, "execution_s")]:
            m = pattern.search(line)
            if m:
                current[key] = float(m.group(1))
    return records


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("input", type=Path)
    ap.add_argument("output", type=Path)
    args = ap.parse_args()
    if not args.input.is_file():
        ap.error(f"Missing solver log: {args.input}")
    if args.output.exists():
        ap.error(f"Refusing overwrite: {args.output}")
    records = extract(args.input.read_text(encoding="utf-8", errors="replace").splitlines())
    if not records:
        ap.error("No Time = ... records found")
    with args.output.open("x", newline="", encoding="utf-8") as f:
        w = csv.DictWriter(f, fieldnames=("time_s", "max_Co", "max_mesh_Co", "execution_s"))
        w.writeheader()
        w.writerows(records)
    print(f"Wrote {len(records)} rows to {args.output}")

if __name__ == "__main__":
    main()

Two facts of the code change how the CSV is read. First: in foamRun.C the “Courant Number” line is printed in preSolve(), before “Time =”; it describes the step just completed and the script, which attaches lines to the previous “Time =”, correctly assigns it to the instant at which that step ended. Second: “Mesh Courant Number” appears only with checkMeshCourantNo yes (section 2), and the tutorial does not enable it: the max_mesh_Co column stays empty, exactly as the guide says not to read it as zero. The package’s synthetic_log.txt, synthetic data for testing the parser, has the lines in a different order (Courant after “Time =”):

Note: SYNTHETIC teaching-only input, not actual CFD output.
Time = 0.00001
Courant Number mean: 0.03 max: 0.15
Mesh Courant Number mean: 0.01 max: 0.04
ExecutionTime = 0.15 s  ClockTime = 0 s
Time = 0.00002
Courant Number mean: 0.05 max: 0.22
Mesh Courant Number mean: 0.02 max: 0.07
ExecutionTime = 0.32 s  ClockTime = 0 s
Time = 0.00003
Courant Number mean: 0.07 max: 0.28
ExecutionTime = 0.50 s  ClockTime = 0 s
End

$ python3 extract_log_metrics.py examples/synthetic_log.txt syn.csv
Wrote 3 rows to syn.csv
$ cat syn.csv
time_s,max_Co,max_mesh_Co,execution_s
1e-05,0.15,0.04,0.15
2e-05,0.22,0.07,0.32
3e-05,0.28,,0.5

Extended version, which adds the step used (deltaT =), the “Mapping to mesh time” events attributed to the instant at which they happen (in preSolve(), before time advances), the mean Courant number and the ClockTime, and prints a summary. Missing quantities stay empty:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Per-step diagnostics of an OpenFOAM 14 foamRun log with moving and mapped meshes, to CSV.

Extends the package's extract_log_metrics.py. For every 'Time = t' step it records:
  deltaT            the 'deltaT = ...' line printed just BEFORE 'Time = t' (step size used to reach t);
  Co_mean, Co_max   the 'Courant Number' line printed before 'Time = t': foamRun computes it in
                    preSolve() from the flux and time step of the step that ENDED at the previous time,
                    so it is stored on the previous row (prevCo_* columns keep the raw association);
  meshCo_mean/max   'Mesh Courant Number' (printed only with PIMPLE checkMeshCourantNo yes);
  mapped            1 on the row of the time at which 'Mapping to mesh time X' was printed: the mapping is
                    done in preSolve(), before the time is advanced, so it happens at the time of the
                    previous row (the instant the step size was adjusted to reach); X in map_mesh;
  exec_s, clock_s   ExecutionTime and ClockTime.
Missing quantities stay empty: nothing is reconstructed. A summary with the mapping events,
the extreme values and the number of steps is printed. The output file must not exist.
"""
import argparse
import csv
import re
from pathlib import Path

F = r'[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?'
RX = {
    'time': re.compile(r'^Time\s*=\s*(' + F + r')\s*$'),
    'dt': re.compile(r'^deltaT\s*=\s*(' + F + r')\s*$'),
    'co': re.compile(r'^Courant Number mean:\s*(' + F + r')\s+max:\s*(' + F + r')'),
    'mco': re.compile(r'^Mesh Courant Number mean:\s*(' + F + r')\s+max:\s*(' + F + r')'),
    'map': re.compile(r'^Mapping to mesh time\s+(\S+)'),
    'exe': re.compile(r'^ExecutionTime\s*=\s*(' + F + r')\s*s(?:\s+ClockTime\s*=\s*(' + F + r')\s*s)?'),
}
COLS = ('time_s', 'deltaT_s', 'Co_mean', 'Co_max', 'meshCo_mean', 'meshCo_max', 'mapped', 'map_mesh', 'exec_s', 'clock_s', 'prevCo_mean', 'prevCo_max')


def parse(lines):
    rows, cur = [], None
    pend = {}                       # quantities printed before the next 'Time =' line
    for raw in lines:
        s = raw.strip()
        m = RX['co'].match(s)
        if m:
            pend['co'] = (float(m.group(1)), float(m.group(2)))
            if rows:                # belongs to the step that ended at the previous time
                rows[-1]['Co_mean'], rows[-1]['Co_max'] = pend['co']
            continue
        m = RX['map'].match(s)
        if m:
            if rows:
                rows[-1]['mapped'], rows[-1]['map_mesh'] = 1, m.group(1)
            else:
                print(f'note: mapping to mesh {m.group(1)} before the first time step')
            continue
        m = RX['dt'].match(s)
        if m:
            pend['dt'] = float(m.group(1))
            continue
        m = RX['time'].match(s)
        if m:
            cur = dict.fromkeys(COLS, '')
            cur['time_s'] = float(m.group(1))
            cur['deltaT_s'] = pend.get('dt', '')
            if 'co' in pend:
                cur['prevCo_mean'], cur['prevCo_max'] = pend['co']
            cur['mapped'] = 0
            rows.append(cur)
            pend = {}
            continue
        if cur is None:
            continue
        m = RX['mco'].match(s)
        if m:
            cur['meshCo_mean'], cur['meshCo_max'] = float(m.group(1)), float(m.group(2))
            continue
        m = RX['exe'].match(s)
        if m:
            cur['exec_s'] = float(m.group(1))
            cur['clock_s'] = float(m.group(2)) if m.group(2) else ''
    return rows


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('log', type=Path)
    ap.add_argument('output', type=Path)
    a = ap.parse_args()
    if not a.log.is_file():
        ap.error(f'missing log: {a.log}')
    if a.output.exists():
        ap.error(f'refusing to overwrite {a.output}')
    rows = parse(a.log.read_text(encoding='utf-8', errors='replace').splitlines())
    if not rows:
        ap.error('no "Time = ..." lines found')
    with a.output.open('x', newline='', encoding='utf-8') as f:
        w = csv.DictWriter(f, fieldnames=COLS)
        w.writeheader()
        w.writerows(rows)
    print(f'{len(rows)} steps written to {a.output}')
    num = lambda k: [r[k] for r in rows if r[k] != '']
    for k in ('deltaT_s', 'Co_max', 'meshCo_max'):
        v = num(k)
        print(f'  {k:11s}: ' + (f'min {min(v):.4g}  max {max(v):.4g}  ({len(v)} values)' if v else 'not present in the log'))
    maps = [(r['time_s'], r['map_mesh']) for r in rows if r['mapped'] == 1]
    print(f'  mapping events: {len(maps)}' + ''.join(f'\n    t = {t:.6g} s -> mesh {m}' for t, m in maps))
    ex = num('exec_s')
    if ex:
        print(f'  last ExecutionTime: {ex[-1]:.4g} s')


if __name__ == '__main__':
    main()

python3 extract_log_metrics_plus.py "$HOME/OpenFOAM/guide16_cases/baseline/log.foamRun" \
    "$HOME/OpenFOAM/guide16_cases/baseline/diagnostics_plus.csv"

On a real tutorial log the summary must report six remapping events, at the instants of the section 7 table. The shortest steps of the run fall just before the events (the step is shortened to hit them) and in the stretch of closest approach, where the exit velocity from the gap raises the fluid Courant number.

• The ParaView OpenFOAM reader rereads the polyMesh of every time: the cell count goes from 1690 to 1580 at 1.5 ms and to 1675 at 3.0 ms, and so on. Filters that depend on cell indices (for instance a cell selection) do not survive an event.
• The case is a 5° wedge: for a three-dimensional view use the Rotational Extrusion filter (Resolution 72, Angle 360) after extracting the boundary, or work on the z = 0 plane, which the tutorial already exports as postProcessing/cutPlane/<time>/zNormal.vtk.
• To see the mesh motion without fields: Surface With Edges representation, solid colour, and animation with “Snap To TimeSteps”. The six events can be recognised by the jump in the distribution of cells in front of and behind the cone.
• The probe and flow-rate plots (postProcessing/gapProbes, postProcessing/Qleft) can be opened directly in ParaView as text tables, or downloaded with FileZilla and read with the guide’s scripts.

The classical reference is A. H. Shapiro, M. Y. Jaffrin, S. L. Weinberg, “Peristaltic pumping with long wavelengths at low Reynolds number”, Journal of Fluid Mechanics 37, 799–825, 1969, doi:10.1017/S0022112069000899, with the review M. Y. Jaffrin, A. H. Shapiro, “Peristaltic pumping”, Annual Review of Fluid Mechanics 3, 13–37, 1971, doi:10.1146/annurev.fl.03.010171.000305. In the frame moving with the wave the flow is steady and, under the long-wave and low-Reynolds assumptions, lubrication gives for the whole gap H (one or both walls deforming vertically):

dp/dx = −12 μ (q + c H) / H³                      q = flow rate in the wave frame
Δp_λ = −12 μ λ (q ⟨H⁻³⟩ + c ⟨H⁻²⟩)               pressure rise over one wavelength
Q   = q + c ⟨H⟩                                   time-mean flow in the laboratory frame
H = H₀(1 + φ sin):  Θ = Q/(c H₀) = 3φ²/(2 + φ²)   at Δp = 0

Two consequences useful for design. The pump curve is linear: Q decreases linearly with back-pressure, from Θ at Δp = 0 to zero at the “shut-off” pressure rise. And, in the lubrication limit, only the gap H(x) matters: two opposed walls with amplitudes a and b are equivalent to a single wall with amplitude a + b. The following script computes these quantities for a sinusoidal wave, checking the numerical integration against the closed form, or for the CSV profiles of peristaltic_gap_demo.py, treated as one wavelength of a periodic wave:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Long-wavelength, low-Reynolds peristaltic pumping of a 2-D channel (lubrication theory, Shapiro et al. 1969).

A wave of gap H(x - c t) travels at speed c. In the wave frame the flow is steady, both walls move at -c,
and lubrication theory gives for the full gap (one or both walls deforming vertically):
    dp/dx = -12 mu (q + c H) / H^3,     q = wave-frame flow rate per unit depth (constant).
Over one wavelength lambda the pressure rise is
    Delta p = -12 mu [ q I3 + c I2 ] lambda,   I2 = <H^-2>, I3 = <H^-3>  (averages over one wavelength),
and the time-mean flow in the laboratory frame is Q = q + c <H>. Hence the pump is linear:
    Q = c <H> - c I2/I3 - Delta p / (12 mu lambda I3).
For a sinusoidal gap H = H0 (1 + phi sin) the zero-pressure flow is Theta = Q/(c H0) = 3 phi^2 / (2 + phi^2),
the classical result; the script checks its numerical integration against this formula.
Input: either --sine PHI, or a CSV of the package's peristaltic_gap_demo.py (columns mode,time_s,x_m,gap_m),
treated as ONE wavelength of a periodic wave (lambda = x range). Output: Theta, the pressure rise per wavelength at zero
net flow (12 (<H> I3 - I2) mu c lambda) and, with --mu and --c, dimensional values. Planar 2-D theory: no tube wall mechanics, no inertia,
no leakage at contact; it is a design-stage estimate, not a CFD result. Read-only.
"""
import argparse
import csv
import math
from pathlib import Path


def averages(H):
    n = len(H)
    return sum(H) / n, sum(h ** -2 for h in H) / n, sum(h ** -3 for h in H) / n


def from_csv(path, time):
    rows = [r for r in csv.DictReader(path.open(encoding='utf-8')) if abs(float(r['time_s']) - time) < 1e-12]
    if not rows:
        raise SystemExit(f'no rows at time_s = {time}')
    rows.sort(key=lambda r: float(r['x_m']))
    x = [float(r['x_m']) for r in rows]
    H = [float(r['gap_m']) for r in rows]
    # periodic: drop the last point if it repeats the first position of the next wavelength
    return x[-1] - x[0], H[:-1] if len(H) > 2 else H, rows[0]['mode']


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    g = ap.add_mutually_exclusive_group(required=True)
    g.add_argument('--sine', type=float, metavar='PHI', help='sinusoidal gap with amplitude ratio 0 <= phi < 1')
    g.add_argument('--csv', type=Path, help='CSV from peristaltic_gap_demo.py')
    ap.add_argument('--time', type=float, default=0.0, help='time_s of the CSV profile to use')
    ap.add_argument('--n', type=int, default=20000, help='samples per wavelength for --sine')
    ap.add_argument('--mu', type=float, help='dynamic viscosity [Pa s] for dimensional output')
    ap.add_argument('--c', type=float, help='wave speed [m/s] for dimensional output')
    a = ap.parse_args()
    if a.sine is not None:
        if not 0 <= a.sine < 1:
            ap.error('phi must be in [0, 1)')
        lam, H0 = 1.0, 1.0
        H = [H0 * (1 + a.sine * math.sin(2 * math.pi * k / a.n)) for k in range(a.n)]
        label = f'sine phi = {a.sine}'
    else:
        lam, H, mode = from_csv(a.csv, a.time)
        label = f'{a.csv.name} ({mode}, t = {a.time} s), lambda = {lam:.4g} m'
    Hm, I2, I3 = averages(H)
    theta0 = (Hm - I2 / I3) / Hm               # Q/(c <H>) at zero pressure rise
    print(label)
    print(f'<H> = {Hm:.6g}, min H = {min(H):.6g}, occlusion 1 - Hmin/<H> = {1 - min(H)/Hm:.4f}')
    print(f'Theta = Q/(c <H>) at zero pressure rise = {theta0:.6f}')
    # pressure rise per wavelength at zero time-mean flow: Q = 0 -> q = -c <H>
    dp0_nd = 12 * (Hm * I3 - I2) * Hm ** 2        # multiply by mu c lambda / <H>^2
    if a.sine is not None:
        exact = 3 * a.sine ** 2 / (2 + a.sine ** 2)
        print(f'closed form 3 phi^2/(2+phi^2)       = {exact:.6f}   (difference {abs(theta0-exact):.2e})')
    if a.mu and a.c:
        dp_block = 12 * a.mu * a.c * lam * (Hm * I3 - I2)
        print(f'flow per unit depth at zero Delta p: Q = {theta0*a.c*Hm:.4g} m2/s')
        print(f'pressure rise per wavelength at zero net flow (blocked outlet): {dp_block:.4g} Pa')
    else:
        print(f'pressure rise per wavelength at zero net flow = {dp0_nd:.4g} x mu c lambda / <H>^2')


if __name__ == '__main__':
    main()

$ python3 peristaltic_lubrication.py --sine 0.4
sine phi = 0.4
<H> = 1, min H = 0.6, occlusion 1 - Hmin/<H> = 0.4000
Theta = Q/(c <H>) at zero pressure rise = 0.222222
closed form 3 phi^2/(2+phi^2)       = 0.222222   (difference 2.78e-17)
pressure rise per wavelength at zero net flow = 4.453 x mu c lambda / <H>^2
$ python3 peristaltic_lubrication.py --sine 0.7
sine phi = 0.7
<H> = 1, min H = 0.3, occlusion 1 - Hmin/<H> = 0.7000
Theta = Q/(c <H>) at zero pressure rise = 0.590361
closed form 3 phi^2/(2+phi^2)       = 0.590361   (difference 1.11e-16)
pressure rise per wavelength at zero net flow = 47.48 x mu c lambda / <H>^2
$ python3 peristaltic_lubrication.py --csv examples/upper_only_demo.csv --mu 1e-3 --c 0.03
upper_only_demo.csv (upper_only, t = 0.0 s), lambda = 0.1 m
<H> = 0.00112911, min H = 0.0007, occlusion 1 - Hmin/<H> = 0.3800
Theta = Q/(c <H>) at zero pressure rise = 0.083046
flow per unit depth at zero Delta p: Q = 2.813e-06 m2/s
pressure rise per wavelength at zero net flow (blocked outlet): 2.747 Pa
$ python3 peristaltic_lubrication.py --csv examples/opposed_demo.csv --mu 1e-3 --c 0.03
opposed_demo.csv (opposed, t = 0.0 s), lambda = 0.1 m
<H> = 0.00107949, min H = 0.00035, occlusion 1 - Hmin/<H> = 0.6758
Theta = Q/(c <H>) at zero pressure rise = 0.415427
flow per unit depth at zero Delta p: Q = 1.345e-05 m2/s
pressure rise per wavelength at zero net flow (blocked outlet): 32.39 Pa

With the package profiles the opposed configuration, at equal wave speed, gives about five times more flow because the minimum gap drops from 0.70 to 0.35 mm (occlusion from 0.38 to 0.68). The values are lubrication estimates on a planar channel: they do not include tube mechanics, contact, leakage at closure or inertia (with c = 0.03 m/s and H ≈ 1 mm the Reynolds number for water is of order 30), and they serve to set up the CFD campaign, not to replace it.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Generate purely geometric, non-CFD traveling-occlusion profiles as CSV."""
import argparse
import csv
from pathlib import Path
import math


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("output", type=Path)
    ap.add_argument("--mode", choices=("upper_only", "opposed"), default="upper_only")
    args = ap.parse_args()
    x0, length, h0, speed, width = 0.0, 0.10, 0.0012, 0.03, 0.008
    a_top, a_bottom = 0.00050, (0.00035 if args.mode == "opposed" else 0.0)
    rows=[]
    for t in (0.0, 1.0, 2.0):
        center = 0.02 + speed*t
        for i in range(101):
            x = x0 + length*i/100.0
            w = math.exp(-((x-center)/width)**2)
            gap = h0-(a_top+a_bottom)*w
            assert gap>0.0, "Demo must not imply full tube occlusion"
            rows.append((args.mode,t,x,gap))
    with args.output.open("x",newline="",encoding="utf-8") as f:
        wr=csv.writer(f)
        wr.writerow(("mode","time_s","x_m","gap_m"))
        wr.writerows(rows)
    print("Geometric-only CSV written. No fluid or wall mechanics solved.")

if __name__ == "__main__":
    main()

Fixed parameters in the code: length 0.10 m, rest gap h₀ = 1.2 mm, Gaussian width w = 8 mm, speed c = 0.03 m/s, times 0, 1 and 2 s (centre at 0.02, 0.05 and 0.08 m), amplitudes 0.50 mm (upper wall) and 0.35 mm (lower, only in opposed). The assertion prevents complete occlusion. Regenerated with the two section 14 commands, the files are byte-for-byte identical to those of the package. The full content of the two CSV files (geometry only, no simulation):

examples/upper_only_demo.csv

mode,time_s,x_m,gap_m
upper_only,0.0,0.0,0.0011990347729318861
upper_only,0.0,0.001,0.0011982246757213787
upper_only,0.0,0.002,0.001196835142286257
upper_only,0.0,0.0030000000000000005,0.0011945316162446974
upper_only,0.0,0.004,0.0011908421805556328
upper_only,0.0,0.005,0.0011851353918069206
upper_only,0.0,0.006000000000000001,0.0011766146888080205
upper_only,0.0,0.007000000000000001,0.0011643416586511208
upper_only,0.0,0.008,0.0011473003877190678
upper_only,0.0,0.009000000000000001,0.0011245112907720426
upper_only,0.0,0.01,0.001095194306424451
upper_only,0.0,0.011000000000000001,0.0010589685241530922
upper_only,0.0,0.012000000000000002,0.0010160602794142787
upper_only,0.0,0.013000000000000001,0.0009674784059329716
upper_only,0.0,0.014000000000000002,0.0009151085876345382
upper_only,0.0,0.015,0.0008616830769191355
upper_only,0.0,0.016,0.0008105996084642975
upper_only,0.0,0.017,0.0007655924718685782
upper_only,0.0,0.018000000000000002,0.000730293468593262
upper_only,0.0,0.019000000000000003,0.0007077517814972956
upper_only,0.0,0.02,0.0006999999999999999
upper_only,0.0,0.021,0.0007077517814972958
upper_only,0.0,0.022000000000000002,0.000730293468593262
upper_only,0.0,0.023000000000000003,0.0007655924718685783
upper_only,0.0,0.024000000000000004,0.0008105996084642976
upper_only,0.0,0.025,0.0008616830769191355
upper_only,0.0,0.026000000000000002,0.0009151085876345384
upper_only,0.0,0.027000000000000003,0.0009674784059329718
upper_only,0.0,0.028000000000000004,0.0010160602794142789
upper_only,0.0,0.029000000000000005,0.0010589685241530922
upper_only,0.0,0.03,0.001095194306424451
upper_only,0.0,0.031,0.0011245112907720426
upper_only,0.0,0.032,0.0011473003877190678
upper_only,0.0,0.033,0.0011643416586511208
upper_only,0.0,0.034,0.0011766146888080205
upper_only,0.0,0.035,0.0011851353918069206
upper_only,0.0,0.036000000000000004,0.0011908421805556328
upper_only,0.0,0.037000000000000005,0.0011945316162446974
upper_only,0.0,0.038000000000000006,0.001196835142286257
upper_only,0.0,0.03900000000000001,0.0011982246757213787
upper_only,0.0,0.04,0.0011990347729318861
upper_only,0.0,0.04100000000000001,0.0011994913610781925
upper_only,0.0,0.042,0.0011997402126589225
upper_only,0.0,0.043,0.0011998713959405995
upper_only,0.0,0.044000000000000004,0.0011999382950979565
upper_only,0.0,0.045,0.0011999713045556303
upper_only,0.0,0.046000000000000006,0.0011999870659498886
upper_only,0.0,0.047,0.0011999943495319782
upper_only,0.0,0.04800000000000001,0.0011999976074413038
upper_only,0.0,0.049,0.0011999990180958894
upper_only,0.0,0.05,0.0011999996094255296
upper_only,0.0,0.051000000000000004,0.0011999998494201275
upper_only,0.0,0.052000000000000005,0.0011999999437324125
upper_only,0.0,0.053000000000000005,0.0011999999796212303
upper_only,0.0,0.054000000000000006,0.0011999999928463789
upper_only,0.0,0.055,0.0011999999975661028
upper_only,0.0,0.05600000000000001,0.0011999999991973858
upper_only,0.0,0.057,0.001199999999743469
upper_only,0.0,0.05800000000000001,0.0011999999999205303
upper_only,0.0,0.059000000000000004,0.0011999999999761388
upper_only,0.0,0.06,0.001199999999993056
upper_only,0.0,0.061000000000000006,0.0011999999999980412
upper_only,0.0,0.062,0.0011999999999994645
upper_only,0.0,0.063,0.001199999999999858
upper_only,0.0,0.064,0.0011999999999999635
upper_only,0.0,0.065,0.0011999999999999908
upper_only,0.0,0.066,0.0011999999999999977
upper_only,0.0,0.067,0.0011999999999999995
upper_only,0.0,0.068,0.0011999999999999997
upper_only,0.0,0.069,0.0012
upper_only,0.0,0.07,0.0012
upper_only,0.0,0.07100000000000001,0.0012
upper_only,0.0,0.07200000000000001,0.0012
upper_only,0.0,0.07300000000000001,0.0012
upper_only,0.0,0.07400000000000001,0.0012
upper_only,0.0,0.075,0.0012
upper_only,0.0,0.07600000000000001,0.0012
upper_only,0.0,0.077,0.0012
upper_only,0.0,0.07800000000000001,0.0012
upper_only,0.0,0.079,0.0012
upper_only,0.0,0.08,0.0012
upper_only,0.0,0.081,0.0012
upper_only,0.0,0.08200000000000002,0.0012
upper_only,0.0,0.083,0.0012
upper_only,0.0,0.084,0.0012
upper_only,0.0,0.085,0.0012
upper_only,0.0,0.086,0.0012
upper_only,0.0,0.08700000000000001,0.0012
upper_only,0.0,0.08800000000000001,0.0012
upper_only,0.0,0.08900000000000001,0.0012
upper_only,0.0,0.09,0.0012
upper_only,0.0,0.091,0.0012
upper_only,0.0,0.09200000000000001,0.0012
upper_only,0.0,0.09300000000000001,0.0012
upper_only,0.0,0.094,0.0012
upper_only,0.0,0.095,0.0012
upper_only,0.0,0.09600000000000002,0.0012
upper_only,0.0,0.09700000000000002,0.0012
upper_only,0.0,0.098,0.0012
upper_only,0.0,0.099,0.0012
upper_only,0.0,0.1,0.0012
upper_only,1.0,0.0,0.0012
upper_only,1.0,0.001,0.0012
upper_only,1.0,0.002,0.0011999999999999997
upper_only,1.0,0.0030000000000000005,0.0011999999999999995
upper_only,1.0,0.004,0.0011999999999999977
upper_only,1.0,0.005,0.0011999999999999908
upper_only,1.0,0.006000000000000001,0.0011999999999999635
upper_only,1.0,0.007000000000000001,0.001199999999999858
upper_only,1.0,0.008,0.0011999999999994645
upper_only,1.0,0.009000000000000001,0.0011999999999980412
upper_only,1.0,0.01,0.001199999999993056
upper_only,1.0,0.011000000000000001,0.0011999999999761388
upper_only,1.0,0.012000000000000002,0.0011999999999205303
upper_only,1.0,0.013000000000000001,0.001199999999743469
upper_only,1.0,0.014000000000000002,0.0011999999991973858
upper_only,1.0,0.015,0.0011999999975661028
upper_only,1.0,0.016,0.0011999999928463789
upper_only,1.0,0.017,0.0011999999796212303
upper_only,1.0,0.018000000000000002,0.0011999999437324125
upper_only,1.0,0.019000000000000003,0.0011999998494201275
upper_only,1.0,0.02,0.0011999996094255296
upper_only,1.0,0.021,0.0011999990180958894
upper_only,1.0,0.022000000000000002,0.0011999976074413038
upper_only,1.0,0.023000000000000003,0.0011999943495319782
upper_only,1.0,0.024000000000000004,0.0011999870659498886
upper_only,1.0,0.025,0.0011999713045556303
upper_only,1.0,0.026000000000000002,0.0011999382950979565
upper_only,1.0,0.027000000000000003,0.0011998713959405995
upper_only,1.0,0.028000000000000004,0.0011997402126589225
upper_only,1.0,0.029000000000000005,0.0011994913610781925
upper_only,1.0,0.03,0.0011990347729318861
upper_only,1.0,0.031,0.0011982246757213787
upper_only,1.0,0.032,0.001196835142286257
upper_only,1.0,0.033,0.0011945316162446974
upper_only,1.0,0.034,0.0011908421805556328
upper_only,1.0,0.035,0.0011851353918069206
upper_only,1.0,0.036000000000000004,0.0011766146888080205
upper_only,1.0,0.037000000000000005,0.0011643416586511208
upper_only,1.0,0.038000000000000006,0.0011473003877190678
upper_only,1.0,0.03900000000000001,0.0011245112907720424
upper_only,1.0,0.04,0.001095194306424451
upper_only,1.0,0.04100000000000001,0.001058968524153092
upper_only,1.0,0.042,0.0010160602794142787
upper_only,1.0,0.043,0.000967478405932972
upper_only,1.0,0.044000000000000004,0.0009151085876345382
upper_only,1.0,0.045,0.0008616830769191356
upper_only,1.0,0.046000000000000006,0.0008105996084642973
upper_only,1.0,0.047,0.0007655924718685783
upper_only,1.0,0.04800000000000001,0.0007302934685932618
upper_only,1.0,0.049,0.0007077517814972958
upper_only,1.0,0.05,0.0006999999999999999
upper_only,1.0,0.051000000000000004,0.0007077517814972958
upper_only,1.0,0.052000000000000005,0.000730293468593262
upper_only,1.0,0.053000000000000005,0.0007655924718685783
upper_only,1.0,0.054000000000000006,0.0008105996084642976
upper_only,1.0,0.055,0.0008616830769191354
upper_only,1.0,0.05600000000000001,0.0009151085876345387
upper_only,1.0,0.057,0.0009674784059329716
upper_only,1.0,0.05800000000000001,0.001016060279414279
upper_only,1.0,0.059000000000000004,0.0010589685241530922
upper_only,1.0,0.06,0.0010951943064244508
upper_only,1.0,0.061000000000000006,0.0011245112907720426
upper_only,1.0,0.062,0.0011473003877190678
upper_only,1.0,0.063,0.0011643416586511208
upper_only,1.0,0.064,0.0011766146888080205
upper_only,1.0,0.065,0.0011851353918069206
upper_only,1.0,0.066,0.0011908421805556328
upper_only,1.0,0.067,0.0011945316162446974
upper_only,1.0,0.068,0.001196835142286257
upper_only,1.0,0.069,0.0011982246757213787
upper_only,1.0,0.07,0.0011990347729318861
upper_only,1.0,0.07100000000000001,0.0011994913610781925
upper_only,1.0,0.07200000000000001,0.0011997402126589225
upper_only,1.0,0.07300000000000001,0.0011998713959405995
upper_only,1.0,0.07400000000000001,0.0011999382950979565
upper_only,1.0,0.075,0.0011999713045556303
upper_only,1.0,0.07600000000000001,0.0011999870659498886
upper_only,1.0,0.077,0.0011999943495319782
upper_only,1.0,0.07800000000000001,0.0011999976074413038
upper_only,1.0,0.079,0.0011999990180958894
upper_only,1.0,0.08,0.0011999996094255296
upper_only,1.0,0.081,0.0011999998494201275
upper_only,1.0,0.08200000000000002,0.0011999999437324125
upper_only,1.0,0.083,0.0011999999796212303
upper_only,1.0,0.084,0.0011999999928463789
upper_only,1.0,0.085,0.0011999999975661028
upper_only,1.0,0.086,0.0011999999991973858
upper_only,1.0,0.08700000000000001,0.001199999999743469
upper_only,1.0,0.08800000000000001,0.0011999999999205303
upper_only,1.0,0.08900000000000001,0.0011999999999761388
upper_only,1.0,0.09,0.001199999999993056
upper_only,1.0,0.091,0.0011999999999980412
upper_only,1.0,0.09200000000000001,0.0011999999999994645
upper_only,1.0,0.09300000000000001,0.001199999999999858
upper_only,1.0,0.094,0.0011999999999999635
upper_only,1.0,0.095,0.0011999999999999908
upper_only,1.0,0.09600000000000002,0.0011999999999999977
upper_only,1.0,0.09700000000000002,0.0011999999999999995
upper_only,1.0,0.098,0.0011999999999999997
upper_only,1.0,0.099,0.0012
upper_only,1.0,0.1,0.0012
upper_only,2.0,0.0,0.0012
upper_only,2.0,0.001,0.0012
upper_only,2.0,0.002,0.0012
upper_only,2.0,0.0030000000000000005,0.0012
upper_only,2.0,0.004,0.0012
upper_only,2.0,0.005,0.0012
upper_only,2.0,0.006000000000000001,0.0012
upper_only,2.0,0.007000000000000001,0.0012
upper_only,2.0,0.008,0.0012
upper_only,2.0,0.009000000000000001,0.0012
upper_only,2.0,0.01,0.0012
upper_only,2.0,0.011000000000000001,0.0012
upper_only,2.0,0.012000000000000002,0.0012
upper_only,2.0,0.013000000000000001,0.0012
upper_only,2.0,0.014000000000000002,0.0012
upper_only,2.0,0.015,0.0012
upper_only,2.0,0.016,0.0012
upper_only,2.0,0.017,0.0012
upper_only,2.0,0.018000000000000002,0.0012
upper_only,2.0,0.019000000000000003,0.0012
upper_only,2.0,0.02,0.0012
upper_only,2.0,0.021,0.0012
upper_only,2.0,0.022000000000000002,0.0012
upper_only,2.0,0.023000000000000003,0.0012
upper_only,2.0,0.024000000000000004,0.0012
upper_only,2.0,0.025,0.0012
upper_only,2.0,0.026000000000000002,0.0012
upper_only,2.0,0.027000000000000003,0.0012
upper_only,2.0,0.028000000000000004,0.0012
upper_only,2.0,0.029000000000000005,0.0012
upper_only,2.0,0.03,0.0012
upper_only,2.0,0.031,0.0012
upper_only,2.0,0.032,0.0011999999999999997
upper_only,2.0,0.033,0.0011999999999999995
upper_only,2.0,0.034,0.0011999999999999977
upper_only,2.0,0.035,0.0011999999999999908
upper_only,2.0,0.036000000000000004,0.0011999999999999635
upper_only,2.0,0.037000000000000005,0.001199999999999858
upper_only,2.0,0.038000000000000006,0.0011999999999994645
upper_only,2.0,0.03900000000000001,0.0011999999999980412
upper_only,2.0,0.04,0.001199999999993056
upper_only,2.0,0.04100000000000001,0.0011999999999761388
upper_only,2.0,0.042,0.0011999999999205303
upper_only,2.0,0.043,0.001199999999743469
upper_only,2.0,0.044000000000000004,0.0011999999991973858
upper_only,2.0,0.045,0.0011999999975661028
upper_only,2.0,0.046000000000000006,0.0011999999928463789
upper_only,2.0,0.047,0.0011999999796212303
upper_only,2.0,0.04800000000000001,0.0011999999437324125
upper_only,2.0,0.049,0.0011999998494201275
upper_only,2.0,0.05,0.0011999996094255296
upper_only,2.0,0.051000000000000004,0.0011999990180958894
upper_only,2.0,0.052000000000000005,0.0011999976074413038
upper_only,2.0,0.053000000000000005,0.0011999943495319782
upper_only,2.0,0.054000000000000006,0.0011999870659498886
upper_only,2.0,0.055,0.0011999713045556303
upper_only,2.0,0.05600000000000001,0.0011999382950979565
upper_only,2.0,0.057,0.0011998713959405995
upper_only,2.0,0.05800000000000001,0.0011997402126589225
upper_only,2.0,0.059000000000000004,0.0011994913610781925
upper_only,2.0,0.06,0.0011990347729318861
upper_only,2.0,0.061000000000000006,0.0011982246757213787
upper_only,2.0,0.062,0.001196835142286257
upper_only,2.0,0.063,0.0011945316162446974
upper_only,2.0,0.064,0.0011908421805556328
upper_only,2.0,0.065,0.0011851353918069206
upper_only,2.0,0.066,0.0011766146888080205
upper_only,2.0,0.067,0.0011643416586511208
upper_only,2.0,0.068,0.0011473003877190678
upper_only,2.0,0.069,0.0011245112907720424
upper_only,2.0,0.07,0.0010951943064244508
upper_only,2.0,0.07100000000000001,0.001058968524153092
upper_only,2.0,0.07200000000000001,0.0010160602794142784
upper_only,2.0,0.07300000000000001,0.0009674784059329713
upper_only,2.0,0.07400000000000001,0.0009151085876345379
upper_only,2.0,0.075,0.0008616830769191356
upper_only,2.0,0.07600000000000001,0.0008105996084642969
upper_only,2.0,0.077,0.0007655924718685783
upper_only,2.0,0.07800000000000001,0.0007302934685932617
upper_only,2.0,0.079,0.0007077517814972958
upper_only,2.0,0.08,0.0006999999999999999
upper_only,2.0,0.081,0.0007077517814972958
upper_only,2.0,0.08200000000000002,0.0007302934685932624
upper_only,2.0,0.083,0.0007655924718685783
upper_only,2.0,0.084,0.0008105996084642976
upper_only,2.0,0.085,0.0008616830769191356
upper_only,2.0,0.086,0.0009151085876345379
upper_only,2.0,0.08700000000000001,0.000967478405932972
upper_only,2.0,0.08800000000000001,0.001016060279414279
upper_only,2.0,0.08900000000000001,0.0010589685241530924
upper_only,2.0,0.09,0.0010951943064244508
upper_only,2.0,0.091,0.0011245112907720424
upper_only,2.0,0.09200000000000001,0.001147300387719068
upper_only,2.0,0.09300000000000001,0.001164341658651121
upper_only,2.0,0.094,0.0011766146888080205
upper_only,2.0,0.095,0.0011851353918069206
upper_only,2.0,0.09600000000000002,0.0011908421805556328
upper_only,2.0,0.09700000000000002,0.0011945316162446974
upper_only,2.0,0.098,0.001196835142286257
upper_only,2.0,0.099,0.0011982246757213787
upper_only,2.0,0.1,0.0011990347729318861

examples/opposed_demo.csv

mode,time_s,x_m,gap_m
opposed,0.0,0.0,0.0011983591139842064
opposed,0.0,0.001,0.0011969819487263438
opposed,0.0,0.002,0.001194619741886637
opposed,0.0,0.0030000000000000005,0.0011907037476159856
opposed,0.0,0.004,0.001184431706944576
opposed,0.0,0.005,0.001174730166071765
opposed,0.0,0.006000000000000001,0.0011602449709736348
opposed,0.0,0.007000000000000001,0.0011393808197069056
opposed,0.0,0.008,0.0011104106591224151
opposed,0.0,0.009000000000000001,0.0010716691943124725
opposed,0.0,0.01,0.0010218303209215667
opposed,0.0,0.011000000000000001,0.0009602464910602568
opposed,0.0,0.012000000000000002,0.0008873024750042737
opposed,0.0,0.013000000000000001,0.0008047132900860518
opposed,0.0,0.014000000000000002,0.0007156845989787151
opposed,0.0,0.015,0.0006248612307625304
opposed,0.0,0.016,0.0005380193343893057
opposed,0.0,0.017,0.0004615072021765831
opposed,0.0,0.018000000000000002,0.00040149889660854535
opposed,0.0,0.019000000000000003,0.0003631780285454026
opposed,0.0,0.02,0.00034999999999999983
opposed,0.0,0.021,0.0003631780285454027
opposed,0.0,0.022000000000000002,0.00040149889660854546
opposed,0.0,0.023000000000000003,0.0004615072021765833
opposed,0.0,0.024000000000000004,0.000538019334389306
opposed,0.0,0.025,0.0006248612307625304
opposed,0.0,0.026000000000000002,0.0007156845989787154
opposed,0.0,0.027000000000000003,0.0008047132900860523
opposed,0.0,0.028000000000000004,0.0008873024750042742
opposed,0.0,0.029000000000000005,0.000960246491060257
opposed,0.0,0.03,0.0010218303209215667
opposed,0.0,0.031,0.0010716691943124725
opposed,0.0,0.032,0.0011104106591224151
opposed,0.0,0.033,0.0011393808197069056
opposed,0.0,0.034,0.0011602449709736348
opposed,0.0,0.035,0.001174730166071765
opposed,0.0,0.036000000000000004,0.001184431706944576
opposed,0.0,0.037000000000000005,0.0011907037476159856
opposed,0.0,0.038000000000000006,0.001194619741886637
opposed,0.0,0.03900000000000001,0.0011969819487263438
opposed,0.0,0.04,0.0011983591139842064
opposed,0.0,0.04100000000000001,0.0011991353138329274
opposed,0.0,0.042,0.0011995583615201683
opposed,0.0,0.043,0.0011997813730990194
opposed,0.0,0.044000000000000004,0.0011998951016665262
opposed,0.0,0.045,0.0011999512177445714
opposed,0.0,0.046000000000000006,0.0011999780121148107
opposed,0.0,0.047,0.0011999903942043631
opposed,0.0,0.04800000000000001,0.0011999959326502165
opposed,0.0,0.049,0.0011999983307630121
opposed,0.0,0.05,0.0011999993360234003
opposed,0.0,0.051000000000000004,0.001199999744014217
opposed,0.0,0.052000000000000005,0.0011999999043451014
opposed,0.0,0.053000000000000005,0.0011999999653560915
opposed,0.0,0.054000000000000006,0.0011999999878388443
opposed,0.0,0.055,0.001199999995862375
opposed,0.0,0.05600000000000001,0.001199999998635556
opposed,0.0,0.057,0.0011999999995638974
opposed,0.0,0.05800000000000001,0.0011999999998649016
opposed,0.0,0.059000000000000004,0.001199999999959436
opposed,0.0,0.06,0.001199999999988195
opposed,0.0,0.061000000000000006,0.0011999999999966703
opposed,0.0,0.062,0.0011999999999990896
opposed,0.0,0.063,0.0011999999999997588
opposed,0.0,0.064,0.0011999999999999379
opposed,0.0,0.065,0.0011999999999999845
opposed,0.0,0.066,0.0011999999999999962
opposed,0.0,0.067,0.001199999999999999
opposed,0.0,0.068,0.0011999999999999997
opposed,0.0,0.069,0.0012
opposed,0.0,0.07,0.0012
opposed,0.0,0.07100000000000001,0.0012
opposed,0.0,0.07200000000000001,0.0012
opposed,0.0,0.07300000000000001,0.0012
opposed,0.0,0.07400000000000001,0.0012
opposed,0.0,0.075,0.0012
opposed,0.0,0.07600000000000001,0.0012
opposed,0.0,0.077,0.0012
opposed,0.0,0.07800000000000001,0.0012
opposed,0.0,0.079,0.0012
opposed,0.0,0.08,0.0012
opposed,0.0,0.081,0.0012
opposed,0.0,0.08200000000000002,0.0012
opposed,0.0,0.083,0.0012
opposed,0.0,0.084,0.0012
opposed,0.0,0.085,0.0012
opposed,0.0,0.086,0.0012
opposed,0.0,0.08700000000000001,0.0012
opposed,0.0,0.08800000000000001,0.0012
opposed,0.0,0.08900000000000001,0.0012
opposed,0.0,0.09,0.0012
opposed,0.0,0.091,0.0012
opposed,0.0,0.09200000000000001,0.0012
opposed,0.0,0.09300000000000001,0.0012
opposed,0.0,0.094,0.0012
opposed,0.0,0.095,0.0012
opposed,0.0,0.09600000000000002,0.0012
opposed,0.0,0.09700000000000002,0.0012
opposed,0.0,0.098,0.0012
opposed,0.0,0.099,0.0012
opposed,0.0,0.1,0.0012
opposed,1.0,0.0,0.0012
opposed,1.0,0.001,0.0012
opposed,1.0,0.002,0.0011999999999999997
opposed,1.0,0.0030000000000000005,0.001199999999999999
opposed,1.0,0.004,0.0011999999999999962
opposed,1.0,0.005,0.0011999999999999845
opposed,1.0,0.006000000000000001,0.0011999999999999379
opposed,1.0,0.007000000000000001,0.0011999999999997588
opposed,1.0,0.008,0.0011999999999990896
opposed,1.0,0.009000000000000001,0.0011999999999966703
opposed,1.0,0.01,0.001199999999988195
opposed,1.0,0.011000000000000001,0.001199999999959436
opposed,1.0,0.012000000000000002,0.0011999999998649016
opposed,1.0,0.013000000000000001,0.0011999999995638974
opposed,1.0,0.014000000000000002,0.001199999998635556
opposed,1.0,0.015,0.001199999995862375
opposed,1.0,0.016,0.0011999999878388443
opposed,1.0,0.017,0.0011999999653560915
opposed,1.0,0.018000000000000002,0.0011999999043451014
opposed,1.0,0.019000000000000003,0.001199999744014217
opposed,1.0,0.02,0.0011999993360234003
opposed,1.0,0.021,0.0011999983307630121
opposed,1.0,0.022000000000000002,0.0011999959326502165
opposed,1.0,0.023000000000000003,0.0011999903942043631
opposed,1.0,0.024000000000000004,0.0011999780121148107
opposed,1.0,0.025,0.0011999512177445714
opposed,1.0,0.026000000000000002,0.0011998951016665262
opposed,1.0,0.027000000000000003,0.0011997813730990194
opposed,1.0,0.028000000000000004,0.0011995583615201683
opposed,1.0,0.029000000000000005,0.0011991353138329274
opposed,1.0,0.03,0.0011983591139842064
opposed,1.0,0.031,0.0011969819487263438
opposed,1.0,0.032,0.001194619741886637
opposed,1.0,0.033,0.0011907037476159856
opposed,1.0,0.034,0.001184431706944576
opposed,1.0,0.035,0.001174730166071765
opposed,1.0,0.036000000000000004,0.0011602449709736348
opposed,1.0,0.037000000000000005,0.0011393808197069056
opposed,1.0,0.038000000000000006,0.0011104106591224151
opposed,1.0,0.03900000000000001,0.0010716691943124723
opposed,1.0,0.04,0.001021830320921567
opposed,1.0,0.04100000000000001,0.0009602464910602564
opposed,1.0,0.042,0.000887302475004274
opposed,1.0,0.043,0.0008047132900860525
opposed,1.0,0.044000000000000004,0.0007156845989787151
opposed,1.0,0.045,0.0006248612307625306
opposed,1.0,0.046000000000000006,0.0005380193343893055
opposed,1.0,0.047,0.0004615072021765833
opposed,1.0,0.04800000000000001,0.00040149889660854514
opposed,1.0,0.049,0.0003631780285454027
opposed,1.0,0.05,0.00034999999999999983
opposed,1.0,0.051000000000000004,0.0003631780285454027
opposed,1.0,0.052000000000000005,0.00040149889660854546
opposed,1.0,0.053000000000000005,0.0004615072021765833
opposed,1.0,0.054000000000000006,0.000538019334389306
opposed,1.0,0.055,0.0006248612307625301
opposed,1.0,0.05600000000000001,0.0007156845989787158
opposed,1.0,0.057,0.0008047132900860518
opposed,1.0,0.05800000000000001,0.0008873024750042744
opposed,1.0,0.059000000000000004,0.0009602464910602568
opposed,1.0,0.06,0.0010218303209215665
opposed,1.0,0.061000000000000006,0.0010716691943124725
opposed,1.0,0.062,0.0011104106591224151
opposed,1.0,0.063,0.0011393808197069056
opposed,1.0,0.064,0.0011602449709736348
opposed,1.0,0.065,0.001174730166071765
opposed,1.0,0.066,0.001184431706944576
opposed,1.0,0.067,0.0011907037476159856
opposed,1.0,0.068,0.001194619741886637
opposed,1.0,0.069,0.0011969819487263438
opposed,1.0,0.07,0.0011983591139842064
opposed,1.0,0.07100000000000001,0.0011991353138329274
opposed,1.0,0.07200000000000001,0.0011995583615201683
opposed,1.0,0.07300000000000001,0.0011997813730990194
opposed,1.0,0.07400000000000001,0.0011998951016665262
opposed,1.0,0.075,0.0011999512177445714
opposed,1.0,0.07600000000000001,0.0011999780121148107
opposed,1.0,0.077,0.0011999903942043631
opposed,1.0,0.07800000000000001,0.0011999959326502165
opposed,1.0,0.079,0.0011999983307630121
opposed,1.0,0.08,0.0011999993360234003
opposed,1.0,0.081,0.001199999744014217
opposed,1.0,0.08200000000000002,0.0011999999043451014
opposed,1.0,0.083,0.0011999999653560915
opposed,1.0,0.084,0.0011999999878388443
opposed,1.0,0.085,0.001199999995862375
opposed,1.0,0.086,0.001199999998635556
opposed,1.0,0.08700000000000001,0.0011999999995638974
opposed,1.0,0.08800000000000001,0.0011999999998649016
opposed,1.0,0.08900000000000001,0.001199999999959436
opposed,1.0,0.09,0.001199999999988195
opposed,1.0,0.091,0.0011999999999966703
opposed,1.0,0.09200000000000001,0.0011999999999990896
opposed,1.0,0.09300000000000001,0.0011999999999997588
opposed,1.0,0.094,0.0011999999999999379
opposed,1.0,0.095,0.0011999999999999845
opposed,1.0,0.09600000000000002,0.0011999999999999962
opposed,1.0,0.09700000000000002,0.001199999999999999
opposed,1.0,0.098,0.0011999999999999997
opposed,1.0,0.099,0.0012
opposed,1.0,0.1,0.0012
opposed,2.0,0.0,0.0012
opposed,2.0,0.001,0.0012
opposed,2.0,0.002,0.0012
opposed,2.0,0.0030000000000000005,0.0012
opposed,2.0,0.004,0.0012
opposed,2.0,0.005,0.0012
opposed,2.0,0.006000000000000001,0.0012
opposed,2.0,0.007000000000000001,0.0012
opposed,2.0,0.008,0.0012
opposed,2.0,0.009000000000000001,0.0012
opposed,2.0,0.01,0.0012
opposed,2.0,0.011000000000000001,0.0012
opposed,2.0,0.012000000000000002,0.0012
opposed,2.0,0.013000000000000001,0.0012
opposed,2.0,0.014000000000000002,0.0012
opposed,2.0,0.015,0.0012
opposed,2.0,0.016,0.0012
opposed,2.0,0.017,0.0012
opposed,2.0,0.018000000000000002,0.0012
opposed,2.0,0.019000000000000003,0.0012
opposed,2.0,0.02,0.0012
opposed,2.0,0.021,0.0012
opposed,2.0,0.022000000000000002,0.0012
opposed,2.0,0.023000000000000003,0.0012
opposed,2.0,0.024000000000000004,0.0012
opposed,2.0,0.025,0.0012
opposed,2.0,0.026000000000000002,0.0012
opposed,2.0,0.027000000000000003,0.0012
opposed,2.0,0.028000000000000004,0.0012
opposed,2.0,0.029000000000000005,0.0012
opposed,2.0,0.03,0.0012
opposed,2.0,0.031,0.0012
opposed,2.0,0.032,0.0011999999999999997
opposed,2.0,0.033,0.001199999999999999
opposed,2.0,0.034,0.0011999999999999962
opposed,2.0,0.035,0.0011999999999999845
opposed,2.0,0.036000000000000004,0.0011999999999999379
opposed,2.0,0.037000000000000005,0.0011999999999997588
opposed,2.0,0.038000000000000006,0.0011999999999990896
opposed,2.0,0.03900000000000001,0.0011999999999966703
opposed,2.0,0.04,0.001199999999988195
opposed,2.0,0.04100000000000001,0.001199999999959436
opposed,2.0,0.042,0.0011999999998649016
opposed,2.0,0.043,0.0011999999995638974
opposed,2.0,0.044000000000000004,0.001199999998635556
opposed,2.0,0.045,0.001199999995862375
opposed,2.0,0.046000000000000006,0.0011999999878388443
opposed,2.0,0.047,0.0011999999653560915
opposed,2.0,0.04800000000000001,0.0011999999043451014
opposed,2.0,0.049,0.001199999744014217
opposed,2.0,0.05,0.0011999993360234003
opposed,2.0,0.051000000000000004,0.0011999983307630121
opposed,2.0,0.052000000000000005,0.0011999959326502165
opposed,2.0,0.053000000000000005,0.0011999903942043631
opposed,2.0,0.054000000000000006,0.0011999780121148107
opposed,2.0,0.055,0.0011999512177445714
opposed,2.0,0.05600000000000001,0.0011998951016665262
opposed,2.0,0.057,0.0011997813730990194
opposed,2.0,0.05800000000000001,0.0011995583615201683
opposed,2.0,0.059000000000000004,0.0011991353138329274
opposed,2.0,0.06,0.0011983591139842064
opposed,2.0,0.061000000000000006,0.0011969819487263438
opposed,2.0,0.062,0.001194619741886637
opposed,2.0,0.063,0.0011907037476159856
opposed,2.0,0.064,0.001184431706944576
opposed,2.0,0.065,0.001174730166071765
opposed,2.0,0.066,0.0011602449709736348
opposed,2.0,0.067,0.0011393808197069056
opposed,2.0,0.068,0.0011104106591224151
opposed,2.0,0.069,0.0010716691943124723
opposed,2.0,0.07,0.0010218303209215665
opposed,2.0,0.07100000000000001,0.0009602464910602564
opposed,2.0,0.07200000000000001,0.0008873024750042733
opposed,2.0,0.07300000000000001,0.0008047132900860513
opposed,2.0,0.07400000000000001,0.0007156845989787145
opposed,2.0,0.075,0.0006248612307625306
opposed,2.0,0.07600000000000001,0.0005380193343893048
opposed,2.0,0.077,0.0004615072021765833
opposed,2.0,0.07800000000000001,0.0004014988966085448
opposed,2.0,0.079,0.0003631780285454027
opposed,2.0,0.08,0.00034999999999999983
opposed,2.0,0.081,0.0003631780285454027
opposed,2.0,0.08200000000000002,0.0004014988966085461
opposed,2.0,0.083,0.0004615072021765833
opposed,2.0,0.084,0.000538019334389306
opposed,2.0,0.085,0.0006248612307625306
opposed,2.0,0.086,0.0007156845989787145
opposed,2.0,0.08700000000000001,0.0008047132900860525
opposed,2.0,0.08800000000000001,0.0008873024750042744
opposed,2.0,0.08900000000000001,0.0009602464910602572
opposed,2.0,0.09,0.0010218303209215665
opposed,2.0,0.091,0.0010716691943124723
opposed,2.0,0.09200000000000001,0.0011104106591224155
opposed,2.0,0.09300000000000001,0.0011393808197069058
opposed,2.0,0.094,0.0011602449709736348
opposed,2.0,0.095,0.001174730166071765
opposed,2.0,0.09600000000000002,0.001184431706944576
opposed,2.0,0.09700000000000002,0.0011907037476159858
opposed,2.0,0.098,0.001194619741886637
opposed,2.0,0.099,0.0011969819487263438
opposed,2.0,0.1,0.0011983591139842064

• Contact and closure. An occlusive pump closes the tube: with a deforming mesh the gap cannot go to zero. movingCone shows the way (remap onto meshes with fewer cells in the gap), but a real closure requires a conventional minimum gap, declared and subjected to a sensitivity test.
• Deformable tube. The tube shape depends on the rollers and the material. G. Formato and co-authors, “Fluid–structure interaction modeling applied to peristaltic pump flow simulations”, Machines 7(3), 50, 2019, doi:10.3390/machines7030050 (open access), couple an ALE fluid to a neo-Hookean hyperelastic tube squeezed by a roller and compare the outlet pressure with experiments (difference below 5 %); they also report a time step dropping to 10⁻⁶ s when both ends of the tube are compressed.
• Prescribed motion as a first step. Before FSI a model with prescribed tube shape is advisable (for instance from the guide’s h(x,t) formula or from a measurement), which in OpenFOAM 14 is built with a point displacement field and one of the section 1 motion solvers; the lubrication theory of the previous box provides the expected flow and back-pressure values to check it.
• Measurements to plan. Instantaneous and mean flow rate, reflux (negative flow during part of the cycle), upstream and downstream pressure, and the volume displaced per cycle to compare with the geometric volume between two rollers.

The direct GCL test (section 2), i.e. a fluid at rest on a moving mesh, does not apply to movingCone, where the cone motion is precisely what drives the flow. The most effective substitute test is the exact balance on left, which checks both the GCL and the flux correction after the remappings.

cd "$HOME/OpenFOAM"
foamVersion > guide16_cases/OPENFOAM_VERSION.txt 2>&1
tar czf guide16_cases_$(date +%Y%m%d).tar.gz guide16_cases
sha256sum guide16_cases_*.tar.gz > guide16_cases_SHA256.txt

The archive must also contain constant/meshes and system/meshes: without the remapping meshes a restart from an intermediate time is not possible. No file is removed.

• “ERROR: Existing logs/results” on a freshly prepared case: it is the run_case.sh defect described in section 9; use extra/run_case_fixed.sh.
• Error reading constant/meshes/0.0015/polyMesh at the first event: the blockMesh -mesh commands were not run (or were run in another folder).
• Negative cell volumes during motion: the gap closed more than expected (amplitude or frequency changed without redesigning the meshes, section 4) or the diffusivity does not protect the cells; check with movingcone_schedule.py and the checkMesh function.
• Fields with extreme values just after an event: look at the initial residual of pcorr and p in the first step after “Mapping to mesh time”; if it stays high, reduce the step around the events in a copy.
• No “Mesh Courant Number” line in the log: it is the default behaviour; enable checkMeshCourantNo yes (enable_monitors.py).
• In parallel: the remapping meshes must be decomposed one by one, as Allrun-parallel does; decomposing only the initial mesh is not enough.

• C. W. Hirt, A. A. Amsden, J. L. Cook, J. Comput. Phys. 14, 227–253, 1974, doi:10.1016/0021-9991(74)90051-5 — the ALE method;
• P. D. Thomas, C. K. Lombard, AIAA J. 17, 1030–1037, 1979, doi:10.2514/3.61273 — the geometric conservation law;
• I. Demirdžić, M. Perić, Int. J. Numer. Methods Fluids 8, 1037–1050, 1988, doi:10.1002/fld.1650080906 — the space conservation law in finite volumes;
• C. Farhat, P. Geuzaine, C. Grandmont, J. Comput. Phys. 174, 669–694, 2001, doi:10.1006/jcph.2001.6932 — discrete GCL and stability of ALE schemes;
• H. Jasak, Ž. Tuković, Transactions of FAMENA 30(2), 1–20, 2006, CROSBI — automatic mesh motion for unstructured finite volumes;
• H. Jasak, «Dynamic mesh handling in OpenFOAM», AIAA 2009-341, doi:10.2514/6.2009-341;
• P. E. Farrell, J. R. Maddison, Comput. Methods Appl. Mech. Eng. 200, 89–100, 2011, doi:10.1016/j.cma.2010.07.015 — conservative interpolation between meshes;
• A. H. Shapiro, M. Y. Jaffrin, S. L. Weinberg, J. Fluid Mech. 37, 799–825, 1969, doi:10.1017/S0022112069000899; M. Y. Jaffrin, A. H. Shapiro, Annu. Rev. Fluid Mech. 3, 13–37, 1971, doi:10.1146/annurev.fl.03.010171.000305 — peristaltic pumping;
• G. Formato et al., Machines 7(3), 50, 2019, doi:10.3390/machines7030050 — FSI of a peristaltic pump, open access;
• J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, 4th ed., Springer, 2020, doi:10.1007/978-3-319-99693-6 — chapter on moving grids.

• tutorials/incompressibleFluid/movingCone — the complete case
• meshToMesh_fvMeshTopoChanger.C — remapping, cycle, repeat, timeDelta
• velocityComponentLaplacian_fvMotionSolver.C — motion solver
• directionalDiffusivity.C — directional diffusivity
• SquareI.H — square function
• movingWallVelocityFvPatchVectorField.C — moving wall condition
• fvMesh.C — meshPhi from swept volumes
• moveMesh.C, basicFluidSolver.C — correctPhi, Courant and mesh Courant numbers
• foamRun.C — order of operations and of log lines
• checkMesh.H — mesh check during the run