PyIMR

Fast, validated solvers for inertial microcavitation rheometry (IMR) with a typed, dimensional material API.

The package is designed for inference campaigns that require many forward evaluations. It retains closed-form hot paths for common constitutive laws and also supports composable hyperelastic and generalized-Newtonian materials, plus distributed Giesekus and linear PTT memory.

python -m pip install -e ".[test]"
python -m pytest                  # everything, including numerical validation
python -m pytest -m "not slow"    # skip the high-resolution convergence studies

pytest prints a table of the measured deviations after the run, not only pass/fail — a check that still passes but has moved an order of magnitude is worth seeing.

Contents

Solver scope

pyimr.simulate supports:

Unsupported option values and inconsistent thermal/mass-transfer combinations fail during configuration. Integration failures and material-domain violations raise SimulationError.

Public API

All inputs are dimensional. The material is explicit and typed; there is no integer constitutive selector or shared bag of material parameters.

import numpy as np
from pyimr import (
    NeoHookeanKelvinVoigt,
    SimulationConfig,
    prepare,
    simulate,
)

t = np.linspace(0.0, 120e-6, 300)
config = SimulationConfig(
    R0=225e-6,
    Req=37.5e-6,
    material=NeoHookeanKelvinVoigt(
        shear_modulus_pa=2500.0,
        viscosity_pa_s=0.1,
    ),
)
result = simulate(t, config)

result.time_s
result.radius_ratio
result.radius_m
result.wall_velocity_m_s
result.internal_pressure_pa
result.stress_integral_pa
result.stats

The returned arrays are read-only. Thermal configurations also return gas and liquid temperature fields and, when enabled, vapor mass fraction. Materials with memory expose their internal stress state. Inactive fields are None.

For repeated solves with one configuration, preparation hoists constant work such as state layout, grids, finite-difference operators, constitutive quadrature, and Jacobian sparsity:

problem = prepare(config)
first = problem.solve(t)
second = problem.solve(t)  # immutable setup is reused; solve state is fresh

The same prepared problem can integrate any set of continuous parameter directions with the state:

sensitivity = problem.solve_with_sensitivities(
    t,
    (
        "R0",
        "material.shear_modulus_pa",
        "material.viscosity_pa_s",
        "physics.polytropic_exponent",
    ),
)

sensitivity.simulation
sensitivity.radius_m       # shape: (time, parameter)
sensitivity.state          # complete internal-state derivatives

Parameter paths follow the frozen configuration objects. Returned derivatives are with respect to dimensional parameter values. The sensitivity result also contains wall-velocity, pressure, stress, gas-temperature, medium-temperature, and vapor-fraction derivatives when those fields are active.

Composable instantaneous materials

An InstantaneousMaterial can contain an elastic law, a viscous law, or both:

from pyimr import CarreauYasuda, Gent, InstantaneousMaterial

material = InstantaneousMaterial(
    elastic=Gent(
        shear_modulus_pa=2500.0,
        extensibility=250.0,
    ),
    viscous=CarreauYasuda(
        zero_shear_viscosity_pa_s=0.5,
        infinite_shear_viscosity_pa_s=0.02,
        time_constant_s=20e-6,
        transition_exponent=2.0,
        power_index=0.45,
    ),
)
config = SimulationConfig(R0=225e-6, Req=37.5e-6, material=material)

Elastic laws:

• NeoHookean
• MooneyRivlin
• Yeoh
• Fung
• Gent
• ArrudaBoyce
• Ogden

Generalized-Newtonian laws:

• Newtonian
• PowerLaw
• CarreauYasuda
• Cross
• PowellEyring
• ModifiedPowellEyring
• HerschelBulkley
• Bingham

Carreau is CarreauYasuda(transition_exponent=2); simplified Cross is Cross(transition_exponent=1).

Ogden takes matched tuples and is the only elastic law here that depends on the principal stretches rather than on I1 alone, so exponents may be negative or fractional:

from pyimr import Ogden

material = Ogden(shear_moduli_pa=(1800.0, 600.0, -300.0), exponents=(1.3, 4.0, -2.0))

A single term with exponents=(2.0,) is neo-Hookean, and reduces to it exactly rather than asymptotically. The small-strain shear modulus is sum(shear_moduli_pa * exponents) / 2, which must be positive; individual moduli may be negative, as in the example above.

BlatzKo is deliberately absent. The standard Blatz-Ko strain energy is distinguished by its dependence on I3, and this solver assumes an incompressible spherical deformation with stretches (l^-2, l, l), so I3 = 1 identically. In that limit Blatz-Ko is MooneyRivlin(c10=0.0, c01=mu/2) -- verified equal to machine precision -- so a separate class would be an alias, not a new capability. A genuinely compressible Blatz-Ko needs the incompressibility assumption relaxed throughout the radial dynamics.

