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=True, 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='', pressure_enrichment=False, surface_tension_model='csf', sigma_ij=None, sigma=0.0, surface_tension_theta=0.0)

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 drag in stabilisation term

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.

pressure_enrichmentbool

Pressure continuous in each material but discontinuous between them: every node gets one pressure dof per material of its elements (set_materials()). 2D and P1P1 only.

surface_tension_modelstr

csf (the capillary force of set_two_fluid_properties()) or membrane (a tension on the edges between materials, requires pressure_enrichment)

sigma_ijNDArray[np.float64]

Membrane only: surface tensions, a symmetric (n_materials+1, n_materials+1) matrix with a zero diagonal, the materials then a solid (last index, see set_wall_tension()). Fixes n_materials.

sigmafloat

Membrane only, when sigma_ij is not given: the tension between the two materials of a two-material problem (sigma_ij = [[0, sigma, 0], [sigma, 0, 0], [0, 0, 0]])

surface_tension_thetafloat

Membrane only, see set_surface_tension_theta()

Pointer

alias of LP__Structure

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

Adapt the mesh, projecting the current solution onto the new mesh. With pressure_enrichment only the velocity is projected: the materials of the new mesh are given by materials (or by set_materials() before the next solve), and the pressure restarts from 0.

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

materialsNDArray[np.float64]

pressure_enrichment only: the material of each new element, see set_materials()

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(solution, solution_old, dt)

Assemble the linear system of a Newton step of implicit_euler() at a given state, without changing the problem: the residual R and its jacobian J, as the solver gets them (strong boundaries and constraints included, one row per constraint after the n_dof rows of the solution). The stabilisation parameters are those of the state. The state of the problem (solution, old solution) is restored afterwards.

Return type:

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

Parameters

solutionNDArray[np.float64]

State, shape (n_dof,)

solution_oldNDArray[np.float64]

State at the beginning of the time step, shape (n_dof,)

dtfloat

Time step

Returns

row_ptrNDArray[np.int32]

CSR row pointers of J, shape (n_rows + 1,)

columnsNDArray[np.int32]

CSR column indices of J, shape (nnz,)

valuesNDArray[np.float64]

CSR values of J, shape (nnz,)

residualNDArray[np.float64]

R, shape (n_rows,)

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

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

dof_map_version()

Return the number of times the dof maps have been rebuilt (a new mesh, new materials): it changes exactly when the numbering of the solution, and so the linear system, does. Returns ——-

versionint

Dof map version

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: the volume weighted mean of the element gradients around each node. With pressure_enrichment the pressure gradient at a node between materials mixes those of both. 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

get_bodies_csr_force()

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

forcesRealDeviceArray

CSR body force array

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: the value of every pressure dof, in the order of pressure_index(). With pressure_enrichment these are not nodal values: a node between materials has one per material, see pressure_dg(). 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

global_map()

Return the solution-vector index of every local dof of every element, fields one after the other (the velocity components, then the pressure). Returns ——-

mapNDArray[np.int32]

Global dof map, shape (n_elements, local_size)

Return type:

ndarray[Any, dtype[int32]]

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)

materials()

Return a copy of the material of every element (see set_materials()). Returns ——-

labelsNDArray[np.int32]

Material of each element, shape (n_elements,)

Return type:

ndarray[Any, dtype[int32]]

membrane_energy()

Return the membrane energy, the sum of sigma L over the edges between materials and the edges of the boundaries of set_wall_tension(), on the current coordinates. Returns ——-

energyfloat

Membrane energy

Return type:

float

membrane_energy_gradient()

Return the derivative of membrane_energy() with respect to the position of every node: the explicit membrane term of the velocity equations (the residual, not a force: its opposite pulls the nodes). Returns ——-

gradientNDArray[np.float64]

Energy gradient, shape (n_vel_nodes, dim)

Return type:

ndarray[Any, dtype[float64]]

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_materials()

Return the number of materials the labels of set_materials() range over (from sigma_ij, 2 otherwise). Returns ——-

nint

Number of materials

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 prescribed pressure jump field, on P1DG nodes. Not available with pressure_enrichment, where the jump is an unknown. 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]]

pressure_dg()

Return the pressure at the nodes of every element, which keeps the jump between materials of pressure_enrichment. Returns ——-

pressureNDArray[np.float64]

Element pressure, shape (n_elements, local pressure size)

Return type:

ndarray[Any, dtype[float64]]

pressure_enrichment()

Return whether the pressure is enriched (see create()). Returns ——-

