migflow.fluid

C++ Python bindings for the MigFlow fluid solver.

class migflow.fluid.FluidProblem2(g=None, mu=0.0, rho=0.0, coeff_stab=0.01, volume_drag=0.0, quadratic_drag=0.0, drag_in_stab=False, drag_coefficient_factor=1.0, temporal=True, advection=True, flag_div_us=True, p2p1=False, full_implicit=False, model_b=False, density_element='', viscosity_element='', solver=None, solver_options='')

Bases: object

MigFlow incompressible fluid solver (2-D or 3-D). Exposes a high-level C++/Python API. Call set boundary conditions, then iterate with implicit_euler().

Construct a FluidProblem. Mirrors Python FluidProblem.__init__().

Parameters

gUnion[None, NDArray[np.float64], List[float]]

Gravity vector, length dim. Default: zero vector.

mufloat

Dynamic viscosity, the initial uniform value of the viscosity field

rhofloat

Density, the initial uniform value of the density field

coeff_stabfloat

Stabilisation coefficient

volume_dragfloat

Volume drag coefficient

quadratic_dragfloat

Quadratic drag coefficient

drag_in_stabbool

Include the coupling terms in the SUPG/PSPG residual. In the mixture model this is the drag term. In the UNFITTED model it is the Babuska penalty (trick 2) and the reconstructed-trace band (trick 3): the stabilisation block is applied after them, so the flag selects between the Galerkin-only snapshot and the full residual. True reproduces what the legacy mixture kernel did, and is required under model=’legacy’, where beta = gamma_D/(mu V_b) is a physical drag coefficient – excluding it leaves the pressure row blind to the drag and a fluidised bed collapses (2d_bubbling_fluidized_bed, drag/weight saturating at 0.61). False restores the earlier unfitted-only rule, under which betap*pspg = gamma_B/4 is mesh-independent so a stabilised Babuska penalty never refines away – keep it False when measuring convergence orders of the pure unfitted model. DEFAULT IS FALSE, and the two regimes genuinely disagree, so this is a modelling choice the case must make rather than inherit. Measured 2026-08-21 on one build: flag 2d_stokes_noslip_mms force error 2d_bubbling_fluidized_bed 0 244% -> 63% -> 16% (converging) bed COLLAPSES 4.67 -> 2.24 cm 1 1023% -> 1368% -> 1643% (diverging) bed FLUIDISES 4.67 -> 13.2 cm A sub-cell grain carries a physical drag that belongs in the pressure row; a resolved grain carries a mesh penalty that must stay out of it. The blend beta = s*gamma_B + (1-s)*gamma_D already interpolates between those two regimes per entry, so the principled successor to this global flag is a per-entry weight of (1-s) – not yet implemented.

drag_coefficient_factorfloat

Factor multiplying the drag coefficient

temporalbool

Enable temporal (d/dt) term

advectionbool

Enable advective terms (Navier-Stokes vs Stokes)

flag_div_usbool

Enable div(u_solid) term

p2p1bool

Use P2P1 (Taylor-Hood) elements instead of stabilised P1P1

full_implicitbool

Use fully implicit nonlinear scheme

model_bbool

Enable model B

density_elementstr

Density discretisation element (empty = default)

viscosity_elementstr

Viscosity discretisation element (empty = default)

solverUnion[None, str, ILinearSystem]

Optional plug-in linear solver object or solver name: petsc, pardiso, cudss, amgx, or default.

solver_optionsstr

Options forwarded to the named solver: a PETSc options string for petsc/default, an AmgX JSON configuration or config file path for amgx.

Pointer

alias of LP__Structure

adapt_mesh(nodes, elements, element_tags, boundaries, boundary_tags, boundary_names, periodic=None)

Adapt the mesh, projecting the current solution onto the new mesh.

Parameters

nodesNDArray[np.float64]

Node coordinates, shape (n_nodes, 3)

elementsNDArray[np.int32]

Element connectivity, shape (n_elements, dim+1)

element_tagsNDArray[np.int32]

Element tags, shape (n_elements,)

boundariesNDArray[np.int32]

Boundary edge nodes, shape (n_bnd, dim)

boundary_tagsNDArray[np.int32]

Boundary tags, shape (n_bnd,)

boundary_namesList[str]

Physical group names

periodicNDArray[np.int32]

Parent DOF index for each node (n_nodes,), or None

add_constraint(dofs, weights, rhs)

Add the linear constraint sum_i weights[i] * solution[dofs[i]] = rhs.

Parameters

dofsNDArray[np.int32]

Solution-vector indices

weightsNDArray[np.float64]

One weight per index

rhsfloat

Right-hand side

advance_concentration(dt)

Advance the concentration field by one time step (two-fluid problems).

Parameters

dtfloat

Time step

assemble_local_system(dt, localv, localm)

Assemble the element-local residual and Jacobian at the CURRENT solution, without solving. DIAGNOSTIC ENTRY POINT. It exists so a finite-difference check can ask whether the assembled Jacobian really is the derivative of the assembled residual: perturb one dof, reassemble, and compare the difference quotient against the matrix column. Without it the two can only be compared through a solve, which conflates a wrong derivative with an ill-conditioned one. Call it AFTER at least one implicit_euler, which is what sets the mesh metrics, the boundary conditions and the old-solution snapshot; this method deliberately refreshes none of those, so sold stays fixed while the current solution is perturbed.

Parameters

dtfloat

Time step

localvRealDeviceArray

Output local right-hand side, (n_elements, local_size)

localmRealDeviceArray

Output local matrix, (n_elements, local_size, local_size)

boundary_force_by_contributions(tag)

Return the total forces at a named boundary split into pressure and viscous. Shape: (2*dim,) = [p_x,..,p_dim, v_x,..,v_dim].

Return type:

ndarray[Any, dtype[float64]]

Parameters

tagstr

Boundary name

Returns

forcesNDArray[np.float64]

Force array, size 2*dim

boundary_forces()

Return the accumulated weak-boundary force array. Shape: (n_boundary_edges, dim). Returns ——-

forcesRealDeviceArray

Per-boundary-edge force array

Return type:

RealDeviceArray

bulk_force()

Return the bulk force field on the P1 nodes, shape (ndof_p1, dim). A body force per unit volume, added to the momentum equation. Note that compute_capillary_force() OVERWRITES it whenever sigma is nonzero, so a caller setting it directly must either leave surface tension off or set it again after each call. Returns ——-

bulk_forceRealDeviceArray

Bulk force field

Return type:

RealDeviceArray

clear_constraints()

Drop every constraint added by add_constraint(). A caller that rebuilds its constraints at each time step must call this first, otherwise they accumulate.

compute_node_force(dt=-1)

Compatibility alias for get_forces_on_bodies().

Return type:

ndarray[Any, dtype[float64]]

Parameters

dtfloat

Deprecated and ignored.

Returns

forcesNDArray[np.float64]

Body force array

concentration_dg()

Return the discontinuous concentration field on P1DG nodes. Returns ——-

concentrationRealDeviceArray

Concentration array

Return type:

RealDeviceArray

concentration_dg_grad()

Return the continuous concentration gradient on P1 nodes. Returns ——-

gradientRealDeviceArray

Concentration gradient array

Return type:

RealDeviceArray

coordinates()

Return node coordinates as a (n_nodes, 3) array. Returns ——-

coordinatesNDArray[np.float64]

Node coordinates

Return type:

ndarray[Any, dtype[float64]]

coordinates_fields()

Return spatial coordinates for every DOF in the solution vector, shape (n_sol_dofs, dim). Returns ——-

coordsNDArray[np.float64]

Coordinate array

Return type:

ndarray[Any, dtype[float64]]

density()

Return the density field. Returns ——-

densityRealDeviceArray

Density array

Return type:

RealDeviceArray

dimension()

Return the spatial dimension (2 or 3). Returns ——-

dimint

Spatial dimension

Return type:

int

element_size()

Return the representative size of each element, shape (n_elements,). Returns ——-

hRealDeviceArray

Element size array

Return type:

RealDeviceArray

element_tags()

Return element tags. Returns ——-

tagsNDArray[np.int32]

Element tags

Return type:

ndarray[Any, dtype[int32]]

elements()

Return the element connectivity. Returns ——-

elementsNDArray[np.int32]

Element array, shape (n_elements, dim + 1)

Return type:

ndarray[Any, dtype[int32]]

enable_stability(enable_pspg, enable_supg)

Enable or disable PSPG and SUPG stabilisation terms.

Parameters

enable_pspgint

Enable PSPG (pressure-stabilising Petrov-Galerkin): 1=yes

enable_supgint

Enable SUPG (streamline-upwind Petrov-Galerkin): 1=yes

field_indices(ifield)

Return solution-vector indices for one field.

Return type:

ndarray[Any, dtype[int32]]

Parameters

ifieldint

Field number

Returns

idxNDArray[np.int32]

Solution-vector indices

fields_gradient()

Return the continuous gradient of each field at P1 nodes. Shape: (ndof_p1, n_fields, dim). Returns ——-

gradNDArray[np.float64]

Gradient array

Return type:

ndarray[Any, dtype[float64]]

full_implicit_euler(dt, tol=1e-06, reduced_gravity=0, stab_param=0.0, itermax=100)

Advance by one time step using a fully implicit Newton scheme.

Parameters

dtfloat

Time step

tolfloat

Newton convergence tolerance

reduced_gravityint

Use reduced gravity formulation

stab_paramfloat

Additional stabilization parameter

itermaxint

Maximum Newton iterations

g()

Return the gravity vector the problem was constructed with. The coupling closure and the Patankar datum both need it (the datum carries g*dt), and a scheme that builds them has no other way to ask. Returns ——-

gNDArray[np.float64]

Gravity vector, shape (dim,)

Return type:

ndarray[Any, dtype[float64]]

get_bodies_csr_force()

Return fluid/body force at body CSR integration points. Returns ——-

forcesRealDeviceArray

CSR body force array

Return type:

RealDeviceArray

get_bodies_csr_force_penalty()

Return the PENALTY-REACTION part of the per-entry body force. get_bodies_csr_force() returns the total, which also carries the -int_G grad p buoyancy – the traction-like part, which samples the pressure inside the cut band where it is noisiest. This one is the variationally consistent route: the momentum exchange int beta(u-g) plus the trace band, read from the assembled equations. Returns ——-

forceRealDeviceArray

Per-entry penalty reaction, shape (nnz, dim)

Return type:

RealDeviceArray

get_default_export()

Return the default fields for write_mig as {name: (data, element)}. Mirrors Python FluidProblem.get_default_export(). Returns ——-

fieldsDict[str,Tuple[NDArray[np.float64], str]]

Default export fields

Return type:

Dict[str, Tuple[ndarray[Any, dtype[float64]], str]]

get_density_element()

Return the element type string for the density field. Returns ——-

elementstr

Element type string

Return type:

str

get_forces_on_bodies()

