migflow.time_integration

exception migflow.time_integration.CouplingUnstable

Bases: RuntimeError

A runtime guard aborted the run: continuing would only burn CPU.

exception migflow.time_integration.CouplingWarning

Bases: UserWarning

A runtime guard fired: the coupled step looks unstable or unphysical.

migflow.time_integration.iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-08, external_particles_forces=None, fixed_grains=False, after_sub_iter=None, max_nsub=None, check_residual_norm=-1, use_predictor_corrector=False, mu=None, rho_f=None, gravity=None)

Migflow solver: the solver type depends on the fluid and particles given as arguments.

Parameters:

  • fluid – fluid structure

  • particles – particles structure

  • dt (float) – time step

  • min_nsub (int) – minimal nsub for the particles iterations

  • contact_tol (float) – tolerance for the contact solver

  • external_particles_forces (ndarray) – vector of external forces applied on the particles

  • fixed_grains (bool) – boolean variable specifying if the grains are fixed in the fluid

  • after_sub_iter (callable) – callback to execute once a sub iteration has been made

  • max_nsub (int) – maximum number of times the time step can be further split if conergence is not reached

  • check_residual_norm (float) – check if the fluid solver has converged to the specified norm

  • use_predictor_corrector (bool) – boolean variable specifying if the predictor-corrector scheme is used

  • mu (float) – the coupling closure’s parameters; None reads them off the fluid

  • rho_f (float) – the coupling closure’s parameters; None reads them off the fluid

  • gravity (ndarray) – the coupling closure’s parameters; None reads them off the fluid

Raises:

  • ValueError – fluid and particles cannot be both None

  • ValueError – external_particles_forces must have shape (number of particles,dimension)

migflow.time_integration.iterate_iqn(fluid, particles, dt, min_nsub=1, contact_tol=1e-08, external_particles_forces=None, after_sub_iter=None, max_nsub=None, check_residual_norm=-1, max_iter=10, tol=0.001, omega=0.5, depth=12, reuse=4, fixed_grains=False, set_bodies=None, get_forces=None, coupler=None, set_geometry=None, set_closure=None, geometry=None, mu=None, rho_f=None)

IQN-ILS accelerated fluid-grain step: quasi-Newton on the grain velocities.

The contact-dominated counterpart to iterate_patankar. It drives the fixed point v = G(v) where G(v) is “the velocity the DEM returns when the fluid force is computed with the grains moving at v” – contacts and all – and accelerates it with the IQN-ILS secant method in body_coupling.IQN. Because G contains the contact solve, an impulsive impact is reconciled with the fluid over the outer iterations rather than lagged; this is the scheme to reach for on deposits and impacts, where the explicit-contact Patankar scheme is fragile.

RETURNS THE COUPLER for reuse across steps. IQN-ILS reuses secant columns from the last reuse steps; clearing that history every step degrades it to barely-accelerated Picard (measured). Pass the returned coupler back in on the next step:

cpl = None for step in …:

cpl = iterate_iqn(fluid, particles, dt, …, coupler=cpl)

On the first call coupler is None and one is built here from max_iter/tol/omega/depth/reuse; thereafter the supplied coupler’s own settings govern and those keyword arguments are ignored.

BODY MASK. The coupler operates on ALL particles of nonzero radius: mask = particles.r().ravel() > 0. Walls have zero radius (and zero mass, so a mass division would blow up) and are excluded; the velocity guess is written back only into the masked slots. This keeps the indexing generic – no testcase-specific wall/boundary count leaks into the library.

Hooks for the unfitted path – EITHER set_bodies (one hook per iteration) OR the split pair set_geometry/set_closure (tier A once per step, tier B per iteration); get_forces is required with both:

  • set_bodies(fluid, particles, v_full=None) – hand the solver the body CSR for the grains moving at v_full (the current IQN velocity guess, a full p-indexed array). Geometry stays caller-side (architectural rule).

  • get_forces(fluid, particles) – the TOTAL non-contact force for the DEM, p-indexed, walls zero (fluid + gravity + external), as in iterate_patankar.

THE SPLIT HOOKS (optional, opt-in): set_geometry / set_closure

set_bodies is called once per OUTER ITERATION, so the whole body CSR – including the overlap build, which is the expensive part – is rebuilt every iteration. It does not have to be: restore_state() at the top of the loop puts the bodies back where they started, so WITHIN A STEP the positions and radii are identical across iterations and only the velocity guess moves. Only the closure (which reads the guess, hence the slip) is genuinely per-iteration.

A caller who can split its geometry from its closure says so by passing BOTH

  • set_geometry(fluid, particles) -> geo – built ONCE PER STEP, before the loop, at the start-of-step positions. geo is opaque to the library: whatever the caller’s closure hook needs (volume_coupling returns its tier-A dict, and has already handed it to the kernel).

  • set_closure(fluid, particles, geo, v_full) -> anything – called per iteration in set_bodies’ place, with the current velocity guess.

particles is passed to both because the closure legitimately reads particle state the geometry does not carry (contact forces); geo is passed through untouched.

geometry= hands the library a geo the caller built ITSELF, outside the time loop, in place of set_geometry: that is the held-bed case (darcy/2d_discs), where the bodies never move and tier A is invariant for the WHOLE RUN, not merely within a step. It is accepted only with fixed_grains=True – the iterating branch advances the grains, so a prebuilt geometry would be stale from the second step on, silently.