enrichedbool

True with pressure_enrichment

Return type:

bool

pressure_index()

Return the solution-vector index of every pressure dof, in increasing order: the node order for the first n_nodes of a P1 pressure, then, with pressure_enrichment, the copies at the nodes between materials. Returns ——-

idxNDArray[np.int32]

Pressure dof indices, shape (n_pressure_dofs,)

Return type:

ndarray[Any, dtype[int32]]

pressure_node_volume()

Return the integration weight of every pressure dof, in the order of pressure_index(): the integral of its shape function over the elements that use it (with pressure_enrichment, those of its material). With a continuous P1 pressure it is node_volume(). These are the weights of the mean pressure constraint. Returns ——-

volumeNDArray[np.float64]

Weight of each pressure dof, shape (n_pressure_dofs,)

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.

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_inflow(tag, velocity, concentration=None, superficial_velocity=False)

An inflow: the velocity is imposed (Nitsche, with the viscous traction) and its flux u_G . n enters the continuity equation. set_weak_boundary(tag, ‘data’, velocity, nitsche=True, viscous_traction=True, concentration=concentration, superficial_velocity=superficial_velocity).

Parameters

tagstr

Boundary physical-group name

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

Velocity data, list of dim components (Real or BoundaryComponentFn) or nodal values

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

Concentration/alpha of the fluid entering, or None

superficial_velocitybool

If True, velocity is the superficial velocity (fluid fraction times the fluid velocity, the flux per unit area)

set_materials(labels)

Set the material of every element, the only definition of the phases with pressure_enrichment: the pressure is continuous inside a material and discontinuous across the mesh edges between two of them. New labels rebuild the pressure dofs: the velocity is kept, the pressure is reset to 0 (a time step gives it back) and the linear system is rebuilt; unchanged labels change nothing. They are reset to 0 by set_mesh() and must be given again after adapt_mesh().

Parameters

labelsNDArray[np.float64]

Material of each element, integer values in [0, n_materials), shape (n_elements,); any other value (fractional, NaN, inf, out of range) is refused before anything changes

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_outflow(tag, pressure=0.0, viscous_traction=False)

An outflow, or a free surface: the flux u . n is free and the pressure is imposed. set_weak_boundary(tag, ‘free’, pressure=pressure, viscous_traction=viscous_traction).

Parameters

tagstr

Boundary physical-group name

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

Pressure p_G (Real, nodal values or BoundaryComponentFn)

viscous_tractionbool

Keep the viscous traction -mu (grad u + grad u^T) n (the natural outflow of the viscous term); False: the normal stress is the pressure alone (a free surface)

set_particles(delassus, volume, position, velocity, omega, contact=None)

Set particle/body data from DEM arrays and compute porosity/drag data.

Parameters

delassusNDArray[np.float64]

Particle Delassus operators, shape (n_particles, dim, dim)

volumeNDArray[np.float64]

Particle volumes/areas, shape (n_particles,)

positionNDArray[np.float64]

Particle positions, shape (n_particles, dim)

velocityNDArray[np.float64]

Particle velocities, shape (n_particles, dim)

omegaNDArray[np.float64]

Particle angular velocities, shape (n_particles, dim)

contactNDArray[np.float64]

Particle contact forces, shape (n_particles, dim)

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_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, one value per pressure dof in the order of pressure_index() (so not nodal values with pressure_enrichment).

Parameters

presNDArray[np.float64]

Pressure array

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_surface_tension_theta(theta)

Make the membrane geometry semi-implicit: the tensions act on the position x + theta dt u (the mesh is then moved with u) instead of x, linearised about x: the residual gets theta dt H u and the jacobian theta dt H, with H the Hessian of sigma L, sigma/L (I - t t^T) on the two diagonal blocks of an edge and its opposite off them. 0 is explicit, bounded by dt < sqrt(rho h^3/sigma); from 1/2 on it is stable for any dt, a capillary mode of pulsation omega losing about theta pi omega dt per period. 1/2 damps the least but leaves the mesh-scale modes undamped (remeshing noise can then grow), 1 damps them.

Parameters

thetafloat

Implicitness of the membrane geometry, >= 0

set_surface_tensions(sigma_ij)

Change the membrane tensions. Every edge between materials i and j has the energy sigma_ij L (L its length); the velocity equations of its two nodes get the derivative of that energy with respect to their position, -+ sigma_ij t (t the unit tangent). At equilibrium the pressure jump across a closed interface balances it: the Laplace law, without any curvature. The last row and column are the tensions with a solid, which the boundaries declared by set_wall_tension() use.