Return the fluid force applied on each body, shape (n_bodies, dim). Returns ——-

forcesNDArray[np.float64]

Body force array

Return type:

ndarray[Any, dtype[float64]]

get_mapping(etype)

Return the mapping associated with an element type.

Return type:

ndarray[Any, dtype[int32]]

Parameters

etypestr

Element type, e.g. P1

Returns

mappingNDArray[np.int32]

Element mapping

get_p1_element()

Return the default P1 element name. Returns ——-

elementstr

Element name

Return type:

str

get_p1_mapping()

Return the P1 DOF index for each mesh node. Returns ——-

mappingNDArray[np.int32]

Integer array of length n_nodes

Return type:

ndarray[Any, dtype[int32]]

get_pressure()

Return the pressure field as a (n_p_nodes,) array. Returns ——-

pressureNDArray[np.float64]

Pressure array

Return type:

ndarray[Any, dtype[float64]]

get_pressure_element()

Return the element type string for pressure degrees of freedom. Returns ——-

elementstr

Element type string

Return type:

str

get_velocity_element()

Return the element type string for velocity degrees of freedom. Returns ——-

elementstr

Element type string

Return type:

str

implicit_euler(dt, check_residual_norm=-1, reduced_gravity=False, stab_param=0.0)

Advance the solution one time step with the implicit Euler scheme.

Parameters

dtfloat

Time-step size

check_residual_normfloat

If > 0, throw if residual norm exceeds this value after solve

reduced_gravitybool

Use reduced-gravity formulation

stab_paramfloat

If non-zero, use as pressure-Laplacian stabilisation coefficient instead of PSPG/SUPG

interpolate(solution=None, velocity=None, velocity_x=None, velocity_y=None, velocity_z=None, pressure=None)

Assign solution, velocity, velocity components, or pressure from arrays or callbacks.

Parameters

solutionUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Full solution values, or callback on coordinates_fields()

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

Velocity values, or callback on velocity DOF coordinates

velocity_xUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

X velocity values/callback

velocity_yUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Y velocity values/callback

velocity_zUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Z velocity values/callback

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure values/callback

local_boundary_force_by_contribution(tag)

Return per-edge boundary forces split into pressure and viscous contributions. Shape: (n_edges, 2*dim).

Return type:

ndarray[Any, dtype[float64]]

Parameters

tagstr

Boundary name

Returns

forcesNDArray[np.float64]

Per-edge force array, shape (n_edges, 2*dim)

local_size()

Return the number of local degrees of freedom per element. Returns ——-

nint

Local system size

Return type:

int

mesh_boundaries()

Return all mesh boundaries as {name: edge_nodes}. Mirrors Python FluidProblem.mesh_boundaries(). Returns ——-

boundariesDict[str,NDArray[np.int32]]

Map from boundary name to node index array, shape (bsize, dim)

Return type:

Dict[str, ndarray[Any, dtype[int32]]]

mesh_velocity()

Compatibility alias for mesh_velocity_array(). Returns ——-

velRealDeviceArray

Mesh velocity array, shape (n_nodes, 3)

Return type:

RealDeviceArray

mesh_velocity_array()

Return the mesh velocity, in device memory, shape (n_nodes, dim). Returns ——-

velRealDeviceArray

Mesh velocity device array

Return type:

RealDeviceArray

n_elements()

Return the number of mesh elements. Returns ——-

nint

Number of elements

Return type:

int

n_fields()

Return the number of solution fields. Returns ——-

nint

Number of fields

Return type:

int

n_nodes()

Return the number of mesh nodes. Returns ——-

nint

Number of nodes

Return type:

int

node_volume()

Return nodal control volumes. Returns ——-

volumeRealDeviceArray

Node volumes

Return type:

RealDeviceArray

p_jump()

Return the pressure jump field, on P1DG nodes. Returns ——-

p_jumpRealDeviceArray

Pressure jump array

Return type:

RealDeviceArray

particle_volume_intersected()

Return, per body, how much of its volume lies inside the mesh: the sum of its contributions’ volumes. Bodies outside the mesh give zero. Returns ——-

volumeNDArray[np.float64]

Intersected volume per body, shape (n_bodies,)

Return type:

ndarray[Any, dtype[float64]]

porosity()

Compatibility alias for porosity_array(). Returns ——-

porRealDeviceArray

Porosity array, shape (n_nodes, 1)

Return type:

RealDeviceArray

porosity_array()

Return the porosity, in device memory, shape (n_nodes,). Returns ——-

porRealDeviceArray

Porosity device array

Return type:

RealDeviceArray

pressure()

Compatibility alias for get_pressure(). Returns ——-

pressureNDArray[np.float64]

Pressure array

Return type:

ndarray[Any, dtype[float64]]

read_mig(odir, t=-1.0, iteration=-2147483648)

Read a fluid state written by write_mig().

Parameters

odirstr

Output directory

tfloat

Time to read. If omitted, iteration must be provided.

iterationint

Iteration index, or -1 for last.

reset_timings()

Zero the timing counters, to skip a warm-up phase.

set_concentration_cg(concentration)

Set the concentration field from continuous nodal values.

Parameters

concentrationNDArray[np.float64]

Concentration array, shape (n_nodes,)

set_coordinates(x)

Move the mesh nodes.

Parameters

xNDArray[np.float64]

Node coordinates, shape (n_nodes, 3)

set_coupling_closure(csr_m, csr_mgsp, csr_div_m, csr_beta, csr_h, csr_velocity, csr_contact, csr_eta=None, csr_eta_p=None, csr_gamma_t=None, csr_stab_w=None)

TIER B – the model on the current geometry: the unfitted map and its divergence, the penalty, the trace band, the frozen velocity datum and the contact force. Requires a geometry of the current generation.

Parameters

csr_mNDArray[np.float64]

Unfitted map m, shape (nnz, dim)

csr_mgspNDArray[np.float64]

m . grad(u_p), shape (nnz, dim)

csr_div_mNDArray[np.float64]

div m, shape (nnz,)

csr_betaNDArray[np.float64]

Penalty coefficient, shape (nnz,)

csr_hNDArray[np.float64]

Trace length, shape (nnz,)

csr_velocityNDArray[np.float64]

Frozen body-velocity datum, shape (nnz, dim)

csr_contactNDArray[np.float64]

Contact force per entry, shape (nnz, dim)

csr_etaNDArray[np.float64]

DIMENSIONLESS interior-fill coefficient per entry; absent falls back to the global set_eta

csr_eta_pNDArray[np.float64]

DIMENSIONLESS pressure-diffusivity coefficient per entry; absent falls back to the global set_eta_p

csr_gamma_tNDArray[np.float64]

DIMENSIONLESS trace coefficient per entry; absent means the global set_trace_param value. Switch the band off with 0.0 – never by sending csr_h to infinity, which is what a 1e150 sentinel used to mean and put its square one decade from overflow.

csr_stab_wNDArray[np.float64]

Fraction of the coupling coefficient that is PHYSICAL DRAG rather than numerical penalty, per entry, in [0, 1]. It is the weight with which tricks 2 and 3 enter the SUPG/PSPG residual: 1 puts the whole coupling there, as the legacy mixture kernel did with its drag; 0 keeps it out, which a Babuska penalty requires because betap*pspg = gamma_B/4 is mesh-independent and would never refine away. The producer sets it to the drag fraction of its own blend, so a sub-cell grain gets 1 and a resolved one gets 0. Absent falls back to the global drag_in_stab flag.

set_coupling_datum(csr_velocity, csr_dup_ds, csr_dup_dds)

TIER C – the affine datum a scheme supplies so the body velocity responds to the fluid INSIDE the solve: up_eff = velocity + dup_ds . s + dup_dds . ds. Optional; with no datum set for the current generation the closure’s frozen velocity stands, which is what iterate_iqn wants. The velocity given here REPLACES the closure’s (same slot, better value).

Parameters

csr_velocityNDArray[np.float64]

Constant part up0, shape (nnz, dim)

csr_dup_dsNDArray[np.float64]

d u_p / d s_F, shape (nnz, dim, dim + 1)

csr_dup_ddsNDArray[np.float64]

d u_p / d (ds_F/dx), shape (nnz, dim, dim + 1, dim)

set_coupling_geometry(density, volume, csr_ptr, csr_element_id, csr_xi, csr_surface, csr_weight, csr_d, csr_cut=None)

TIER A – the overlap geometry, and the per-body constants that go with it. Pure of state: no velocity, no fluid solution, no dt. Setting it BUMPS THE COUPLING GENERATION and invalidates any closure and datum, so a stale closure sized to a previous geometry, or last step’s datum left behind after the bodies moved, are refused rather than silently used. Setting a geometry and NOT completing it with a closure is refused at solve time too – it used to be a silent no-op that ran the operator against the previous commit. A held bed can call this once for the whole run; IQN once per step rather than once per iteration.

Parameters

densityNDArray[np.float64]

Per-body density, shape (n_bodies,)

volumeNDArray[np.float64]

Per-body volume, shape (n_bodies,)

csr_ptrNDArray[np.int32]

Body -> entry offsets, shape (n_bodies + 1,)

csr_element_idNDArray[np.int32]

Element of each entry, shape (nnz,)

csr_xiNDArray[np.float64]

Reference coordinates, shape (nnz, dim)

csr_surfaceNDArray[np.float64]

Interface measure per entry, shape (nnz,)

csr_weightNDArray[np.float64]

Integration weight per entry, shape (nnz,)

csr_dNDArray[np.float64]

Signed distance per entry, shape (nnz,)

csr_cutNDArray[np.int32]

TOPOLOGICAL per-entry cut flag (1 = the element is only PARTIALLY covered by this body). For a convex body this is exactly nin < dim+1, with no threshold anywhere. Absent falls back to the kernel’s old volume-fraction inference.

set_eta(eta)

Set the interior grain-Laplacian fill coefficient of the unfitted kernel.

Parameters

etafloat

Dimensionless fill coefficient, applied as eta * mu

set_eta_p(eta_p)

Set the covered-pressure Laplacian coefficient of the unfitted kernel. Applied as eta_p * taup/rho on the covered (P,P) block. 1e-3 is the production default; the legacy-retrieval configuration passes 0.

Parameters

eta_pfloat

Dimensionless covered-pressure coefficient

set_g(g)

Set the gravity vector (body force per unit mass). The constructor value is otherwise immutable, and g() hands back a COPY, so a scheme that has to STEER the body force – a periodic cell driven to a target flow rate, where the driver is the unknown and the velocity is imposed – has no way to write it back. This is that way.

Parameters

gNDArray[np.float64]

Gravity vector, shape (dim,)

set_mean_pressure(p)

Constrain the mean pressure to a fixed value.

Parameters

pfloat

Target mean pressure value

set_mesh(nodes, elements, element_tags, boundaries, boundary_tags, boundary_names, periodic=None)

Set the mesh from arrays.

Parameters

