CFD Domain Pack¶
pip install ontic-engine[cfd]
Modules¶
| Module | Description |
|---|---|
ontic.cfd.euler_1d |
1D compressible Euler equations (shock tubes) |
ontic.cfd.euler_2d |
2D compressible Euler (oblique shocks, wedges) |
ontic.cfd.godunov |
Godunov-type Riemann solvers |
ontic.cfd.navier_stokes |
Viscous Navier-Stokes solver |
ontic.cfd.les |
Large Eddy Simulation |
ontic.cfd.hybrid_les |
Hybrid LES methods |
ontic.cfd.multi_objective |
Multi-objective optimization |
ontic.cfd.pure_qtt_ops |
Pure QTT CFD operators |
API Reference¶
Ontic CFD — public re-exports.
Every symbol listed in __all__ is importable directly from
ontic.cfd so that downstream code (tests, experiments, docs,
CUDA bridge) can use short-form imports like::
from ontic.cfd import Euler1D, exact_riemann, sod_shock_tube_ic
Bases: Enum
Boundary condition types for 1D Euler equations.
1D Euler equation solver using tensor networks.
Discretization: Finite volume with N cells on domain [x_min, x_max].
The solution is stored as an MPS where each site represents a cell, and the physical dimension is 3 (for ρ, ρu, E).
For now, we use a simple product state (χ=1) representation, which is equivalent to classical FVM. The tensor network structure enables future extensions: - Entangled representations for multi-scale features - MPO-based flux operators - Automatic adaptivity via bond dimension
Attributes: N: Number of grid cells x_min, x_max: Domain bounds gamma: Ratio of specific heats cfl: CFL number for time stepping dx: Cell width x_cell: Cell center coordinates x_face: Face coordinates state: Current EulerState t: Current simulation time
Example: >>> solver = Euler1D(N=200, x_min=0.0, x_max=1.0) >>> state = sod_shock_tube_ic(solver.x_cell) >>> solver.set_initial_condition(state) >>> for _ in range(100): ... solver.step() >>> print(f"Final time: {solver.t:.4f}")
Raises: ValueError: If state dimensions don't match grid size RuntimeError: If CFL condition is violated
References: .. [1] Toro, E.F. "Riemann Solvers and Numerical Methods for Fluid Dynamics", 3rd ed., Springer, 2009. .. [2] LeVeque, R.J. "Finite Volume Methods for Hyperbolic Problems", Cambridge University Press, 2002.
Set boundary conditions.
Args: left: Left boundary condition type right: Right boundary condition type
Returns: self for method chaining
Compute numerical flux at cell interface using Rusanov (local Lax-Friedrichs).
F_{i+1/2} = ½(F_L + F_R) - ½ λ_max (U_R - U_L)
Args: U_L: Left state (batch, 3) U_R: Right state (batch, 3)
Returns: Flux (batch, 3)
Advance solution by one time step.
Uses first-order forward Euler in time with Rusanov flux in space.
Args: dt: Time step (computed from CFL if None)
Returns: Actual dt used
Solve to final time.
Args: t_final: Final time callback: Optional callback(solver, t) called periodically callback_interval: Interval for callback (0 = every step)
Returns: List of (time, state) snapshots
Convert current state to MPS representation.
Each site is a grid cell with physical dimension 3 (ρ, ρu, E).
For chi_max=1 (product state), this is equivalent to classical FVM. Larger chi enables entanglement for multi-scale representation.
Args: chi_max: Maximum bond dimension
Returns: MPS representation of state
classmethod
¶
Create solver from MPS state.
Args: mps: MPS with physical dimension 3 x_min, x_max: Domain bounds gamma: Ratio of specific heats cfl: CFL number
Returns: Euler1D solver with state from MPS
dataclass
¶
2D Euler equations solver using finite volume method.
Uses dimensional splitting (Strang splitting) to extend 1D schemes to 2D with second-order accuracy in time.
Strang splitting: L(Δt) = Lx(Δt/2) Ly(Δt) Lx(Δt/2)
Attributes: Nx: Number of cells in x-direction Ny: Number of cells in y-direction x: Cell-center x coordinates y: Cell-center y coordinates dx: Grid spacing in x dy: Grid spacing in y gamma: Ratio of specific heats state: Current flow state
set_boundary_conditions(left=BCType.OUTFLOW, right=BCType.OUTFLOW, bottom=BCType.OUTFLOW, top=BCType.OUTFLOW, inflow_state=None)
Configure boundary conditions.
Compute stable time step based on CFL condition.
For 2D: dt ≤ CFL / (|u|/dx + |v|/dy + a*sqrt(1/dx² + 1/dy²))
Advance solution by one time step using Strang splitting.
Strang splitting achieves second-order accuracy: U^{n+1} = Lx(dt/2) Ly(dt) Lx(dt/2) U^n
Args: dt: Time step (computed from CFL if None) cfl: CFL number for automatic dt
Returns: Actual time step used
Integrate to final time.
Args: t_final: Target simulation time cfl: CFL number max_steps: Maximum number of steps verbose: Print progress
Returns: Dictionary with simulation info
Bases: Enum
Available LES subgrid-scale models.
dataclass
¶
State for LES subgrid quantities.
dataclass
¶
State variables for hybrid RANS-LES models.
Bases: Enum
Available hybrid RANS-LES models.
Bases: Enum
Multi-objective optimization algorithms.
dataclass
¶
Configuration for multi-objective optimization.
dataclass
¶
Result from multi-objective optimization.
Multi-objective optimization driver.
Finds Pareto-optimal designs trading off multiple objectives.
Simulated Binary Crossover (SBX).
Creates two offspring from two parents with distribution controlled by η_c parameter.
Binary tournament selection based on rank and crowding distance.
Run weighted sum scalarization.
Samples different weight combinations to approximate Pareto front. Note: Can miss concave parts of the front.
dataclass
¶
Specification for an objective function.
dataclass
¶
A single solution on the Pareto front.
Coupled 2D compressible Navier-Stokes solver.
Combines the inviscid Euler solver with viscous terms using explicit operator splitting.
Features: - HLLC Riemann solver for inviscid fluxes - Sutherland's law for temperature-dependent viscosity - Automatic timestep selection (min of CFL and viscous) - Boundary layer resolution awareness
Example: >>> config = NavierStokes2DConfig(Nx=128, Ny=64, Lx=1.0, Ly=0.5) >>> ns = NavierStokes2D(config) >>> state = flat_plate_ic(config) >>> result = ns.solve(state, t_final=0.1)
Compute stable timestep from CFL and viscous limits.
Args: state: Current flow state
Returns: Stable timestep
Apply viscous update via explicit Euler.
∂U/∂t = ∇·F_v => U^{n+1} = U^n + dt * ∇·F_v
Args: state: Current state dt: Timestep
Returns: Updated state with viscous effects
Advance solution one timestep using operator splitting.
Strang-type splitting: 1. Half viscous step 2. Full inviscid step (Euler) 3. Half viscous step
Args: state: Current state dt: Timestep
Returns: Updated state
Solve Navier-Stokes equations to final time.
Args: initial_state: Initial condition t_final: Final simulation time callback: Optional callback(state, t, step) returning True to stop max_steps: Maximum timesteps
Returns: NavierStokes2DResult with final state and diagnostics
Compute wall heat transfer and skin friction.
Args: state: Current state wall_j: j-index of wall (default 0 = bottom)
Returns: Dictionary with tau_wall, q_wall, Cf, St
dataclass
¶
Configuration for Navier-Stokes solver.
Immersed boundary method for solid bodies in flow domain.
Uses ghost-cell method: cells inside the body are set to mirror the flow outside, with reflected velocity to enforce no-penetration condition.
Accepts any :class:~ontic.cfd.sdf.SDFGeometry implementation
(circle, airfoil, fin-array, etc.) or the legacy
:class:WedgeGeometry for backward compatibility.
Attributes: geometry: Body geometry (SDFGeometry or WedgeGeometry) X, Y: Grid coordinate meshes mask: Boolean tensor, True inside solid
Apply immersed boundary condition.
Sets ghost cell values by reflecting from image points.
Args: U: Conservative variables (4, Ny, Nx)
Returns: Modified conservative variables
Extract flow data on the body surface.
Args: U: Conservative variables (4, Ny, Nx) gamma: Ratio of specific heats
Returns: Dictionary with surface pressure, Mach number, etc.
dataclass
¶
Definition of a 2D wedge body.
The wedge has its leading edge at (x_le, y_le) and extends in the +x direction with symmetric half-angles.
╱
╱ θ (half-angle)
╱──────────────── ╲──────────────── ╲ θ ╲
Attributes: x_leading_edge: x-coordinate of leading edge y_leading_edge: y-coordinate of leading edge (wedge centerline) half_angle: Wedge half-angle in radians length: Wedge length in x-direction
Compute upper and lower surface y-coordinates.
Args: x: x-coordinates (tensor)
Returns: (y_upper, y_lower): Surface coordinates
Check if points are inside the wedge body.
Args: x, y: Coordinate tensors (same shape)
Returns: Boolean tensor, True where (x, y) is inside solid
Compute signed distance to nearest wedge surface.
Positive outside, negative inside.
Args: x, y: Coordinate tensors
Returns: Signed distance tensor
Compute outward-pointing surface normal at nearest surface point.
Args: x, y: Coordinate tensors
Returns: (nx, ny): Normal vector components
Bases: SDFGeometry
Annular region between two concentric circles.
The SDF is negative in the solid regions (inside the inner circle OR outside the outer circle), i.e. the fluid region (the annular gap) has positive SDF.
For immersed-boundary use: the inner cylinder is the solid body; the outer cylinder is the domain boundary. This SDF represents the inner cylinder only.
Bases: SDFGeometry
Axis-aligned ellipse at (cx, cy) with semi-axes a (x) and b (y).
Uses a smooth approximation; exact SDF for an ellipse requires iterative root-finding. The approximation is excellent for immersed-boundary purposes (error < grid spacing).
Bases: SDFGeometry
Array of rectangular fins extending from a base plate.
Fins extend in +y from base_y. The base plate itself is included as part of the SDF.
Bases: SDFGeometry
CSG union of multiple SDFGeometry objects.
The SDF at each point is min(sdf_1, sdf_2, …) — the standard
signed-distance union operator.
Bases: SDFGeometry
NACA 4-digit airfoil profile as an SDF.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
chord
|
float
|
Chord length (metres). |
required |
naca_code
|
str
|
4-digit NACA code, e.g. "0012", "2412". |
'0012'
|
leading_edge_x
|
float
|
Position of the leading edge. |
0.0
|
leading_edge_y
|
float
|
Position of the leading edge. |
0.0
|
aoa_deg
|
float
|
Angle of attack in degrees (positive nose-up). |
0.0
|
Bases: SDFGeometry
2D representation of a pipe with a circular bend.
The pipe centre-line enters horizontally from the left, bends through bend_angle_deg (default 90°), and exits vertically upward. The wall thickness defines an annular region around the bend.
Bases: ABC
Abstract base for 2-D signed-distance-function geometries.
Every concrete subclass must implement sdf. The default
is_inside, normal, and surface_points derive from
sdf automatically but may be overridden for efficiency.
abstractmethod
property
¶
Tight axis-aligned bounding box: (x_min, x_max, y_min, y_max).
abstractmethod
property
¶
Reference length for Reynolds-number computation (metres).
abstractmethod
¶
Signed distance from each (x, y) to the nearest surface.
Returns:
| Type | Description |
|---|---|
Tensor
|
Same shape as x / y. Negative inside, positive outside. |
Outward-pointing unit normal via autograd on sdf.
Falls back to finite-difference if autograd is unavailable.
Sample n approximately equi-spaced points on the surface.
Default implementation marches around the bounding box and uses bisection along radial rays from the centroid.
Bases: SDFGeometry
Backward-facing step geometry.
A solid block occupies the upper portion of the channel upstream of the step location, creating a sudden expansion.
::
████████████████████
████████████████████
████████████
████████████ ← step_height
───────────────────── channel floor (y=0)
^
step_x
Bases: SDFGeometry
Symmetric sharp wedge pointing into the flow (-x direction).
Leading edge at (tip_x, tip_y), extending in +x by length.
sod_shock_tube_ic(N, x_min=0.0, x_max=1.0, x_discontinuity=0.5, gamma=1.4, dtype=torch.float64, device=None)
Sod shock tube initial condition.
Classic test problem with exact solution available.
Left state (x < 0.5): ρ = 1, u = 0, p = 1 Right state (x > 0.5): ρ = 0.125, u = 0, p = 0.1
Features: rarefaction, contact, shock
exact_riemann(rho_L, u_L, p_L, rho_R, u_R, p_R, gamma=1.4, x=None, t=1.0, x0=0.5, tol=1e-08, max_iter=100)
Exact Riemann solver using Newton-Raphson iteration.
Solves for the star-state pressure p* at the contact, then constructs the full solution.
Args: rho_L, u_L, p_L: Left primitive state rho_R, u_R, p_R: Right primitive state gamma: Ratio of specific heats x: Spatial coordinates for sampling solution t: Time at which to sample solution x0: Initial discontinuity location tol: Newton tolerance max_iter: Maximum Newton iterations
Returns: (rho, u, p) sampled at x coordinates
supersonic_wedge_ic(Nx, Ny, M_inf=5.0, p_inf=1.0, rho_inf=1.0, gamma=gamma_default, dtype=torch.float64, device='cpu')
Create initial condition for supersonic flow (uniform freestream).
Args: Nx, Ny: Grid dimensions M_inf: Freestream Mach number p_inf: Freestream pressure rho_inf: Freestream density gamma: Ratio of specific heats
Returns: Euler2DState with uniform supersonic flow in +x direction
Unified interface for computing SGS viscosity.
Args: model: LES model type Velocity gradients delta: Filter width rho: Density u, v: Velocities (needed for dynamic model) **kwargs: Model-specific parameters
Returns: SGS eddy viscosity μ_sgs
Compute LES filter width Δ.
Common choices: - Δ = (dx·dy·dz)^(1/3) for 3D - Δ = (dx·dy)^(1/2) for 2D - Δ = max(dx, dy, dz) for anisotropic grids
Args: dx, dy: Grid spacings dz: Optional z spacing for 3D
Returns: Filter width scalar
Classic Smagorinsky subgrid-scale viscosity.
ν_sgs = (C_s Δ)² |S|
Args: S: Strain rate magnitude [1/s] delta: Filter width [m] rho: Density [kg/m³] C_s: Smagorinsky constant
Returns: SGS eddy viscosity μ_sgs [Pa·s]
strain_rate_magnitude(du_dx, du_dy, dv_dx, dv_dy, du_dz=None, dv_dz=None, dw_dx=None, dw_dy=None, dw_dz=None)
Compute strain rate magnitude |S| = √(2 S_ij S_ij).
S_ij = (1/2)(∂u_i/∂x_j + ∂u_j/∂x_i)
Args: Velocity gradients in 2D or 3D
Returns: |S| strain rate magnitude [1/s]
vreman_viscosity(du_dx, du_dy, dv_dx, dv_dy, delta, rho, C_v=C_V, du_dz=None, dv_dz=None, dw_dx=None, dw_dy=None, dw_dz=None)
Vreman subgrid-scale model (2004).
Based on the first invariant of the velocity gradient tensor:
ν_sgs = C_v √(B_β / (α_ij α_ij))
where α_ij = ∂u_j/∂x_i, β_ij = Δ²_m α_mi α_mj and B_β = β_11 β_22 - β_12² + β_11 β_33 - β_13² + β_22 β_33 - β_23²
Advantages: - Vanishes in many laminar flows - Works well on anisotropic grids - Simple and efficient
Args: Velocity gradients delta: Filter width [m] rho: Density [kg/m³] C_v: Vreman constant
Returns: SGS eddy viscosity μ_sgs [Pa·s]
wale_viscosity(du_dx, du_dy, dv_dx, dv_dy, delta, rho, C_w=C_W, du_dz=None, dv_dz=None, dw_dx=None, dw_dy=None, dw_dz=None)
WALE (Wall-Adapting Local Eddy-viscosity) model.
Based on the traceless symmetric part of the squared velocity gradient tensor:
S^d_ij = (1/2)(g²_ij + g²_ji) - (1/3) δ_ij g²_kk
where g_ij = ∂u_i/∂x_j and g²_ij = g_ik g_kj
Advantages: - Proper near-wall behavior (ν_sgs ~ y³) - No ad-hoc damping functions required - Zero in laminar shear flow
Args: Velocity gradients delta: Filter width [m] rho: Density [kg/m³] C_w: WALE constant
Returns: SGS eddy viscosity μ_sgs [Pa·s]
DDES delay/shielding function f_d.
f_d = 1 - tanh([C_d1 * r_d]^C_d2)
When f_d ≈ 0 (r_d large), RANS mode is preserved. When f_d ≈ 1 (r_d small), LES mode can activate.
Args: r_d: Delay function parameter c_d1, c_d2: Model constants
Returns: Delay function f_d
Original DES length scale.
l_DES = min(l_RANS, C_DES * Δ)
When Δ is smaller (in separated regions), LES mode activates.
Args: l_rans: RANS length scale (κ * d_wall) delta: LES grid scale c_des: DES coefficient
Returns: Hybrid length scale
Compute statistics on RANS vs LES content.
Args: state: Hybrid LES state
Returns: Dict with RANS/LES percentages
Main driver for hybrid RANS-LES computation.
Args: rho: Density field u: Velocity field (Nd, Nx, Ny, [Nz]) d_wall: Wall distance field grid_spacing: Tuple of (dx, dy, [dz]) tensors nu: Molecular viscosity nu_rans: RANS turbulent viscosity field model: Hybrid model to use
Returns: HybridLESState with computed quantities
Check if solution A dominates solution B.
A dominates B if: - A is no worse than B in all objectives - A is strictly better than B in at least one objective
Args: obj_a, obj_b: Objective values for solutions A and B minimize: Dict indicating which objectives to minimize
Returns: True if A dominates B
Fast non-dominated sorting (NSGA-II).
Assigns each solution to a Pareto front (rank). Rank 0 = non-dominated, Rank 1 = dominated by rank 0, etc.
Args: population: List of solutions minimize: Minimization flags per objective
Returns: List of fronts, each containing indices of solutions