Parameters

sigma_ijNDArray[np.float64]

Symmetric (n_materials+1, n_materials+1) matrix, zero diagonal, non-negative

set_symmetry(tag)

A symmetry plane: no flux, no traction (a free-slip wall). set_weak_boundary(tag, ‘zero’).

Parameters

tagstr

Boundary physical-group name

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. Not available with the membrane model, whose materials are set_materials() labels: set the density and viscosity of each material directly.

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_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(tag, velocity=None)

A wall: no fluid crosses it, (u - u_m) . n = 0 (a zero flux on a fixed mesh). With a velocity, the fluid sticks to it (Nitsche, with the viscous traction); without, it slips freely (no traction). A wall takes no pressure: the pressure level comes from set_mean_pressure() or an outflow. set_weak_boundary(tag, ‘zero’, velocity, nitsche=velocity is not None, viscous_traction=velocity is not None).

Parameters

tagstr

Boundary physical-group name

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

Velocity of the wall, list of dim components (Real or BoundaryComponentFn) or nodal values; its normal part must be the one of the mesh, u_m . n; None for a free-slip wall

set_wall_tension(tag, sigma=None)

Declare the energy of the materials along a boundary: each of its edges, in an element of material i, has the energy sigma[i] L, treated as the membranes of set_surface_tensions(). Where an interface between i and j meets the boundary, the difference of the two wall tensions balances the tangential part of sigma_ij: the Young law sigma_ij cos(theta_i) = sigma[j] - sigma[i], with no prescribed angle (the contact line moves only if the tangential velocity is free there). A boundary the interface meets must be declared (zeros for an outflow, which carries no energy); an undeclared one raises an error at the next solve.

Parameters

tagstr

Boundary physical-group name

sigmaNDArray[np.float64]

Tension of each material with this boundary, shape (n_materials,). None: the solid of sigma_ij (its last row).

set_weak_boundary(tag, normal_flux, velocity=None, nitsche=False, pressure=None, viscous_traction=False, concentration=None, superficial_velocity=False)

Set a weak boundary condition, every term chosen explicitly (set_wall(), set_inflow(), set_outflow() and set_symmetry() are the usual cases). Combinations that leave a choice open raise an error saying what is missing.

Parameters

tagstr

Boundary physical-group name

normal_fluxstr

The normal flux in the continuity equation: ‘zero’ (no fluid crosses the boundary, (u - u_m) . n = 0, u_m the mesh velocity), ‘data’ (the flux of the velocity data, u_G . n; velocity required) or ‘free’ (the flux of the solution, u . n, an outflow)

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

Velocity data u_G, list of dim components (Real or BoundaryComponentFn) or nodal values, or None; used by normal_flux=’data’, nitsche and the upwinding of the advection

nitschebool

Impose the velocity u = u_G weakly (Nitsche penalty); velocity required

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

Normal traction: the pressure p_G (Real or BoundaryComponentFn), or None; not with normal_flux=’zero’ (a pressure level comes from set_mean_pressure() or an outflow)

viscous_tractionbool

Keep the viscous traction -mu (grad u + grad u^T) n on the boundary (the natural term of the viscous term); without it the boundary is traction free

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

Concentration/alpha of the fluid entering, for the two-fluid CSF transport, or None; not with normal_flux=’zero’

superficial_velocitybool

If True, velocity is the superficial velocity (fluid fraction times the fluid velocity, the flux per unit area through the boundary) instead of the fluid velocity; velocity required

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): the pressure, then the velocity components, then the concentration, those the boundary was given

sigma_ij()

Return a copy of the membrane tensions, materials then the solid. Returns ——-

sigma_ijNDArray[np.float64]

Tension matrix, shape (n_materials+1, n_materials+1)

Return type:

ndarray[Any, dtype[float64]]

solid_fraction()

Return the solid fraction of the grains on P1DG nodes, shape (n_elements, dim+1): the L2 projection of the grain indicator, whose mean over an element is the exact grain/element overlap over its volume. Returns ——-

solid_fractionRealDeviceArray

Solid fraction array

Return type:

RealDeviceArray

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)

surface_tension_model()

Return the surface tension model, csf or membrane (see create()). Returns ——-

modelstr

Surface tension model

Return type:

str

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=True, 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='', pressure_enrichment=False, surface_tension_model='csf', sigma_ij=None, sigma=0.0, surface_tension_theta=0.0)

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 drag in stabilisation term

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.