nodesNDArray[np.float64]

Node coordinates, shape (n_nodes, 3)

elementsNDArray[np.int32]

Element connectivity, shape (n_elements, dim+1)

element_tagsNDArray[np.int32]

Element tags, shape (n_elements,)

boundariesNDArray[np.int32]

Boundary edge nodes, shape (n_bnd, dim)

boundary_tagsNDArray[np.int32]

Boundary tags, shape (n_bnd,)

boundary_namesList[str]

Physical group names

periodicNDArray[np.int32]

Parent DOF index for each node (n_nodes,), or None

set_open_boundary(tag, velocity=None, pressure=None, viscous_flag=True, porous=True)

Set an open (type=1, viscous_flag=1) boundary condition.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

viscous_flagbool

Compute viscous term at the boundary

porousbool

Weight this boundary’s data by the CURRENT porosity. True imposes a SUPERFICIAL quantity, as the mixture kernel always did (f0(P) = c*unext): a bed sitting on an inlet then admits only c times the nominal velocity. False imposes the datum as given. Measured on 2d_bubbling_fluidized_bed, where the bed sits on the inlet at c = 0.5982: imposing the datum unweighted made the whole solution 1.683x too large – velocity, pressure and grain force alike, and in grain-free regions too, because an inflow-flux error is global. False is what a traction outlet wants: the gradient-form momentum reduces to the imposed traction only at c = 1 (3d_darcy_permeability, c = 0.478, showed p_out = 0.0354 against an imposed 0 and a 10.1% grain-force overshoot at c < 1).

set_particles_body(xp, rp, density, velocity=None, omega=None, contact_forces=None, discretisation='overlap', use_voidage=True)

Set body coupling data from particle centers/radii and particle fields.

Parameters

xpNDArray[np.float64]

Particle positions, shape (n_particles, dim)

rpNDArray[np.float64]

Particle radii, shape (n_particles,)

densityNDArray[np.float64]

Particle densities, shape (n_particles,)

velocityNDArray[np.float64]

Particle velocities, shape (n_particles, dim)

omegaNDArray[np.float64]

Particle angular velocities, shape (n_particles, dim)

contact_forcesNDArray[np.float64]

Particle contact forces, shape (n_particles, dim)

discretisationstr

overlap in 2D or centroid otherwise

use_voidagebool

Include voidage in drag coefficient

set_particles_body_csr(xp, rp, density, csr_ptr, csr_element_id, csr_xi, csr_surface, csr_weight, velocity=None, omega=None, contact_forces=None, use_voidage=True)

Set body coupling data from particle fields and a CALLER-SUPPLIED quadrature CSR – the mixture-path twin of set_unfitted_bodies. The overlap geometry is produced outside the core (overlap, through migflow.overlap_quadrature); the solver only copies. An entry is a quadrature point: csr_weight integrates, csr_surface is the pair-total overlap duplicated on each entry (cut-band detection only).

Parameters

xpNDArray[np.float64]

Particle positions, shape (n_particles, dim)

rpNDArray[np.float64]

Particle radii, shape (n_particles,)

densityNDArray[np.float64]

Particle densities, shape (n_particles,)

csr_ptrNDArray[np.int32]

CSR row pointer, shape (n_particles + 1,)

csr_element_idNDArray[np.int32]

Element of each entry, shape (nnz,)

csr_xiNDArray[np.float64]

Reference coordinates of each point, shape (nnz, dim)

csr_surfaceNDArray[np.float64]

Pair-total overlap area/volume, shape (nnz,)

csr_weightNDArray[np.float64]

Integration weight of each point, shape (nnz,)

velocityNDArray[np.float64]

Particle velocities, shape (n_particles, dim)

omegaNDArray[np.float64]

Particle angular velocities, shape (n_particles, dim)

contact_forcesNDArray[np.float64]

Particle contact forces, shape (n_particles, dim)

use_voidagebool

Include voidage in drag coefficient

set_polygon_bodies(vertices, density, volume, gamma, velocity=None, contact_forces=None)

Set body coupling data from convex polygons rather than discs. The drag coefficient is taken from the caller instead of being derived from a radius, which a polygon does not have.

Parameters

verticesNDArray[np.float64]

Polygon vertices, shape (n_bodies, n_vertices, dim)

densityNDArray[np.float64]

Body densities, shape (n_bodies,)

volumeNDArray[np.float64]

Body volumes, shape (n_bodies,)

gammaNDArray[np.float64]

Drag coefficient of each body, shape (n_bodies,)

velocityNDArray[np.float64]

Body velocities, shape (n_bodies, dim)

contact_forcesNDArray[np.float64]

Body contact forces, shape (n_bodies, dim)

set_pressure(pres)

Set the pressure field from a (n_p_nodes,) array.

Parameters

presNDArray[np.float64]

Pressure array

set_stabilisation_complement(pspg=True, supg=True)

Choose which halves of f_volume’s stabilisation image the unfitted body kernel subtracts over the element-grain overlap. PSPG is needed for equal-order P1/P1 and its complement is uncontroversial. The SUPG complement is not: f_volume’s SUPG image carries the convective term, and over the covered region the advecting field is the extended velocity INSIDE the grain, which is fictitious. Subtracting a convective stabilisation built from it is a modelling choice, invisible at Stokes – which is exactly where this formulation agrees with cutfem.

Parameters

pspgbool

Complement the PSPG image

supgbool

Complement the SUPG image

set_strong_boundary(tag, velocity=None, velocity_x=None, velocity_y=None, velocity_z=None)

Set a strong (Dirichlet) boundary condition. velocity imposes every component at once; velocity_x/y/z impose one, leaving the others free, which is what a slip wall needs.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

velocity_xUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

X velocity value/callback, or None

velocity_yUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Y velocity value/callback, or None

velocity_zUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Z velocity value/callback, or None

set_symmetry_boundary(tag, pressure=None)

Set a symmetry boundary: only the normal flux is imposed. This is a weak boundary with no velocity and viscous_flag left unspecified, which with no velocity imposed resolves to a traction-free surface.

Parameters

tagstr

Boundary physical-group name

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

set_trace_param(gamma_t)

Set the trace-stabilisation coefficient of the unfitted kernel.

Parameters

gamma_tfloat

Dimensionless trace coefficient

set_two_fluid_properties(rho, mu, sigma=0.0)

Derive the density and viscosity fields from the concentration, and optionally the capillary force. The density and viscosity elements must both be P1DG, like the concentration. This helper owns bulk_force: it overwrites it on every call when sigma is non-zero. Add any other bulk force contribution afterwards.

Parameters

rhoList[float]

Density of each of the two fluids

muList[float]

Dynamic viscosity of each of the two fluids

sigmafloat

Surface tension coefficient; 0 leaves bulk_force alone

set_unfitted(unfitted)

Select the unfitted (boundary-free) body model. With it enabled the no-slip datum on a grain boundary is imposed purely through volume integrals over the element-grain overlap – no facet loop, no interface integral – and the nodal porosity no longer weights the volume operators: the geometry enters as the sharp complement int_K (.) - int_{K n G} (.).

Parameters

unfittedbool

Whether to use the unfitted kernel

set_unfitted_bodies(density, volume, csr_ptr, csr_element_id, csr_xi, csr_surface, csr_weight, csr_velocity, csr_contact, csr_m, csr_mgsp, csr_div_m, csr_d, csr_h, csr_beta, csr_dup_ds=None, csr_dup_dds=None, csr_eta=None, csr_eta_p=None, csr_gamma_t=None, csr_cut=None, csr_stab_w=None)

Hand the unfitted overlap table to the solver as a flat CSR. An entry is a QUADRATURE POINT, not a (body, element) pair: csr_element_id repeats within a body, csr_xi is the point in reference coordinates, and csr_weight is its integration weight. csr_surface stays the pair total, duplicated on every entry, and is used only for the geometric cut test and as a diagnostic. The geometry is produced entirely caller-side; the solver only copies. Arguments follow the field order of the internal FluidProblemBodies.

Parameters

densityNDArray[np.float64]

Per-body density, shape (n_bodies,)

volumeNDArray[np.float64]

Per-body volume, shape (n_bodies,)

csr_ptrNDArray[np.int32]

CSR row pointer, shape (n_bodies + 1,)

csr_element_idNDArray[np.int32]

Element of each entry, shape (nnz,)

csr_xiNDArray[np.float64]

Reference coordinates of each point, shape (nnz, dim)

csr_surfaceNDArray[np.float64]

Pair-total overlap area/volume, shape (nnz,)

csr_weightNDArray[np.float64]

Integration weight of each point, shape (nnz,)

csr_velocityNDArray[np.float64]

Imposed velocity datum u_p, shape (nnz, dim). When the sensitivity arrays below are supplied this is the CONSTANT PART of an AFFINE datum: the kernel imposes up_eff = csr_velocity + dup_ds . s + dup_dds . ds, evaluated at each quadrature point’s own state, in every residual and force line. The caller must therefore fold its linearisation reference into csr_velocity analytically – csr_velocity = up_pred - dup_ds . s_ref - dup_dds . ds_ref – so that up_eff reproduces up_pred at the reference state. Left null the datum is frozen and csr_velocity is u_p as-is.

csr_contactNDArray[np.float64]

Contact force at each point, shape (nnz, dim)

csr_mNDArray[np.float64]

Direction field (x-c)/r, shape (nnz, dim)

csr_mgspNDArray[np.float64]

m . grad(u_p) per velocity component, shape (nnz, dim)

csr_div_mNDArray[np.float64]

Divergence of m, shape (nnz,)

csr_dNDArray[np.float64]

Signed distance, negative inside the grain, shape (nnz,)

csr_hNDArray[np.float64]

Penalty length, constant per grain, shape (nnz,)

csr_betaNDArray[np.float64]

Babuska penalty density gamma/h^2, shape (nnz,)

csr_dup_dsNDArray[np.float64]

d u_p,i / d s_F, shape (nnz, dim, dim + 1)

csr_dup_ddsNDArray[np.float64]

d u_p,i / d (ds_F/dx_a), shape (nnz, dim, dim + 1, dim)

csr_etaNDArray[np.float64]

DIMENSIONLESS interior-fill coefficient per entry; absent falls back to the global set_eta

csr_eta_pNDArray[np.float64]

DIMENSIONLESS pressure-diffusivity coefficient per entry; absent falls back to the global set_eta_p

csr_gamma_tNDArray[np.float64]

DIMENSIONLESS trace coefficient per entry; absent means the global set_trace_param value. betat = mu*gamma_t/h^2, so the band goes off with gamma_t = 0.0 – exactly, and with h left a real length – rather than by sending h to infinity.

csr_cutNDArray[np.int32]