The Powell-Eyring pair uses the standard laws, eta_inf + (eta_0 - eta_inf) * asinh(x)/x and its log1p(x)/x variant with x = lambda*|gdot|. Both reduce exactly to Newtonian(eta_0) as lambda -> 0. IMRv2's f_viscosity.m instead uses sinh(x)/x^nc, which is shear-thickening and diverges exponentially, and a log(1+x)/x^nc variant with no finite zero-shear limit unless nc == 1; neither was copied.

Prepared Gauss-Legendre rules evaluate the finite-interval stress integrals. The solver evaluates the stress-rate terms and acceleration coefficient analytically, including the viscosity tangent. No finite-difference derivative is used inside the radial dynamics. The specialized NeoHookeanKelvinVoigt path and the equivalent composable material agree to solver tolerance.

PowerLaw, HerschelBulkley, and Bingham require a positive regularization_rate_per_s. The latter two use a smooth yield-stress regularization so the implicit radial equations retain a finite tangent at zero strain rate.

Gent lock-up is a material-domain error. If a trajectory reaches $I_1-3\ge J_m$, the solve stops and raises SimulationError instead of continuing with nonphysical stress.

Closed-form memory models

The finite-dimensional hot paths are:

• Zener
• QuadraticZener
• OldroydB
• LinearMaxwell

For example:

from pyimr import Zener

material = Zener(
    shear_modulus_pa=2500.0,
    viscosity_pa_s=0.1,
    relaxation_time_s=40e-6,
    retardation_time_s=8e-6,
)

Distributed nonlinear memory

Giesekus and LinearPTT evolve radial and hoop stress on a prepared, wall-clustered Lagrangian grid:

from pyimr import Giesekus

material = Giesekus(
    viscosity_pa_s=0.1,
    relaxation_time_s=40e-6,
    retardation_time_s=8e-6,
    mobility=0.2,
)
config = SimulationConfig(R0=225e-6, Req=37.5e-6, material=material)

result.stress_state contains radial stress followed by hoop stress on result.stress_reference_radius_ratio. Prepared coupled heat/mass-transfer problems use sparse BDF, while non-stiff configurations retain LSODA.

The constitutive equations at each material point are ordinary differential equations -- there are no spatial derivatives -- so the only spatial approximation is the quadrature for the stress integral. quadrature="gauss" (the default, 240 points) places the material points at Gauss-Legendre nodes and converges spectrally; it is about five orders of magnitude more accurate than the trapezoid grid it replaced, and cheaper; the table is in docs/discretization.md. See docs/discretization.md.

At zero mobility or zero extensibility, the distributed models converge to the analytic Oldroyd-B solution. Use OldroydB directly when that closure applies: it is substantially smaller and faster.

Physical parameters, initial state, and forcing

from pyimr import InitialState, PhysicalParameters

config = SimulationConfig(
    R0=225e-6,
    Req=37.5e-6,
    material=NeoHookeanKelvinVoigt(2500.0, 0.1),
    physics=PhysicalParameters(polytropic_exponent=1.47),
    initial=InitialState(
        wall_velocity_m_s=2.0,
        internal_pressure_pa=1.5e5,
    ),
)

Sampled forcing values are pressure perturbations relative to the far-field baseline. A shape-preserving cubic interpolant is used between samples and the perturbation is zero outside their time span.

from pyimr import SampledForcing

config = SimulationConfig(
    R0=225e-6,
    Req=37.5e-6,
    material=NeoHookeanKelvinVoigt(2500.0, 0.1),
    sampled_forcing=SampledForcing(
        time_s=tuple(measured_time_s),
        pressure_pa=tuple(measured_pressure_perturbation_pa),
    ),
)

Collapse-state shooting

Memory materials can initialize from a resolved equilibrium-to-maximum-radius precursor rather than an assumed unstressed state:

from pyimr import CollapseInitialization

config = SimulationConfig(
    R0=225e-6,
    Req=37.5e-6,
    material=Zener(2500.0, 0.1, 40e-6, 8e-6),
    collapse=CollapseInitialization(),
)
problem = prepare(config)
problem.collapse_stats

Preparation brackets the precursor velocity, shoots to R/R0 == 1, and retains the complete memory state at the maximum. CollapseStats records the root, achieved maximum, integration work, and immutable stress state. Sensitivities differentiate the event time and shooting root implicitly. Oldroyd-B and distributed Giesekus/PTT states use the same mechanism; the Zener precursor retains the upstream IMRv2 formulation.

Sensitivities and inference

The tangent-linear solver differentiates the production RHS rather than a reduced surrogate. It covers radial models 1--5, every typed material, thermal and mass-transfer states, distributed nonlinear memory, forcing, geometry, initial conditions, and continuous physical parameters.