pressure_enrichmentbool

Pressure continuous in each material but discontinuous between them: every node gets one pressure dof per material of its elements (set_materials()). 2D and P1P1 only.

surface_tension_modelstr

csf (the capillary force of set_two_fluid_properties()) or membrane (a tension on the edges between materials, requires pressure_enrichment)

sigma_ijNDArray[np.float64]

Membrane only: surface tensions, a symmetric (n_materials+1, n_materials+1) matrix with a zero diagonal, the materials then a solid (last index, see set_wall_tension()). Fixes n_materials.

sigmafloat

Membrane only, when sigma_ij is not given: the tension between the two materials of a two-material problem (sigma_ij = [[0, sigma, 0], [sigma, 0, 0], [0, 0, 0]])

surface_tension_thetafloat

Membrane only, see set_surface_tension_theta()

Pointer

alias of LP__Structure

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

Adapt the mesh, projecting the current solution onto the new mesh. With pressure_enrichment only the velocity is projected: the materials of the new mesh are given by materials (or by set_materials() before the next solve), and the pressure restarts from 0.

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

materialsNDArray[np.float64]

pressure_enrichment only: the material of each new element, see set_materials()

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(solution, solution_old, dt)

Assemble the linear system of a Newton step of implicit_euler() at a given state, without changing the problem: the residual R and its jacobian J, as the solver gets them (strong boundaries and constraints included, one row per constraint after the n_dof rows of the solution). The stabilisation parameters are those of the state. The state of the problem (solution, old solution) is restored afterwards.

Return type:

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

Parameters

solutionNDArray[np.float64]

State, shape (n_dof,)

solution_oldNDArray[np.float64]

State at the beginning of the time step, shape (n_dof,)

dtfloat

Time step

Returns

row_ptrNDArray[np.int32]

CSR row pointers of J, shape (n_rows + 1,)

columnsNDArray[np.int32]

CSR column indices of J, shape (nnz,)

valuesNDArray[np.float64]

CSR values of J, shape (nnz,)

residualNDArray[np.float64]

R, shape (n_rows,)

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

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

dof_map_version()

Return the number of times the dof maps have been rebuilt (a new mesh, new materials): it changes exactly when the numbering of the solution, and so the linear system, does. Returns ——-

versionint

Dof map version

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: the volume weighted mean of the element gradients around each node. With pressure_enrichment the pressure gradient at a node between materials mixes those of both. 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

get_bodies_csr_force()

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

forcesRealDeviceArray

CSR body force array

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: the value of every pressure dof, in the order of pressure_index(). With pressure_enrichment these are not nodal values: a node between materials has one per material, see pressure_dg(). 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

global_map()

Return the solution-vector index of every local dof of every element, fields one after the other (the velocity components, then the pressure). Returns ——-

mapNDArray[np.int32]

Global dof map, shape (n_elements, local_size)

Return type:

ndarray[Any, dtype[int32]]

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)

materials()

Return a copy of the material of every element (see set_materials()). Returns ——-

labelsNDArray[np.int32]

Material of each element, shape (n_elements,)

Return type:

ndarray[Any, dtype[int32]]

membrane_energy()

Return the membrane energy, the sum of sigma L over the edges between materials and the edges of the boundaries of set_wall_tension(), on the current coordinates. Returns ——-

energyfloat

Membrane energy

Return type:

float

membrane_energy_gradient()

Return the derivative of membrane_energy() with respect to the position of every node: the explicit membrane term of the velocity equations (the residual, not a force: its opposite pulls the nodes). Returns ——-

gradientNDArray[np.float64]

Energy gradient, shape (n_vel_nodes, dim)

Return type:

ndarray[Any, dtype[float64]]

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_materials()

Return the number of materials the labels of set_materials() range over (from sigma_ij, 2 otherwise). Returns ——-

nint

Number of materials

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 prescribed pressure jump field, on P1DG nodes. Not available with pressure_enrichment, where the jump is an unknown. 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]]

pressure_dg()

Return the pressure at the nodes of every element, which keeps the jump between materials of pressure_enrichment. Returns ——-

pressureNDArray[np.float64]

Element pressure, shape (n_elements, local pressure size)

Return type:

ndarray[Any, dtype[float64]]

pressure_enrichment()

Return whether the pressure is enriched (see create()). Returns ——-

enrichedbool

True with pressure_enrichment

Return type:

bool

pressure_index()