TOPOLOGICAL per-entry cut flag (1 = the element is only PARTIALLY covered by this body). Absent falls back to the kernel’s old volume-fraction inference, surface < 0.99*|K|. The two sensitivity arrays are how a SEMI-IMPLICIT coupling is expressed without a kernel branch. Left null (the default) the imposed velocity is a frozen Dirichlet datum, every Jacobian contribution below is multiplied by zero, and the assembly is bit-for-bit what it was before they existed – that is the property the explicit and the iterating schemes rely on, where dv = 0 is EXACT for the statement being linearised. Filled, they let the caller declare that u_p depends on the fluid state at the same solve, which is what the legacy mixture kernel did internally with its up_pred: the old du_pred_du is csr_dup_ds[i][U+i] and the old du_pred_dp is csr_dup_dds[i][P][i]. The field index F runs over the NF solution fields, U + 0 .. U + dim - 1 then P. SUPPLYING A SENSITIVITY THAT DOES NOT MATCH THE VELOCITY BESIDE IT IS SILENT AND CATASTROPHIC – implicit_euler is a single Newton step, so a mismatch corrupts the answer outright rather than merely slowing convergence. Measured on the no-slip MMS: a mismatched pair turned a 63% force error into 131512%, and attaching a sensitivity implying a mobility to a residual holding the velocity fixed gave 69479% against 245%. Derive all three from one closure (see body_coupling.get_patankar_pred) and check them with body_coupling.check_sensitivity before trusting a number.

csr_stab_wNDArray[np.float64]

Fraction of the coupling coefficient that is PHYSICAL DRAG rather than numerical penalty, per entry, in [0, 1]. The weight with which tricks 2 and 3 enter the SUPG/PSPG residual; see set_coupling_closure. Absent falls back to the global drag_in_stab.

set_velocity(vel)

Set the velocity field from a (n_vel_nodes, dim) array.

Parameters

velNDArray[np.float64]

Velocity array, shape (n_vel_nodes, dim)

set_wall_boundary(tag, velocity=None, pressure=None, viscous_flag=-1, porous=True)

Set a wall-type (type=0) boundary condition.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

viscous_flagint

Compute viscous term at the boundary: 1 to force it on, 0 to force it off, -1 (default) to compute it only when a velocity is imposed

porousbool

Weight this boundary’s data by the CURRENT porosity. True imposes a SUPERFICIAL quantity, as the mixture kernel always did (f0(P) = c*unext): a bed sitting on an inlet then admits only c times the nominal velocity. False imposes the datum as given. Measured on 2d_bubbling_fluidized_bed, where the bed sits on the inlet at c = 0.5982: imposing the datum unweighted made the whole solution 1.683x too large – velocity, pressure and grain force alike, and in grain-free regions too, because an inflow-flux error is global. False is what a traction outlet wants: the gradient-form momentum reduces to the imposed traction only at c = 1 (3d_darcy_permeability, c = 0.478, showed p_out = 0.0354 against an imposed 0 and a 10.1% grain-force overshoot at c < 1).

set_weak_boundary(tag, velocity=None, pressure=None, concentration=None, viscous_flag=-1, porous=True)

Set a weak (type=0) boundary condition.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

concentrationUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Concentration/alpha value for two-fluid problems, or None

viscous_flagint

Compute viscous term at the boundary: 1 to force it on, 0 to force it off, -1 (default) to enable it only when a velocity or a pressure is imposed

porousbool

Weight this boundary’s data by the CURRENT porosity; see set_open_boundary.

set_weak_boundary_nodal_values(tag, values)

Replace the nodal values of a weak boundary already set by set_weak_boundary(). Unlike set_weak_boundary() this does not add a boundary, so it can be called at every time step.

Parameters

tagstr

Boundary physical-group name

valuesNDArray[np.float64]

One value per field and per boundary node, shape (n_edges, n_closure_dofs, n_values)

solution()

Compatibility alias for solution_array(). Returns ——-

solRealDeviceArray

Solution array of length n_dof

Return type:

RealDeviceArray

solution_array()

Return the full solution, in device memory, as an array of length n_dof. Returns ——-

solRealDeviceArray

Solution device array of length n_dof

Return type:

RealDeviceArray

solution_at_coordinates(x)

Interpolate the solution at arbitrary coordinates.

Return type:

ndarray[Any, dtype[float64]]

Parameters

xNDArray[np.float64]

Coordinates, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

solution_at_reference(eid, xi)

Interpolate the solution at points given by element id and REFERENCE coordinates. No mesh search: use this whenever the caller already knows which element each point is in – an overlap CSR does, and going in by physical position instead costs a locate_point_in_mesh per point and is ambiguous on a face. Evaluated with the field’s OWN shape functions, so it is correct at any element order (a hand-rolled P1 formula in the caller is not).

Return type:

ndarray[Any, dtype[float64]]

Parameters

eidNDArray[np.int32]

Element index per point, shape (n,)

xiNDArray[np.float64]

Reference coordinates in that element, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

solution_grad_at_coordinates(x)

Solution and gradient at arbitrary physical coordinates: locates each point, then evaluates as solution_grad_at_reference does.

Return type:

Tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]]]

Parameters

xNDArray[np.float64]

Coordinates, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

gradientNDArray[np.float64]

Solution gradient, shape (n, n_fields, dim)

solution_grad_at_reference(eid, xi)

Solution AND its gradient at element id + reference coordinates, in one pass over the shape functions.

Return type:

Tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]]]

Parameters

eidNDArray[np.int32]

Element index per point, shape (n,)

xiNDArray[np.float64]

Reference coordinates in that element, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

gradientNDArray[np.float64]

Solution gradient, shape (n, n_fields, dim)

timings()

Return the wall-clock breakdown of implicit_euler, in seconds. Slots, in order: setup, assemble, scatter, factorise, solve, bodies, n_steps. Each mark synchronises the device first, because the assembly kernels are asynchronous and an unsynced timestamp would charge their cost to whichever phase blocks next. That serialisation makes the TOTAL a slight over-estimate; the split is what it is for. Returns ——-

timingsNDArray[np.float64]

Seven values: the six phases and the step count

Return type:

ndarray[Any, dtype[float64]]

u_solid()

Return the solid velocity field at P1 nodes, shape (n_p1_nodes, dim). Returns ——-

usNDArray[np.float64]

Solid velocity array

Return type:

ndarray[Any, dtype[float64]]

update_node_volume()

Rebuild the cached mesh metrics, including the nodal control volumes, after the mesh has been deformed in place through coordinates().

velocity()

Return the velocity field as a (n_vel_nodes, dim) array. Returns ——-

velocityNDArray[np.float64]

Velocity array

Return type:

ndarray[Any, dtype[float64]]

velocity_dof_coordinates()

Return spatial coordinates of each velocity DOF, shape (n_vel_nodes, dim). Returns ——-

coordsNDArray[np.float64]

Coordinate array

Return type:

ndarray[Any, dtype[float64]]

velocity_index()

Return solution-vector indices for velocity DOFs, shape (n_vel_dofs, dim). velocity_index()[i, d] is the solution index of component d at velocity DOF i. Returns ——-

idxNDArray[np.int32]

Index array

Return type:

ndarray[Any, dtype[int32]]

viscosity()

Return the dynamic-viscosity field. Returns ——-

viscosityRealDeviceArray

Viscosity array

Return type:

RealDeviceArray

write_mig(output_dir, t, fields={})

Write output files for post-visualisation. Mirrors Python FluidProblem.write_mig().

Parameters

output_dirstr

Output directory

tfloat

Computational time

fieldsDict[str,Tuple[NDArray[np.float64], str]]

Fields to write as {name: (data, element)}. Uses get_default_export() if empty.

class migflow.fluid.FluidProblem3(g=None, mu=0.0, rho=0.0, coeff_stab=0.01, volume_drag=0.0, quadratic_drag=0.0, drag_in_stab=False, drag_coefficient_factor=1.0, temporal=True, advection=True, flag_div_us=True, p2p1=False, full_implicit=False, model_b=False, density_element='', viscosity_element='', solver=None, solver_options='')

Bases: object

MigFlow incompressible fluid solver (2-D or 3-D). Exposes a high-level C++/Python API. Call set boundary conditions, then iterate with implicit_euler().

Construct a FluidProblem. Mirrors Python FluidProblem.__init__().

Parameters

gUnion[None, NDArray[np.float64], List[float]]

Gravity vector, length dim. Default: zero vector.

mufloat

Dynamic viscosity, the initial uniform value of the viscosity field

rhofloat

Density, the initial uniform value of the density field

coeff_stabfloat

Stabilisation coefficient

volume_dragfloat

Volume drag coefficient

quadratic_dragfloat

Quadratic drag coefficient

drag_in_stabbool

Include the coupling terms in the SUPG/PSPG residual. In the mixture model this is the drag term. In the UNFITTED model it is the Babuska penalty (trick 2) and the reconstructed-trace band (trick 3): the stabilisation block is applied after them, so the flag selects between the Galerkin-only snapshot and the full residual. True reproduces what the legacy mixture kernel did, and is required under model=’legacy’, where beta = gamma_D/(mu V_b) is a physical drag coefficient – excluding it leaves the pressure row blind to the drag and a fluidised bed collapses (2d_bubbling_fluidized_bed, drag/weight saturating at 0.61). False restores the earlier unfitted-only rule, under which betap*pspg = gamma_B/4 is mesh-independent so a stabilised Babuska penalty never refines away – keep it False when measuring convergence orders of the pure unfitted model. DEFAULT IS FALSE, and the two regimes genuinely disagree, so this is a modelling choice the case must make rather than inherit. Measured 2026-08-21 on one build: flag 2d_stokes_noslip_mms force error 2d_bubbling_fluidized_bed 0 244% -> 63% -> 16% (converging) bed COLLAPSES 4.67 -> 2.24 cm 1 1023% -> 1368% -> 1643% (diverging) bed FLUIDISES 4.67 -> 13.2 cm A sub-cell grain carries a physical drag that belongs in the pressure row; a resolved grain carries a mesh penalty that must stay out of it. The blend beta = s*gamma_B + (1-s)*gamma_D already interpolates between those two regimes per entry, so the principled successor to this global flag is a per-entry weight of (1-s) – not yet implemented.

drag_coefficient_factorfloat

Factor multiplying the drag coefficient

temporalbool

Enable temporal (d/dt) term

advectionbool

Enable advective terms (Navier-Stokes vs Stokes)

flag_div_usbool

Enable div(u_solid) term

p2p1bool

Use P2P1 (Taylor-Hood) elements instead of stabilised P1P1

full_implicitbool

Use fully implicit nonlinear scheme

model_bbool

Enable model B

density_elementstr

Density discretisation element (empty = default)

viscosity_elementstr

Viscosity discretisation element (empty = default)

solverUnion[None, str, ILinearSystem]

Optional plug-in linear solver object or solver name: petsc, pardiso, cudss, amgx, or default.

solver_optionsstr

Options forwarded to the named solver: a PETSc options string for petsc/default, an AmgX JSON configuration or config file path for amgx.

Pointer

alias of LP__Structure