These are opt-in and mutually exclusive with set_bodies: given neither, the set_bodies path runs exactly as before.

Machinery per step: save_state; capture the fluid solution; start_step; then until done(): restore_state (the fluid must see the START-of-step geometry, only the velocity guess varying), set the guess into the masked slots, set_bodies, reset the fluid solution, implicit_euler, form the force, _advance_particles, take_target the resulting grain velocities.

migflow.time_integration.iterate_patankar(fluid, particles, dt, min_nsub=1, contact_tol=1e-08, external_particles_forces=None, fixed_grains=False, after_sub_iter=None, max_nsub=None, check_residual_norm=-1, set_bodies=None, get_forces=None, timers=None, prediction_contact_forces=True, patankar_datum=False, gravity=None, mu=None, rho_f=None)

Patankar (mixture-model) fluid-grain step: one fluid solve, one DEM advance.

This is the first-class library form of the Patankar scheme historically run through iterate(). iterate() itself is kept untouched for backward compatibility; this function duplicates its coupled statements verbatim (the set_bodies is None branch below) rather than delegating, and adds the unfitted (caller-supplied body CSR) route.

WHEN TO USE. patankar is the cheap coupling: a single fluid solve per step, no outer iteration. It is FRAGILE AT IMPULSIVE IMPACT – the contact force is carried explicitly into the DEM prediction and never reaches the fluid Jacobian (f00/f01/…), so contact and fluid cannot be reconciled inside the solve. That is inherent to the scheme, not a defect, and is exactly why iterate_iqn exists for contact-dominated problems (deposits, impacts).

THE CONTACT FORCE IS THE DEM SOLVER’S, BY DEFAULT (owner ruling 2026-08-16): the prediction carries the real contact forces, which is the physical coupling and converges to the contact-consistent hydrostatic state at temporal convergence. STABILITY IS THE TIME STEP’S JOB: the lagged force injects O(dt*Fc/m) of velocity per step, so dt must keep that below the flow scale – scale dt with grain density (the steel depot needs dt = 1e-3 where its rho_p = 1500 sibling takes 5e-3; fine-grain beds like vortices r = 2e-5 need dt sized accordingly or the IQN scheme). prediction_contact_forces=False is the measured-legacy parity setting: migflow-main’s contact_forces getter always returned zeros (per-contact native vector + zero-filling shape guard), so the old model never saw a contact force – that is what the legacy-retrieval configuration reproduces, not the physical default.

Two paths, selected by set_bodies:

  • set_bodies is None – LIBRARY PATH: tier A geometry, tier B closure with model="unfitted" (m = (x - c)/r, div_m = dim/r, the Babuska penalty blended on resolution into the Dallavalle law, trace band and both extra diffusivities live) and the tier C Patankar datum, which condenses the grain momentum into the single solve exactly as the old in-kernel prediction did. mu/rho_f/gravity are the closure’s parameters and are read off the fluid when not given.

  • set_bodies given – UNFITTED PATH. Geometry stays caller-side (the core never builds overlap geometry – an architectural rule): set_bodies(fluid, particles, v_full=None) hands the solver its body CSR, with v_full=None meaning “use particles’ own velocity”. get_forces(fluid, particles) must return the TOTAL non-contact force to hand the DEM, p-indexed (walls included, zeros there): fluid + gravity + any external. This mirrors the compute_node_force convention where the caller adds g*m outside – here the whole non-contact force composition is the caller’s responsibility, so the library never assumes where gravity or buoyancy enter. The flow is set_bodies -> implicit_euler -> F = get_forces + external -> advance.

    patankar_datum=True makes THE SCHEME apply its own semi-implicit treatment inside that single solve: set_bodies must return the tier dicts, and the library builds the affine datum from them (needs gravity). get_forces is then just the fluid force plus gravity – no caller-side 1/(1 + gamma dt/m). Without it the branch is explicit and the case must do that correction itself, which is what made “patankar” name two different schemes.

Args mirror iterate(); set_bodies/get_forces are the unfitted hooks.

migflow.time_integration.predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-08, external_particles_forces=None, alpha=0.5, after_sub_iter=None, max_nsub=None, check_residual_norm=-1, mu=None, rho_f=None, gravity=None)

Predictor-corrector scheme to solve fluid and grains.

Parameters:

  • fluid – fluid structure

  • particles – particles structure

  • dt (float) – time step

  • min_nsub (int) – minimal nsub for the particles iterations

  • contact_tol (float) – tolerance for the contact solver

  • external_particles_forces (ndarray) – external forces applied on the particles

  • alpha (float) – parametre of the predictor-corrector scheme [alpha*f(n)+(1-alpha)*f(n+1)]

  • after_sub_iter (callable) – callback to execute once a sub iteration has been made

  • max_split – maximum number of times the time step can be further split if convergence is not reached

  • check_residual_norm (float) – check if the fluid solver has converged to the specified norm

  • mu (float) – the coupling closure’s parameters; None reads them off the fluid

  • rho_f (float) – the coupling closure’s parameters; None reads them off the fluid

  • gravity (ndarray) – the coupling closure’s parameters; None reads them off the fluid

migflow.time_integration.reset_guards()

Reset the guard counters (a new run in the same process).