Return the solution-vector index of every pressure dof, in increasing order: the node order for the first n_nodes of a P1 pressure, then, with pressure_enrichment, the copies at the nodes between materials. Returns ——-

idxNDArray[np.int32]

Pressure dof indices, shape (n_pressure_dofs,)

Return type:

ndarray[Any, dtype[int32]]

pressure_node_volume()

Return the integration weight of every pressure dof, in the order of pressure_index(): the integral of its shape function over the elements that use it (with pressure_enrichment, those of its material). With a continuous P1 pressure it is node_volume(). These are the weights of the mean pressure constraint. Returns ——-

volumeNDArray[np.float64]

Weight of each pressure dof, shape (n_pressure_dofs,)

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.

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_inflow(tag, velocity, concentration=None, superficial_velocity=False)

An inflow: the velocity is imposed (Nitsche, with the viscous traction) and its flux u_G . n enters the continuity equation. set_weak_boundary(tag, ‘data’, velocity, nitsche=True, viscous_traction=True, concentration=concentration, superficial_velocity=superficial_velocity).

Parameters

tagstr

Boundary physical-group name

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

Velocity data, list of dim components (Real or BoundaryComponentFn) or nodal values

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

Concentration/alpha of the fluid entering, or None

superficial_velocitybool

If True, velocity is the superficial velocity (fluid fraction times the fluid velocity, the flux per unit area)

set_materials(labels)

Set the material of every element, the only definition of the phases with pressure_enrichment: the pressure is continuous inside a material and discontinuous across the mesh edges between two of them. New labels rebuild the pressure dofs: the velocity is kept, the pressure is reset to 0 (a time step gives it back) and the linear system is rebuilt; unchanged labels change nothing. They are reset to 0 by set_mesh() and must be given again after adapt_mesh().

Parameters

labelsNDArray[np.float64]

Material of each element, integer values in [0, n_materials), shape (n_elements,); any other value (fractional, NaN, inf, out of range) is refused before anything changes

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_outflow(tag, pressure=0.0, viscous_traction=False)

An outflow, or a free surface: the flux u . n is free and the pressure is imposed. set_weak_boundary(tag, ‘free’, pressure=pressure, viscous_traction=viscous_traction).

Parameters

tagstr

Boundary physical-group name

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

Pressure p_G (Real, nodal values or BoundaryComponentFn)

viscous_tractionbool

Keep the viscous traction -mu (grad u + grad u^T) n (the natural outflow of the viscous term); False: the normal stress is the pressure alone (a free surface)

set_particles(delassus, volume, position, velocity, omega, contact=None)

Set particle/body data from DEM arrays and compute porosity/drag data.

Parameters

delassusNDArray[np.float64]

Particle Delassus operators, shape (n_particles, dim, dim)

volumeNDArray[np.float64]

Particle volumes/areas, shape (n_particles,)

positionNDArray[np.float64]

Particle positions, shape (n_particles, dim)

velocityNDArray[np.float64]

Particle velocities, shape (n_particles, dim)

omegaNDArray[np.float64]

Particle angular velocities, shape (n_particles, dim)

contactNDArray[np.float64]

Particle contact forces, shape (n_particles, dim)

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_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, one value per pressure dof in the order of pressure_index() (so not nodal values with pressure_enrichment).

Parameters

presNDArray[np.float64]

Pressure array

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_surface_tension_theta(theta)

Make the membrane geometry semi-implicit: the tensions act on the position x + theta dt u (the mesh is then moved with u) instead of x, linearised about x: the residual gets theta dt H u and the jacobian theta dt H, with H the Hessian of sigma L, sigma/L (I - t t^T) on the two diagonal blocks of an edge and its opposite off them. 0 is explicit, bounded by dt < sqrt(rho h^3/sigma); from 1/2 on it is stable for any dt, a capillary mode of pulsation omega losing about theta pi omega dt per period. 1/2 damps the least but leaves the mesh-scale modes undamped (remeshing noise can then grow), 1 damps them.

Parameters

thetafloat

Implicitness of the membrane geometry, >= 0

set_surface_tensions(sigma_ij)

Change the membrane tensions. Every edge between materials i and j has the energy sigma_ij L (L its length); the velocity equations of its two nodes get the derivative of that energy with respect to their position, -+ sigma_ij t (t the unit tangent). At equilibrium the pressure jump across a closed interface balances it: the Laplace law, without any curvature. The last row and column are the tensions with a solid, which the boundaries declared by set_wall_tension() use.

Parameters

sigma_ijNDArray[np.float64]