adapt_mesh(nodes, elements, element_tags, boundaries, boundary_tags, boundary_names, periodic=None)

Adapt the mesh, projecting the current solution onto the new mesh.

Parameters

nodesNDArray[np.float64]

Node coordinates, shape (n_nodes, 3)

elementsNDArray[np.int32]

Element connectivity, shape (n_elements, dim+1)

element_tagsNDArray[np.int32]

Element tags, shape (n_elements,)

boundariesNDArray[np.int32]

Boundary edge nodes, shape (n_bnd, dim)

boundary_tagsNDArray[np.int32]

Boundary tags, shape (n_bnd,)

boundary_namesList[str]

Physical group names

periodicNDArray[np.int32]

Parent DOF index for each node (n_nodes,), or None

add_constraint(dofs, weights, rhs)

Add the linear constraint sum_i weights[i] * solution[dofs[i]] = rhs.

Parameters

dofsNDArray[np.int32]

Solution-vector indices

weightsNDArray[np.float64]

One weight per index

rhsfloat

Right-hand side

advance_concentration(dt)

Advance the concentration field by one time step (two-fluid problems).

Parameters

dtfloat

Time step

assemble_local_system(dt, localv, localm)

Assemble the element-local residual and Jacobian at the CURRENT solution, without solving. DIAGNOSTIC ENTRY POINT. It exists so a finite-difference check can ask whether the assembled Jacobian really is the derivative of the assembled residual: perturb one dof, reassemble, and compare the difference quotient against the matrix column. Without it the two can only be compared through a solve, which conflates a wrong derivative with an ill-conditioned one. Call it AFTER at least one implicit_euler, which is what sets the mesh metrics, the boundary conditions and the old-solution snapshot; this method deliberately refreshes none of those, so sold stays fixed while the current solution is perturbed.

Parameters

dtfloat

Time step

localvRealDeviceArray

Output local right-hand side, (n_elements, local_size)

localmRealDeviceArray

Output local matrix, (n_elements, local_size, local_size)

boundary_force_by_contributions(tag)

Return the total forces at a named boundary split into pressure and viscous. Shape: (2*dim,) = [p_x,..,p_dim, v_x,..,v_dim].

Return type:

ndarray[Any, dtype[float64]]

Parameters

tagstr

Boundary name

Returns

forcesNDArray[np.float64]

Force array, size 2*dim

boundary_forces()

Return the accumulated weak-boundary force array. Shape: (n_boundary_edges, dim). Returns ——-

forcesRealDeviceArray

Per-boundary-edge force array

Return type:

RealDeviceArray

bulk_force()

Return the bulk force field on the P1 nodes, shape (ndof_p1, dim). A body force per unit volume, added to the momentum equation. Note that compute_capillary_force() OVERWRITES it whenever sigma is nonzero, so a caller setting it directly must either leave surface tension off or set it again after each call. Returns ——-

bulk_forceRealDeviceArray

Bulk force field

Return type:

RealDeviceArray

clear_constraints()

Drop every constraint added by add_constraint(). A caller that rebuilds its constraints at each time step must call this first, otherwise they accumulate.

compute_node_force(dt=-1)

Compatibility alias for get_forces_on_bodies().

Return type:

ndarray[Any, dtype[float64]]

Parameters

dtfloat

Deprecated and ignored.

Returns

forcesNDArray[np.float64]

Body force array

concentration_dg()

Return the discontinuous concentration field on P1DG nodes. Returns ——-

concentrationRealDeviceArray

Concentration array

Return type:

RealDeviceArray

concentration_dg_grad()

Return the continuous concentration gradient on P1 nodes. Returns ——-

gradientRealDeviceArray

Concentration gradient array

Return type:

RealDeviceArray

coordinates()

Return node coordinates as a (n_nodes, 3) array. Returns ——-

coordinatesNDArray[np.float64]

Node coordinates

Return type:

ndarray[Any, dtype[float64]]

coordinates_fields()

Return spatial coordinates for every DOF in the solution vector, shape (n_sol_dofs, dim). Returns ——-

coordsNDArray[np.float64]

Coordinate array

Return type:

ndarray[Any, dtype[float64]]

density()

Return the density field. Returns ——-

densityRealDeviceArray

Density array

Return type:

RealDeviceArray

dimension()

Return the spatial dimension (2 or 3). Returns ——-

dimint

Spatial dimension

Return type:

int

element_size()

Return the representative size of each element, shape (n_elements,). Returns ——-

hRealDeviceArray

Element size array

Return type:

RealDeviceArray

element_tags()

Return element tags. Returns ——-

tagsNDArray[np.int32]

Element tags

Return type:

ndarray[Any, dtype[int32]]

elements()

Return the element connectivity. Returns ——-

elementsNDArray[np.int32]

Element array, shape (n_elements, dim + 1)

Return type:

ndarray[Any, dtype[int32]]

enable_stability(enable_pspg, enable_supg)

Enable or disable PSPG and SUPG stabilisation terms.

Parameters

enable_pspgint

Enable PSPG (pressure-stabilising Petrov-Galerkin): 1=yes

enable_supgint

Enable SUPG (streamline-upwind Petrov-Galerkin): 1=yes

field_indices(ifield)

Return solution-vector indices for one field.

Return type:

ndarray[Any, dtype[int32]]

Parameters

ifieldint

Field number

Returns

idxNDArray[np.int32]

Solution-vector indices

fields_gradient()

Return the continuous gradient of each field at P1 nodes. Shape: (ndof_p1, n_fields, dim). Returns ——-

gradNDArray[np.float64]

Gradient array

Return type:

ndarray[Any, dtype[float64]]

full_implicit_euler(dt, tol=1e-06, reduced_gravity=0, stab_param=0.0, itermax=100)

Advance by one time step using a fully implicit Newton scheme.

Parameters

dtfloat

Time step

tolfloat

Newton convergence tolerance

reduced_gravityint

Use reduced gravity formulation

stab_paramfloat

Additional stabilization parameter

itermaxint

Maximum Newton iterations

g()

Return the gravity vector the problem was constructed with. The coupling closure and the Patankar datum both need it (the datum carries g*dt), and a scheme that builds them has no other way to ask. Returns ——-

gNDArray[np.float64]

Gravity vector, shape (dim,)

Return type:

ndarray[Any, dtype[float64]]

get_bodies_csr_force()

Return fluid/body force at body CSR integration points. Returns ——-

forcesRealDeviceArray

CSR body force array

Return type:

RealDeviceArray

get_bodies_csr_force_penalty()

Return the PENALTY-REACTION part of the per-entry body force. get_bodies_csr_force() returns the total, which also carries the -int_G grad p buoyancy – the traction-like part, which samples the pressure inside the cut band where it is noisiest. This one is the variationally consistent route: the momentum exchange int beta(u-g) plus the trace band, read from the assembled equations. Returns ——-

forceRealDeviceArray

Per-entry penalty reaction, shape (nnz, dim)

Return type:

RealDeviceArray

get_default_export()

Return the default fields for write_mig as {name: (data, element)}. Mirrors Python FluidProblem.get_default_export(). Returns ——-

fieldsDict[str,Tuple[NDArray[np.float64], str]]

Default export fields

Return type:

Dict[str, Tuple[ndarray[Any, dtype[float64]], str]]

get_density_element()

Return the element type string for the density field. Returns ——-

elementstr

Element type string

Return type:

str

get_forces_on_bodies()

Return the fluid force applied on each body, shape (n_bodies, dim). Returns ——-

forcesNDArray[np.float64]

Body force array

Return type:

ndarray[Any, dtype[float64]]

get_mapping(etype)

Return the mapping associated with an element type.

Return type:

ndarray[Any, dtype[int32]]

Parameters

etypestr

Element type, e.g. P1

Returns

mappingNDArray[np.int32]

Element mapping

get_p1_element()

Return the default P1 element name. Returns ——-

elementstr

Element name

Return type:

str

get_p1_mapping()

Return the P1 DOF index for each mesh node. Returns ——-

mappingNDArray[np.int32]

Integer array of length n_nodes

Return type:

ndarray[Any, dtype[int32]]

get_pressure()

Return the pressure field as a (n_p_nodes,) array. Returns ——-

pressureNDArray[np.float64]

Pressure array

Return type:

ndarray[Any, dtype[float64]]

get_pressure_element()

Return the element type string for pressure degrees of freedom. Returns ——-

elementstr

Element type string

Return type:

str

get_velocity_element()

Return the element type string for velocity degrees of freedom. Returns ——-

elementstr

Element type string

Return type:

str

implicit_euler(dt, check_residual_norm=-1, reduced_gravity=False, stab_param=0.0)

Advance the solution one time step with the implicit Euler scheme.

Parameters

dtfloat

Time-step size

check_residual_normfloat

If > 0, throw if residual norm exceeds this value after solve

reduced_gravitybool

Use reduced-gravity formulation

stab_paramfloat

If non-zero, use as pressure-Laplacian stabilisation coefficient instead of PSPG/SUPG

interpolate(solution=None, velocity=None, velocity_x=None, velocity_y=None, velocity_z=None, pressure=None)

Assign solution, velocity, velocity components, or pressure from arrays or callbacks.

Parameters

solutionUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Full solution values, or callback on coordinates_fields()

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

Velocity values, or callback on velocity DOF coordinates

velocity_xUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

X velocity values/callback

velocity_yUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Y velocity values/callback