The mechanical path, including distributed memory, uses a cached compiled directional kernel. After its one-time compilation, six simultaneous NHKV gradients take about 1.9 times one prepared forward solve on the development machine. Thermal and sampled-forcing branches use the same forward-mode reference implementation and augmented sparse BDF structure where appropriate.

Gradient accuracy

Tangent accuracy is not uniform across configurations. The mechanical path is limited by the finite-difference check it is measured against; the coupled thermal paths are limited by how accurately the augmented state/tangent system is integrated, which is a real bound on anything built on those gradients rather than a defect in the tangent equations. Measured error by configuration, and what tightening tolerance buys, are in docs/accuracy.md.

Prepared inference uses normalized bounded coordinates, dimensional Gaussian radius likelihoods, analytic sensitivity Jacobians, deterministic Latin-hypercube starts, and optional process-parallel batch evaluation:

from pyimr.inference import (
    InferenceParameter,
    RadiusObservation,
    prepare_inference,
)

inference = prepare_inference(
    config,
    RadiusObservation(
        measured_time_s,
        measured_radius_m,
        standard_deviation_m=2e-6,
    ),
    (
        InferenceParameter(
            "material.shear_modulus_pa", 500.0, 5000.0, "log"
        ),
        InferenceParameter(
            "material.viscosity_pa_s", 0.01, 1.0, "log"
        ),
    ),
)

batch = inference.evaluate_batch(unit_parameter_matrix, workers=4)
fit = inference.fit_multistart(64, seed=7, workers=4)
fit.endpoints  # every successful and unsuccessful endpoint is retained
fit.best

Unit parameter vectors always lie in [0, 1]; the configured linear or logarithmic transform maps them to physical bounds. Multistart results never discard alternative basins.

Thermal discretization

thermal="spectral" -- the default -- puts both the gas and liquid grids on Chebyshev collocation. thermal="fd" is the cheaper second-order finite difference, and is what every pinned IMRv2 trajectory was generated against, so the pinned suite sets it explicitly regardless of the default.

On the fully coupled model, spectral at Nt = 25 matches finite difference at Nt = 200 for a ninth of the cost. Neither is converged there, and choosing the scheme is a separate decision from choosing the resolution -- the default makes only the first. Both cost/accuracy tables are in docs/discretization.md.

Choosing solver tolerances

Observables do not converge at the same rate, so the right tolerance depends on which one you fit -- internal pressure is roughly two orders behind radius at the same setting. rtol=1e-6, atol=1e-8 is ample for likelihood evaluation against experimental radius data; keep 1e-10, 1e-12 for sensitivities and 1e-9 or tighter for validation. The per-observable table is in docs/accuracy.md.

Tolerance does not bound how long one solve can take. max_steps (default 1_000_000) turns a trajectory that will not finish into a SimulationError at a point of your choosing, which is what makes a grid sweep affordable.

Choosing a resolution

Nt and tolerance requirements depend on record length, material stiffness and which observable is fitted, so a setting that is adequate for one collapse can be badly wrong over five. pyimr.resolution measures on your own problem rather than guessing.

from pyimr.resolution import choose_resolution

setting = choose_resolution(config, times, target=1e-3, field="radius_ratio")
# Resolution(thermal='fd', Nt=5, rtol=1e-06, atol=1e-08,
#            achieved=3.9e-05, seconds=0.0037)

config = setting.apply(config)

It builds a reference and checks it is converged, searches both spectral and fd for the cheapest grid meeting target, then loosens tolerance as far as that grid allows. Roughly 18 solves, which is worth paying before a sampling or sensitivity campaign and not worth paying for a single run.

target is relative to each field's own peak magnitude, and field accepts several names. That matters because observables do not converge together: at identical settings relative error was 3.4e-07 for radius and 2.8e-05 for internal pressure.

It raises rather than guessing if the reference is not converged or the target is out of reach, since a number built on either is indistinguishable from a real answer.

Model selection

Constitutive models for soft matter nest -- NHKV is qKV at zero strain stiffening and SLS at zero relaxation time -- so comparing best fits always favours the flexible ones. pyimr.selection scores them by evidence instead.

from pyimr.selection import STANDARD_MODELS, compare, log_evidence, redundancy_over_grid, solve_grid

evidences = {}
for candidate in STANDARD_MODELS.values():
    points, normalized, radii, stresses = solve_grid(candidate, solve, count=12)
    redundancies = redundancy_over_grid(candidate, STANDARD_MODELS, points, stresses, solve)
    evidences[candidate.name], _ = log_evidence(
        radii, normalized, redundancies, observed, deviations, dimension=candidate.dimension
    )
posterior = compare(evidences)

pyimr.noise supplies the strain-rate weighting and the marginalized noise scale; pyimr.prior the redundancy and Occam penalties. Use ONE grid count for every model compared -- mixed resolutions let grid luck decide which lands nearest the truth.