Symmetric (n_materials+1, n_materials+1) matrix, zero diagonal, non-negative

set_symmetry(tag)

A symmetry plane: no flux, no traction (a free-slip wall). set_weak_boundary(tag, ‘zero’).

Parameters

tagstr

Boundary physical-group name

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. Not available with the membrane model, whose materials are set_materials() labels: set the density and viscosity of each material directly.

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_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(tag, velocity=None)

A wall: no fluid crosses it, (u - u_m) . n = 0 (a zero flux on a fixed mesh). With a velocity, the fluid sticks to it (Nitsche, with the viscous traction); without, it slips freely (no traction). A wall takes no pressure: the pressure level comes from set_mean_pressure() or an outflow. set_weak_boundary(tag, ‘zero’, velocity, nitsche=velocity is not None, viscous_traction=velocity is not None).

Parameters

tagstr

Boundary physical-group name

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

Velocity of the wall, list of dim components (Real or BoundaryComponentFn) or nodal values; its normal part must be the one of the mesh, u_m . n; None for a free-slip wall

set_wall_tension(tag, sigma=None)

Declare the energy of the materials along a boundary: each of its edges, in an element of material i, has the energy sigma[i] L, treated as the membranes of set_surface_tensions(). Where an interface between i and j meets the boundary, the difference of the two wall tensions balances the tangential part of sigma_ij: the Young law sigma_ij cos(theta_i) = sigma[j] - sigma[i], with no prescribed angle (the contact line moves only if the tangential velocity is free there). A boundary the interface meets must be declared (zeros for an outflow, which carries no energy); an undeclared one raises an error at the next solve.

Parameters

tagstr

Boundary physical-group name

sigmaNDArray[np.float64]

Tension of each material with this boundary, shape (n_materials,). None: the solid of sigma_ij (its last row).

set_weak_boundary(tag, normal_flux, velocity=None, nitsche=False, pressure=None, viscous_traction=False, concentration=None, superficial_velocity=False)

Set a weak boundary condition, every term chosen explicitly (set_wall(), set_inflow(), set_outflow() and set_symmetry() are the usual cases). Combinations that leave a choice open raise an error saying what is missing.

Parameters

tagstr

Boundary physical-group name

normal_fluxstr

The normal flux in the continuity equation: ‘zero’ (no fluid crosses the boundary, (u - u_m) . n = 0, u_m the mesh velocity), ‘data’ (the flux of the velocity data, u_G . n; velocity required) or ‘free’ (the flux of the solution, u . n, an outflow)

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

Velocity data u_G, list of dim components (Real or BoundaryComponentFn) or nodal values, or None; used by normal_flux=’data’, nitsche and the upwinding of the advection

nitschebool

Impose the velocity u = u_G weakly (Nitsche penalty); velocity required

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

Normal traction: the pressure p_G (Real or BoundaryComponentFn), or None; not with normal_flux=’zero’ (a pressure level comes from set_mean_pressure() or an outflow)

viscous_tractionbool

Keep the viscous traction -mu (grad u + grad u^T) n on the boundary (the natural term of the viscous term); without it the boundary is traction free

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

Concentration/alpha of the fluid entering, for the two-fluid CSF transport, or None; not with normal_flux=’zero’

superficial_velocitybool

If True, velocity is the superficial velocity (fluid fraction times the fluid velocity, the flux per unit area through the boundary) instead of the fluid velocity; velocity required

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): the pressure, then the velocity components, then the concentration, those the boundary was given

sigma_ij()

Return a copy of the membrane tensions, materials then the solid. Returns ——-

sigma_ijNDArray[np.float64]

Tension matrix, shape (n_materials+1, n_materials+1)

Return type:

ndarray[Any, dtype[float64]]

solid_fraction()

Return the solid fraction of the grains on P1DG nodes, shape (n_elements, dim+1): the L2 projection of the grain indicator, whose mean over an element is the exact grain/element overlap over its volume. Returns ——-

solid_fractionRealDeviceArray

Solid fraction array

Return type:

RealDeviceArray

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)

surface_tension_model()

Return the surface tension model, csf or membrane (see create()). Returns ——-

modelstr

Surface tension model

Return type:

str

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(iterative_refinement=0)

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

Return type:

ILinearSystem

Parameters

iterative_refinementint

Maximum number of iterative refinement steps after the factorisation (MKL iparm[7]); 0 keeps MKL’s default, two steps when pivots were perturbed. The membrane cases of dev used 5.

Returns

solverILinearSystem

Solver object for FluidProblem.create()

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