velocity_zUnion[None, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Z velocity values/callback

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure values/callback

local_boundary_force_by_contribution(tag)

Return per-edge boundary forces split into pressure and viscous contributions. Shape: (n_edges, 2*dim).

Return type:

ndarray[Any, dtype[float64]]

Parameters

tagstr

Boundary name

Returns

forcesNDArray[np.float64]

Per-edge force array, shape (n_edges, 2*dim)

local_size()

Return the number of local degrees of freedom per element. Returns ——-

nint

Local system size

Return type:

int

mesh_boundaries()

Return all mesh boundaries as {name: edge_nodes}. Mirrors Python FluidProblem.mesh_boundaries(). Returns ——-

boundariesDict[str,NDArray[np.int32]]

Map from boundary name to node index array, shape (bsize, dim)

Return type:

Dict[str, ndarray[Any, dtype[int32]]]

mesh_velocity()

Compatibility alias for mesh_velocity_array(). Returns ——-

velRealDeviceArray

Mesh velocity array, shape (n_nodes, 3)

Return type:

RealDeviceArray

mesh_velocity_array()

Return the mesh velocity, in device memory, shape (n_nodes, dim). Returns ——-

velRealDeviceArray

Mesh velocity device array

Return type:

RealDeviceArray

n_elements()

Return the number of mesh elements. Returns ——-

nint

Number of elements

Return type:

int

n_fields()

Return the number of solution fields. Returns ——-

nint

Number of fields

Return type:

int

n_nodes()

Return the number of mesh nodes. Returns ——-

nint

Number of nodes

Return type:

int

node_volume()

Return nodal control volumes. Returns ——-

volumeRealDeviceArray

Node volumes

Return type:

RealDeviceArray

p_jump()

Return the pressure jump field, on P1DG nodes. Returns ——-

p_jumpRealDeviceArray

Pressure jump array

Return type:

RealDeviceArray

particle_volume_intersected()

Return, per body, how much of its volume lies inside the mesh: the sum of its contributions’ volumes. Bodies outside the mesh give zero. Returns ——-

volumeNDArray[np.float64]

Intersected volume per body, shape (n_bodies,)

Return type:

ndarray[Any, dtype[float64]]

porosity()

Compatibility alias for porosity_array(). Returns ——-

porRealDeviceArray

Porosity array, shape (n_nodes, 1)

Return type:

RealDeviceArray

porosity_array()

Return the porosity, in device memory, shape (n_nodes,). Returns ——-

porRealDeviceArray

Porosity device array

Return type:

RealDeviceArray

pressure()

Compatibility alias for get_pressure(). Returns ——-

pressureNDArray[np.float64]

Pressure array

Return type:

ndarray[Any, dtype[float64]]

read_mig(odir, t=-1.0, iteration=-2147483648)

Read a fluid state written by write_mig().

Parameters

odirstr

Output directory

tfloat

Time to read. If omitted, iteration must be provided.

iterationint

Iteration index, or -1 for last.

reset_timings()

Zero the timing counters, to skip a warm-up phase.

set_concentration_cg(concentration)

Set the concentration field from continuous nodal values.

Parameters

concentrationNDArray[np.float64]

Concentration array, shape (n_nodes,)

set_coordinates(x)

Move the mesh nodes.

Parameters

xNDArray[np.float64]

Node coordinates, shape (n_nodes, 3)

set_coupling_closure(csr_m, csr_mgsp, csr_div_m, csr_beta, csr_h, csr_velocity, csr_contact, csr_eta=None, csr_eta_p=None, csr_gamma_t=None, csr_stab_w=None)

TIER B – the model on the current geometry: the unfitted map and its divergence, the penalty, the trace band, the frozen velocity datum and the contact force. Requires a geometry of the current generation.

Parameters

csr_mNDArray[np.float64]

Unfitted map m, shape (nnz, dim)

csr_mgspNDArray[np.float64]

m . grad(u_p), shape (nnz, dim)

csr_div_mNDArray[np.float64]

div m, shape (nnz,)

csr_betaNDArray[np.float64]

Penalty coefficient, shape (nnz,)

csr_hNDArray[np.float64]

Trace length, shape (nnz,)

csr_velocityNDArray[np.float64]

Frozen body-velocity datum, shape (nnz, dim)

csr_contactNDArray[np.float64]

Contact force per entry, shape (nnz, dim)

csr_etaNDArray[np.float64]

DIMENSIONLESS interior-fill coefficient per entry; absent falls back to the global set_eta

csr_eta_pNDArray[np.float64]

DIMENSIONLESS pressure-diffusivity coefficient per entry; absent falls back to the global set_eta_p

csr_gamma_tNDArray[np.float64]

DIMENSIONLESS trace coefficient per entry; absent means the global set_trace_param value. Switch the band off with 0.0 – never by sending csr_h to infinity, which is what a 1e150 sentinel used to mean and put its square one decade from overflow.

csr_stab_wNDArray[np.float64]

Fraction of the coupling coefficient that is PHYSICAL DRAG rather than numerical penalty, per entry, in [0, 1]. It is the weight with which tricks 2 and 3 enter the SUPG/PSPG residual: 1 puts the whole coupling there, as the legacy mixture kernel did with its drag; 0 keeps it out, which a Babuska penalty requires because betap*pspg = gamma_B/4 is mesh-independent and would never refine away. The producer sets it to the drag fraction of its own blend, so a sub-cell grain gets 1 and a resolved one gets 0. Absent falls back to the global drag_in_stab flag.

set_coupling_datum(csr_velocity, csr_dup_ds, csr_dup_dds)

TIER C – the affine datum a scheme supplies so the body velocity responds to the fluid INSIDE the solve: up_eff = velocity + dup_ds . s + dup_dds . ds. Optional; with no datum set for the current generation the closure’s frozen velocity stands, which is what iterate_iqn wants. The velocity given here REPLACES the closure’s (same slot, better value).

Parameters

csr_velocityNDArray[np.float64]

Constant part up0, shape (nnz, dim)

csr_dup_dsNDArray[np.float64]

d u_p / d s_F, shape (nnz, dim, dim + 1)

csr_dup_ddsNDArray[np.float64]

d u_p / d (ds_F/dx), shape (nnz, dim, dim + 1, dim)

set_coupling_geometry(density, volume, csr_ptr, csr_element_id, csr_xi, csr_surface, csr_weight, csr_d, csr_cut=None)

TIER A – the overlap geometry, and the per-body constants that go with it. Pure of state: no velocity, no fluid solution, no dt. Setting it BUMPS THE COUPLING GENERATION and invalidates any closure and datum, so a stale closure sized to a previous geometry, or last step’s datum left behind after the bodies moved, are refused rather than silently used. Setting a geometry and NOT completing it with a closure is refused at solve time too – it used to be a silent no-op that ran the operator against the previous commit. A held bed can call this once for the whole run; IQN once per step rather than once per iteration.

Parameters

densityNDArray[np.float64]

Per-body density, shape (n_bodies,)

volumeNDArray[np.float64]

Per-body volume, shape (n_bodies,)

csr_ptrNDArray[np.int32]

Body -> entry offsets, shape (n_bodies + 1,)

csr_element_idNDArray[np.int32]

Element of each entry, shape (nnz,)

csr_xiNDArray[np.float64]

Reference coordinates, shape (nnz, dim)

csr_surfaceNDArray[np.float64]

Interface measure per entry, shape (nnz,)

csr_weightNDArray[np.float64]

Integration weight per entry, shape (nnz,)

csr_dNDArray[np.float64]

Signed distance per entry, shape (nnz,)

csr_cutNDArray[np.int32]

TOPOLOGICAL per-entry cut flag (1 = the element is only PARTIALLY covered by this body). For a convex body this is exactly nin < dim+1, with no threshold anywhere. Absent falls back to the kernel’s old volume-fraction inference.

set_eta(eta)

Set the interior grain-Laplacian fill coefficient of the unfitted kernel.

Parameters

etafloat

Dimensionless fill coefficient, applied as eta * mu

set_eta_p(eta_p)

Set the covered-pressure Laplacian coefficient of the unfitted kernel. Applied as eta_p * taup/rho on the covered (P,P) block. 1e-3 is the production default; the legacy-retrieval configuration passes 0.

Parameters

eta_pfloat

Dimensionless covered-pressure coefficient

set_g(g)

Set the gravity vector (body force per unit mass). The constructor value is otherwise immutable, and g() hands back a COPY, so a scheme that has to STEER the body force – a periodic cell driven to a target flow rate, where the driver is the unknown and the velocity is imposed – has no way to write it back. This is that way.

Parameters

gNDArray[np.float64]

Gravity vector, shape (dim,)

set_mean_pressure(p)

Constrain the mean pressure to a fixed value.

Parameters

pfloat

Target mean pressure value

set_mesh(nodes, elements, element_tags, boundaries, boundary_tags, boundary_names, periodic=None)

Set the mesh from arrays.

Parameters

nodesNDArray[np.float64]

Node coordinates, shape (n_nodes, 3)

elementsNDArray[np.int32]

Element connectivity, shape (n_elements, dim+1)

element_tagsNDArray[np.int32]

Element tags, shape (n_elements,)

boundariesNDArray[np.int32]

Boundary edge nodes, shape (n_bnd, dim)

boundary_tagsNDArray[np.int32]

Boundary tags, shape (n_bnd,)

boundary_namesList[str]

Physical group names

periodicNDArray[np.int32]

Parent DOF index for each node (n_nodes,), or None

set_open_boundary(tag, velocity=None, pressure=None, viscous_flag=True, porous=True)

Set an open (type=1, viscous_flag=1) boundary condition.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

viscous_flagbool

Compute viscous term at the boundary

porousbool

Weight this boundary’s data by the CURRENT porosity. True imposes a SUPERFICIAL quantity, as the mixture kernel always did (f0(P) = c*unext): a bed sitting on an inlet then admits only c times the nominal velocity. False imposes the datum as given. Measured on 2d_bubbling_fluidized_bed, where the bed sits on the inlet at c = 0.5982: imposing the datum unweighted made the whole solution 1.683x too large – velocity, pressure and grain force alike, and in grain-free regions too, because an inflow-flux error is global. False is what a traction outlet wants: the gradient-form momentum reduces to the imposed traction only at c = 1 (3d_darcy_permeability, c = 0.478, showed p_out = 0.0354 against an imposed 0 and a 10.1% grain-force overshoot at c < 1).

set_particles_body(xp, rp, density, velocity=None, omega=None, contact_forces=None, discretisation='overlap', use_voidage=True)

Set body coupling data from particle centers/radii and particle fields.

Parameters

xpNDArray[np.float64]

Particle positions, shape (n_particles, dim)

rpNDArray[np.float64]

Particle radii, shape (n_particles,)

densityNDArray[np.float64]

Particle densities, shape (n_particles,)

velocityNDArray[np.float64]

Particle velocities, shape (n_particles, dim)

omegaNDArray[np.float64]

Particle angular velocities, shape (n_particles, dim)

contact_forcesNDArray[np.float64]

Particle contact forces, shape (n_particles, dim)

discretisationstr

overlap in 2D or centroid otherwise

use_voidagebool

Include voidage in drag coefficient

set_particles_body_csr(xp, rp, density, csr_ptr, csr_element_id, csr_xi, csr_surface, csr_weight, velocity=None, omega=None, contact_forces=None, use_voidage=True)

Set body coupling data from particle fields and a CALLER-SUPPLIED quadrature CSR – the mixture-path twin of set_unfitted_bodies. The overlap geometry is produced outside the core (overlap, through migflow.overlap_quadrature); the solver only copies. An entry is a quadrature point: csr_weight integrates, csr_surface is the pair-total overlap duplicated on each entry (cut-band detection only).

Parameters

xpNDArray[np.float64]

Particle positions, shape (n_particles, dim)

rpNDArray[np.float64]

Particle radii, shape (n_particles,)

densityNDArray[np.float64]

Particle densities, shape (n_particles,)

csr_ptrNDArray[np.int32]

CSR row pointer, shape (n_particles + 1,)

csr_element_idNDArray[np.int32]

Element of each entry, shape (nnz,)

csr_xiNDArray[np.float64]

Reference coordinates of each point, shape (nnz, dim)

csr_surfaceNDArray[np.float64]

Pair-total overlap area/volume, shape (nnz,)

csr_weightNDArray[np.float64]

Integration weight of each point, shape (nnz,)

velocityNDArray[np.float64]

Particle velocities, shape (n_particles, dim)

omegaNDArray[np.float64]

Particle angular velocities, shape (n_particles, dim)

contact_forcesNDArray[np.float64]

Particle contact forces, shape (n_particles, dim)

use_voidagebool

Include voidage in drag coefficient

set_polygon_bodies(vertices, density, volume, gamma, velocity=None, contact_forces=None)

Set body coupling data from convex polygons rather than discs. The drag coefficient is taken from the caller instead of being derived from a radius, which a polygon does not have.

Parameters

verticesNDArray[np.float64]

Polygon vertices, shape (n_bodies, n_vertices, dim)

densityNDArray[np.float64]

Body densities, shape (n_bodies,)

volumeNDArray[np.float64]

Body volumes, shape (n_bodies,)

gammaNDArray[np.float64]

Drag coefficient of each body, shape (n_bodies,)

velocityNDArray[np.float64]

Body velocities, shape (n_bodies, dim)

contact_forcesNDArray[np.float64]

Body contact forces, shape (n_bodies, dim)

set_pressure(pres)

Set the pressure field from a (n_p_nodes,) array.

Parameters

presNDArray[np.float64]

Pressure array

set_stabilisation_complement(pspg=True, supg=True)

Choose which halves of f_volume’s stabilisation image the unfitted body kernel subtracts over the element-grain overlap. PSPG is needed for equal-order P1/P1 and its complement is uncontroversial. The SUPG complement is not: f_volume’s SUPG image carries the convective term, and over the covered region the advecting field is the extended velocity INSIDE the grain, which is fictitious. Subtracting a convective stabilisation built from it is a modelling choice, invisible at Stokes – which is exactly where this formulation agrees with cutfem.

Parameters

pspgbool

Complement the PSPG image

supgbool

Complement the SUPG image

set_strong_boundary(tag, velocity=None, velocity_x=None, velocity_y=None, velocity_z=None)

Set a strong (Dirichlet) boundary condition. velocity imposes every component at once; velocity_x/y/z impose one, leaving the others free, which is what a slip wall needs.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

velocity_xUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

X velocity value/callback, or None

velocity_yUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Y velocity value/callback, or None

velocity_zUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Z velocity value/callback, or None

set_symmetry_boundary(tag, pressure=None)

Set a symmetry boundary: only the normal flux is imposed. This is a weak boundary with no velocity and viscous_flag left unspecified, which with no velocity imposed resolves to a traction-free surface.

Parameters

tagstr

Boundary physical-group name

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

set_trace_param(gamma_t)

Set the trace-stabilisation coefficient of the unfitted kernel.

Parameters

gamma_tfloat

Dimensionless trace coefficient

set_two_fluid_properties(rho, mu, sigma=0.0)

Derive the density and viscosity fields from the concentration, and optionally the capillary force. The density and viscosity elements must both be P1DG, like the concentration. This helper owns bulk_force: it overwrites it on every call when sigma is non-zero. Add any other bulk force contribution afterwards.

Parameters

rhoList[float]

Density of each of the two fluids

muList[float]

Dynamic viscosity of each of the two fluids

sigmafloat

Surface tension coefficient; 0 leaves bulk_force alone

set_unfitted(unfitted)

Select the unfitted (boundary-free) body model. With it enabled the no-slip datum on a grain boundary is imposed purely through volume integrals over the element-grain overlap – no facet loop, no interface integral – and the nodal porosity no longer weights the volume operators: the geometry enters as the sharp complement int_K (.) - int_{K n G} (.).

Parameters

unfittedbool

Whether to use the unfitted kernel

set_unfitted_bodies(density, volume, csr_ptr, csr_element_id, csr_xi, csr_surface, csr_weight, csr_velocity, csr_contact, csr_m, csr_mgsp, csr_div_m, csr_d, csr_h, csr_beta, csr_dup_ds=None, csr_dup_dds=None, csr_eta=None, csr_eta_p=None, csr_gamma_t=None, csr_cut=None, csr_stab_w=None)

Hand the unfitted overlap table to the solver as a flat CSR. An entry is a QUADRATURE POINT, not a (body, element) pair: csr_element_id repeats within a body, csr_xi is the point in reference coordinates, and csr_weight is its integration weight. csr_surface stays the pair total, duplicated on every entry, and is used only for the geometric cut test and as a diagnostic. The geometry is produced entirely caller-side; the solver only copies. Arguments follow the field order of the internal FluidProblemBodies.

Parameters

densityNDArray[np.float64]

Per-body density, shape (n_bodies,)

volumeNDArray[np.float64]

Per-body volume, shape (n_bodies,)

csr_ptrNDArray[np.int32]

CSR row pointer, shape (n_bodies + 1,)

csr_element_idNDArray[np.int32]

Element of each entry, shape (nnz,)

csr_xiNDArray[np.float64]

Reference coordinates of each point, shape (nnz, dim)

csr_surfaceNDArray[np.float64]

Pair-total overlap area/volume, shape (nnz,)

csr_weightNDArray[np.float64]

Integration weight of each point, shape (nnz,)

csr_velocityNDArray[np.float64]

Imposed velocity datum u_p, shape (nnz, dim). When the sensitivity arrays below are supplied this is the CONSTANT PART of an AFFINE datum: the kernel imposes up_eff = csr_velocity + dup_ds . s + dup_dds . ds, evaluated at each quadrature point’s own state, in every residual and force line. The caller must therefore fold its linearisation reference into csr_velocity analytically – csr_velocity = up_pred - dup_ds . s_ref - dup_dds . ds_ref – so that up_eff reproduces up_pred at the reference state. Left null the datum is frozen and csr_velocity is u_p as-is.

csr_contactNDArray[np.float64]

Contact force at each point, shape (nnz, dim)

csr_mNDArray[np.float64]

Direction field (x-c)/r, shape (nnz, dim)

csr_mgspNDArray[np.float64]

m . grad(u_p) per velocity component, shape (nnz, dim)

csr_div_mNDArray[np.float64]

Divergence of m, shape (nnz,)

csr_dNDArray[np.float64]

Signed distance, negative inside the grain, shape (nnz,)

csr_hNDArray[np.float64]

Penalty length, constant per grain, shape (nnz,)

csr_betaNDArray[np.float64]

Babuska penalty density gamma/h^2, shape (nnz,)

csr_dup_dsNDArray[np.float64]

d u_p,i / d s_F, shape (nnz, dim, dim + 1)

csr_dup_ddsNDArray[np.float64]

d u_p,i / d (ds_F/dx_a), shape (nnz, dim, dim + 1, dim)

csr_etaNDArray[np.float64]

DIMENSIONLESS interior-fill coefficient per entry; absent falls back to the global set_eta

csr_eta_pNDArray[np.float64]

DIMENSIONLESS pressure-diffusivity coefficient per entry; absent falls back to the global set_eta_p

csr_gamma_tNDArray[np.float64]

DIMENSIONLESS trace coefficient per entry; absent means the global set_trace_param value. betat = mu*gamma_t/h^2, so the band goes off with gamma_t = 0.0 – exactly, and with h left a real length – rather than by sending h to infinity.

csr_cutNDArray[np.int32]

TOPOLOGICAL per-entry cut flag (1 = the element is only PARTIALLY covered by this body). Absent falls back to the kernel’s old volume-fraction inference, surface < 0.99*|K|. The two sensitivity arrays are how a SEMI-IMPLICIT coupling is expressed without a kernel branch. Left null (the default) the imposed velocity is a frozen Dirichlet datum, every Jacobian contribution below is multiplied by zero, and the assembly is bit-for-bit what it was before they existed – that is the property the explicit and the iterating schemes rely on, where dv = 0 is EXACT for the statement being linearised. Filled, they let the caller declare that u_p depends on the fluid state at the same solve, which is what the legacy mixture kernel did internally with its up_pred: the old du_pred_du is csr_dup_ds[i][U+i] and the old du_pred_dp is csr_dup_dds[i][P][i]. The field index F runs over the NF solution fields, U + 0 .. U + dim - 1 then P. SUPPLYING A SENSITIVITY THAT DOES NOT MATCH THE VELOCITY BESIDE IT IS SILENT AND CATASTROPHIC – implicit_euler is a single Newton step, so a mismatch corrupts the answer outright rather than merely slowing convergence. Measured on the no-slip MMS: a mismatched pair turned a 63% force error into 131512%, and attaching a sensitivity implying a mobility to a residual holding the velocity fixed gave 69479% against 245%. Derive all three from one closure (see body_coupling.get_patankar_pred) and check them with body_coupling.check_sensitivity before trusting a number.

csr_stab_wNDArray[np.float64]

Fraction of the coupling coefficient that is PHYSICAL DRAG rather than numerical penalty, per entry, in [0, 1]. The weight with which tricks 2 and 3 enter the SUPG/PSPG residual; see set_coupling_closure. Absent falls back to the global drag_in_stab.

set_velocity(vel)

Set the velocity field from a (n_vel_nodes, dim) array.

Parameters

velNDArray[np.float64]

Velocity array, shape (n_vel_nodes, dim)

set_wall_boundary(tag, velocity=None, pressure=None, viscous_flag=-1, porous=True)

Set a wall-type (type=0) boundary condition.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

viscous_flagint

Compute viscous term at the boundary: 1 to force it on, 0 to force it off, -1 (default) to compute it only when a velocity is imposed

porousbool

Weight this boundary’s data by the CURRENT porosity. True imposes a SUPERFICIAL quantity, as the mixture kernel always did (f0(P) = c*unext): a bed sitting on an inlet then admits only c times the nominal velocity. False imposes the datum as given. Measured on 2d_bubbling_fluidized_bed, where the bed sits on the inlet at c = 0.5982: imposing the datum unweighted made the whole solution 1.683x too large – velocity, pressure and grain force alike, and in grain-free regions too, because an inflow-flux error is global. False is what a traction outlet wants: the gradient-form momentum reduces to the imposed traction only at c = 1 (3d_darcy_permeability, c = 0.478, showed p_out = 0.0354 against an imposed 0 and a 10.1% grain-force overshoot at c < 1).

set_weak_boundary(tag, velocity=None, pressure=None, concentration=None, viscous_flag=-1, porous=True)

Set a weak (type=0) boundary condition.

Parameters

tagstr

Boundary physical-group name

velocityUnion[None, List[Union[float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]], NDArray[np.float64]]

List of dim components (Real or BoundaryComponentFn), or None

pressureUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Pressure value (Real or BoundaryComponentFn), or None

concentrationUnion[None, float, NDArray[np.float64], Callable[[NDArray[np.float64]], NDArray[np.float64]]]

Concentration/alpha value for two-fluid problems, or None

viscous_flagint

Compute viscous term at the boundary: 1 to force it on, 0 to force it off, -1 (default) to enable it only when a velocity or a pressure is imposed

porousbool

Weight this boundary’s data by the CURRENT porosity; see set_open_boundary.

set_weak_boundary_nodal_values(tag, values)

Replace the nodal values of a weak boundary already set by set_weak_boundary(). Unlike set_weak_boundary() this does not add a boundary, so it can be called at every time step.

Parameters

tagstr

Boundary physical-group name

valuesNDArray[np.float64]

One value per field and per boundary node, shape (n_edges, n_closure_dofs, n_values)

solution()

Compatibility alias for solution_array(). Returns ——-

solRealDeviceArray

Solution array of length n_dof

Return type:

RealDeviceArray

solution_array()

Return the full solution, in device memory, as an array of length n_dof. Returns ——-

solRealDeviceArray

Solution device array of length n_dof

Return type:

RealDeviceArray

solution_at_coordinates(x)

Interpolate the solution at arbitrary coordinates.

Return type:

ndarray[Any, dtype[float64]]

Parameters

xNDArray[np.float64]

Coordinates, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

solution_at_reference(eid, xi)

Interpolate the solution at points given by element id and REFERENCE coordinates. No mesh search: use this whenever the caller already knows which element each point is in – an overlap CSR does, and going in by physical position instead costs a locate_point_in_mesh per point and is ambiguous on a face. Evaluated with the field’s OWN shape functions, so it is correct at any element order (a hand-rolled P1 formula in the caller is not).

Return type:

ndarray[Any, dtype[float64]]

Parameters

eidNDArray[np.int32]

Element index per point, shape (n,)

xiNDArray[np.float64]

Reference coordinates in that element, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

solution_grad_at_coordinates(x)

Solution and gradient at arbitrary physical coordinates: locates each point, then evaluates as solution_grad_at_reference does.

Return type:

Tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]]]