Always report the best chi-squared per sample alongside the posterior. Model selection only means something where some candidate actually fits; otherwise the winner is the least-bad member of an inadequate set, and the posteriors look just as confident.

Worked studies are in examples/.

Trace estimators

pyimr.data covers the step before inference: getting from a measured R(t) history to the quantities a fit needs.

from pyimr import data

Req = data.equilibrium_radius(R0_m, initial_gas_pressure_pa)
omega_n, beta = data.natural_frequency(R0_m, Req, 2500.0, 0.1)
collapse_times_s, peak_radii_m, peak_times_s = data.collapse_features(
    measured_time_s, measured_radius_m
)
data.resolution_convergence(config, times_s, [10, 20, 40])

equilibrium_radius inverts the solver's own pressure/radius relation exactly. natural_frequency linearises Rayleigh-Plesset about Req in a Kelvin-Voigt medium; it reproduces Minnaert exactly in the gas-only limit and matches the simulated rebound frequency closely. collapse_features locates interior extrema with sub-sample parabolic refinement, replacing the manual index windows of IMR-vanilla calc_3tmins_3Rmaxs. resolution_convergence reports a table for a ladder you supply, where pyimr.resolution searches for a setting; both scale deviations by the field's own peak and both move Mt with Nt only when the medium is actually solved. Pass (Nt, Mt) pairs to set both.

IMR-vanilla's calc_omega_N is deliberately not ported: it is a scratch script whose formula treats the gas pressure at Rmax as the equilibrium value, inflating the stiffness by alpha**(-3*kappa); see docs/upstream.md. Video processing (calcRofT/) is also out of scope -- that is image analysis, and scikit-image covers it.

Validation

The suite pins IMRv2 trajectories across radial equations, forcing, vapor, heat transfer, mass transfer and the specialized constitutive models, and separately checks closed forms, reduction limits, and every analytic tangent against independent centered differences.

Two statistics are reported per pinned case, because they measure different things (#23): the pointwise maximum sits at a collapse in every case and is dominated by sub-nanosecond integrator phase, while the median carries no such sensitivity and is what the suite bounds tightly. The per-case deviation tables and the argument behind that split are in docs/validation.md.

Boundaries

• PhysicalParameters defaults reproduce the pinned reference trajectories. The default polytropic exponent is 1.4; IMRv2 itself ships with 1.47.
• SimulationResult.stress_state contains nondimensional internal variables. Public dimensional outputs carry units in their names or documentation.
• radial = 6 (Gilmore/Mie-Gruneisen) is supported here, and is the one configuration IMRv2 cannot run at all -- upstream returns complex radii without raising. The cause is a wrong root of the Mie-Gruneisen density quadratic; see docs/upstream.md.
• radial = 5 and radial = 6 deliberately diverge from IMRv2. Upstream's Mie-Gruneisen branch is physically wrong; the corrections are validated against the independent Tait branches and the weakly-compressible limit rather than against upstream. tests/ref_radial5.csv is retained as a record of upstream behaviour, not as a target.
• Collapse shooting requires a material with memory and cannot be combined with an explicit initial stress state or nonzero observed wall velocity. IMRv2 does permit collapse for memoryless materials; see PLAN.md W8.

Reference implementation

PyIMR reproduces IMRv2 except where upstream is wrong. One divergence is numerical rather than a defect fix: the Zener acceleration coefficient, which makes three Zener reference trajectories regression pins rather than cross-checks (#174).

PyIMR diverges from IMRv2 in several places, always deliberately. Eight defects were found at dea31cd, each reproduced with MATLAB R2025a via tools/gen_imrv2_cases.m, and each correction validated against something other than upstream -- a closed form, an independent equation of state, or a reduction limit. The wrong Mie-Gruneisen root, the non-functional non-Newtonian viscosity suite, the stubbed collapse initialization and the rest are in docs/upstream.md.

Those defects are why several PyIMR models are validated by reduction limit rather than against a pinned upstream trajectory: for those models, no working upstream implementation exists to pin against.

Tangent equations

Forward sensitivities integrate:

$$ \frac{ds_k}{dt}=J_y s_k+\frac{\partial f}{\partial c_k}. $$

All requested parameter directions share one augmented integration. Prepared parameter scaling keeps error control dimensionless; public derivatives are converted back to dimensional parameter units.

Citation

• Estrada, Barajas, Henann, Johnsen & Franck, High strain-rate soft material characterization via inertial cavitation, JMPS (2018). https://doi.org/10.1016/j.jmps.2017.12.006
• Warnez & Johnsen, Numerical modeling of bubble dynamics in viscoelastic media with relaxation, Physics of Fluids 27, 063103 (2015). https://doi.org/10.1063/1.4922598
• IMRv2: https://github.com/InertialMicrocavitationRheometry/IMRv2

License

MIT — see LICENSE.