Parameters

xNDArray[np.float64]

Coordinates, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

gradientNDArray[np.float64]

Solution gradient, shape (n, n_fields, dim)

solution_grad_at_reference(eid, xi)

Solution AND its gradient at element id + reference coordinates, in one pass over the shape functions.

Return type:

Tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[float64]]]

Parameters

eidNDArray[np.int32]

Element index per point, shape (n,)

xiNDArray[np.float64]

Reference coordinates in that element, shape (n, dim)

Returns

valuesNDArray[np.float64]

Interpolated solution, shape (n, n_fields)

gradientNDArray[np.float64]

Solution gradient, shape (n, n_fields, dim)

timings()

Return the wall-clock breakdown of implicit_euler, in seconds. Slots, in order: setup, assemble, scatter, factorise, solve, bodies, n_steps. Each mark synchronises the device first, because the assembly kernels are asynchronous and an unsynced timestamp would charge their cost to whichever phase blocks next. That serialisation makes the TOTAL a slight over-estimate; the split is what it is for. Returns ——-

timingsNDArray[np.float64]

Seven values: the six phases and the step count

Return type:

ndarray[Any, dtype[float64]]

u_solid()

Return the solid velocity field at P1 nodes, shape (n_p1_nodes, dim). Returns ——-

usNDArray[np.float64]

Solid velocity array

Return type:

ndarray[Any, dtype[float64]]

update_node_volume()

Rebuild the cached mesh metrics, including the nodal control volumes, after the mesh has been deformed in place through coordinates().

velocity()

Return the velocity field as a (n_vel_nodes, dim) array. Returns ——-

velocityNDArray[np.float64]

Velocity array

Return type:

ndarray[Any, dtype[float64]]

velocity_dof_coordinates()

Return spatial coordinates of each velocity DOF, shape (n_vel_nodes, dim). Returns ——-

coordsNDArray[np.float64]

Coordinate array

Return type:

ndarray[Any, dtype[float64]]

velocity_index()

Return solution-vector indices for velocity DOFs, shape (n_vel_dofs, dim). velocity_index()[i, d] is the solution index of component d at velocity DOF i. Returns ——-

idxNDArray[np.int32]

Index array

Return type:

ndarray[Any, dtype[int32]]

viscosity()

Return the dynamic-viscosity field. Returns ——-

viscosityRealDeviceArray

Viscosity array

Return type:

RealDeviceArray

write_mig(output_dir, t, fields={})

Write output files for post-visualisation. Mirrors Python FluidProblem.write_mig().

Parameters

output_dirstr

Output directory

tfloat

Computational time

fieldsDict[str,Tuple[NDArray[np.float64], str]]

Fields to write as {name: (data, element)}. Uses get_default_export() if empty.

class migflow.fluid.ILinearSystem

Bases: object

Plug-in linear solver interface. Use ILinearSystem.create_from_callbacks() to create a custom solver and pass it to FluidProblem.create() as the solver= argument.

Pointer

alias of LP__Structure

static create_amgx(config='')

Create an AmgX iterative solver on the GPU (linked in, CUDA builds only). Measured on the fluid system it converges but loses to PETSc-LU at every size, its iteration count growing linearly with the unknowns; use create_cudss() instead. Cannot solve a problem with Lagrange constraints.

Return type:

ILinearSystem

Parameters

configstr

AmgX JSON configuration, or the path to a file holding one. Empty means the default configuration.

Returns

solverILinearSystem

Solver object for FluidProblem.create()

static create_cudss()

Create a cuDSS direct solver running on the GPU. The only solver that keeps the matrix on the device: nothing is copied back to the host between the assembly and the solution. It matches PARDISO to 1e-15 and is about 7x faster than PETSc-LU; its solve overtakes PARDISO’s around 90000 unknowns, while its factorisation never does. Worth it on a CUDA build from mid-sized problems up, and for many solves per factorisation. Requires a CUDA build (-DAVA_TARGET=CUDA); raises otherwise. Returns ——-

solverILinearSystem

Solver object for FluidProblem.create()

Return type:

ILinearSystem

static create_default(options='')

Try PARDISO first, then PETSc; raise if neither library is available.

Return type:

ILinearSystem

Parameters

optionsstr

PETSc options string forwarded to create_petsc() if used

Returns

solverILinearSystem

Solver object for FluidProblem.create()

static create_from_callbacks(on_create, on_solve)

Create a callback-based linear solver.

Return type:

ILinearSystem

Parameters

on_createCallable[[NDArray[np.uint32], NDArray[np.uint32], int], None]

Called once per sparsity pattern with (row_offsets, columns, n_rows). Pass None to skip.

on_solveCallable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]]

Called for each solve: (values, rhs) -> solution (length >= n_rows)

Returns

solverILinearSystem

Solver object for FluidProblem.create()

static create_pardiso()

Create an MKL PARDISO direct solver (loaded at runtime from libmkl_rt). Returns ——-

solverILinearSystem

Solver object for FluidProblem.create()

Return type:

ILinearSystem

static create_petsc(options='')

Create a PETSc iterative solver (loaded at runtime from libpetsc).

Return type:

ILinearSystem

Parameters

optionsstr

PETSc options string passed to PetscOptionsInsertString

Returns

solverILinearSystem

Solver object for FluidProblem.create()

class migflow.fluid.RealDeviceArray

Bases: object

Scalar array on the device memory (GPU)

Pointer

alias of LP__Structure

property cap: int

the total capacity

static create_empty(shape)

create a newly allocated (and uninitialized) ava device array of scalars

Return type:

RealDeviceArray

Parameters

shapeList[int]

The shape of the array

Returns

new_device_arrayRealDeviceArray

The newly created device array

static create_from_host(host_array)

create a newly allocated ava device array of scalars and set it from a host array

Return type:

RealDeviceArray

Parameters

host_arrayNDArray[np.float64]

The host array to copy the data from.

Returns

new_device_arrayRealDeviceArray

The newly created device array

static create_zeros(shape)

create a zero-initialized ava device array of scalars

Return type:

RealDeviceArray

Parameters

shapeList[int]

The shape of the array

Returns

new_device_arrayRealDeviceArray

The newly created device array

duplicate()

create a copy of a DeviceArray Returns ——-

new_device_arrayRealDeviceArray

The newly created device array

Return type:

RealDeviceArray

fill(value, stream=None)

Fill a DeviceArray with a value, akin to std::fill

Parameters

valuefloat

The scalar value to set

streamStream|None

The stream to use

get()

get a ava::HostArray on the CPU from a DeviceArray on the GPU Returns ——-

new_host_arrayNDArray[np.float64]

The newly created host array. The data is copied from the device to the host

Return type:

ndarray[Any, dtype[float64]]

has_nan(stream=None)

Return type:

bool

Parameters

streamStream|None

The stream to use

Returns

hasnanbool

True if the array contains any NaN

memset(value=0)

set all bytes of a DeviceArray to value

Parameters

valueint

The byte with which to fill the DeviceArray

property name: str

the name of the array

resize(shape)

Resize an existing ava device array, allowing to reuse previous memory. If the total size of new_shape is greater or equal to the current, the underlying data is unchanged. If it is smaller, the data is truncated.

Parameters

shapeList[int]

The new shape of the array

resize_no_copy(shape)

Resize an existing ava device array, allowing to reuse previous memory. The underlying data is not maintained.

Parameters

shapeList[int]

The new shape of the array

set(host_array)

set a DeviceArray on the GPU from a ava::HostArray on the CPU alias for set_from_host (deprecated, for backward compatibility)

Parameters

host_arrayNDArray[np.float64]

The host array. The data is copied from the host to the device

set_from_device(device_array)

set a DeviceArray on the GPU from another DeviceArray on the GPU

Parameters

device_arrayRealDeviceArray

The device array. The data is copied from the device to the device

set_from_host(host_array)

set a DeviceArray on the GPU from a HostArray on the CPU

Parameters

host_arrayNDArray[np.float64]

The host array. The data is copied from the host to the device

set_name(name)

set the name of a DeviceArray

Parameters

namestr

The name of the array

property shape: List[int]

the shape of the array

property size: int

the total number of elements

transpose(stream=None)

transpose an array in a new array (only for 3-dimensional arrays)

Return type:

RealDeviceArray

Parameters

streamStream|None

The stream to use

Returns

destRealDeviceArray

Transposed array

class migflow.fluid.Stream

Bases: object

RAII wrapper around a backend stream handle (cudaStream_t / hipStream_t / void* on CPU). Created through the create() factory and handed around as Stream::Ptr. The stream is destroyed when the last Ptr is released. On the CPU backend stream_create / stream_destroy / stream_synchronize are no-ops.

Create a newly allocated stream.

Pointer

alias of LP__Structure

synchronize()

Block the host until all work previously enqueued on this stream completes.

migflow.fluid.real_epsilon()

Return std::numeric_limits<Real>::epsilon() for this build. Returns ——-

epsilonfloat

Machine epsilon for Real

Return type:

float

migflow.fluid.real_size()

Return the size in bytes of the scalar Real type used by this build. Returns ——-

sizeint

Scalar size in bytes

Return type:

int

migflow.fluid.real_type_name()

Return the scalar Real type name used by this build. Returns ——-

namestr

Scalar type name

Return type:

str