Models

One essential part of a numerical simulation is to define the models and conditions for it. There are multiple features of Nassu that require some configuration like this. These configurations are gathered in the models field and define the behavior of the simulation in multiple ways.

Some of the features available in it are:

• Numerical models

• Boundary conditions

• Initialization

• Devices and precision to use

Note

For further details on theorical aspects of the models presented here, check the theory documentation.

Numerical models

When running a simulation, one of the most important aspects is the numerical models that are used. Nassu combines LBM with other models to achieve its capabilities and each model has its own configurations.

simulations:
  - name: example
    models:
      # Multiblock options
      multiblock:
        # Refinement-interface communication operator: value_interp (default,
        # Lagrava-2012 value interpolation) or conservative (conservative
        # moment-transfer pair). Defaults to value_interp
        interface_coupling: value_interp
        # Number of overlap nodes in Fine-to-coarse communication
        # Defaults to 1
        overlap_F2C: 2
        # Custom overlaps of F2C for specific levels
        # Use as [lvl_from]: overlap to use
        # Defaults to {} (no custom)
        custom_overlap_F2C:
          1: 3
          2: 3
          5: 1
        # Marks nodes inside coarse blocks used for F2C communication as unused,
        # preventing undefined behavior from invalid node values.
        # Defaults to true
        mark_nodes_as_unused: true
        # Optional adaptive checkerboard dissipation for near-inviscid
        # refinement interfaces (issue #1138). Requires models.LES. Omit to
        # disable (generated kernels byte-identical).
        interface_stabilization:
          strength: 0.2
          threshold: 0.5
      # Large eddy simulation options
      LES:
        # Model to use (Smagorinsky)
        model: Smagorinsky
        # Subgrid Smagorinsky constant
        sgs_cte: 0.17
      # LBM options
      LBM:
        # Global force to use
        F: { x: 0, y: 0, z: 0 }
        # Collision operator
        # Options are: RBGK, RRBGK, HRRBGK
        coll_oper: RRBGK
        coll_oper_params:
          # Constant sigma to use for HRRBGK, between 0.95 and 1. 
          # Defaults to 0.99
          sigma_hrrbgk: 0.99
          # Mode for HRRBGK
          # Options are: dynamic, constant
          # Defaults to constant
          mode_hrrbgk: dynamic
        # True to activate thermal model and variable theta
        thermal_model: true
        # Relaxation time (tau) value
        tau: 0.5000008125
        # Optional independent bulk (volume) viscosity zeta (lattice units).
        # When set, the trace of the non-equilibrium second moment relaxes at
        # omega_bulk = 1 / (zeta / cs^2 + 1/2) instead of the shear rate omega,
        # damping the low-Mach reduced-pressure checkerboard mode. Omit (the
        # default) to keep the near-zero intrinsic bulk viscosity; the kernels
        # are then byte-identical.
        bulk_viscosity: 0.1
        # Velocity set
        # Options are: D2Q9 for 2D, D3Q15, D3Q19, D3Q27 for 3D
        vel_set: D3Q27
      # Immersed boundary method options
      IBM:
        # Dirac delta to use
        # Options are: 3_points, 4_points
        # Defaults to 3_points
        dirac_delta: 3_points
        # For all IBM interpolation or spread kernels, the minimun value of dirac for the operation
        # to be done. 
        # Use high value to avoid running IBM on domain borders or multiblock transitions
        # Use low value to allow for IBM to run everywhere
        # Defaults to 0.99
        min_dirac_sum: 0.99
        # Number of time steps to use for force accomodation
        # The force applied linearly increases from 0% to 100% in this interval
        # Defaults to 0
        forces_accomodate_time: 500
        # Limiter to force spreading on IBM iteration. This is the limit between two iterations. 
        # A reasonable value is 1e-3 for limiting (force ~ rho*u^2, So u~0.03 yields this)
        # Defaults to 1000 (in practice, no limit)
        forces_spread_limit: 1e-3
        # Whether to reset forces between time steps or use the force from the previous
        # time step as reference
        # Defaults to true
        reset_forces: true
        # IBM configurations for bodies
        body_cfgs:
          default: {}
          # Configuration name to reference
          building_cfg:
            # Number of iterations
            # Defaults to 5
            n_iterations: 5
            # Force factor to use for force spreading.
            # The force is multiplied by this factor before being applied.
            # Defaults to 1
            forces_factor: 1
            # Debug-only model that applies a constant force to drive the body to a
            # constant velocity. A tool for exercising the IBM force path, not a
            # physical model; leave unset in production runs. Defaults to null
            # (fields explained in `constant_force_cfg`)
            debug_cte_force: null
            # Wall model to use
            wall_model:
              # Name of wall model
              # Availables one are: EqLog, EqTBL, NonEqTBL
              name: NonEqTBL
              # Reference distance for tangential velocity interpolation
              dist_ref: 2
              # Shell distance for spreading forces
              dist_shell: 0.25
              # Step to start applying wall model. Before this step
              # conventional IBM is applied
              # Defaults to 1000
              start_step: 1000
              # Parameters required by the chosen wall model
              params:
                # Roughness length
                z0: 0.0001
                # Max error for TDMA
                TDMA_max_error: 5e-06
                # Max iterations for TDMA 
                TDMA_max_iters: 50
                # Number of divisions in TDMA (max, min)
                # Must be an odd number
                TDMA_max_div: 25
                TDMA_min_div: 21
                # Target y+ for nearest wall point in TDMA discretization
                # defaults to 0.2
                TDMA_yp_target: 0.2
                # Exponential wall-clustering strength (>= 0) for the TBL TDMA
                # grid: n = dist_interp*(exp(beta*xi)-1)/(exp(beta)-1) on a
                # uniform xi. Clusters nodes toward the wall for better near-
                # wall resolution at the same division count. Default 2.0;
                # 0.0 gives a uniform grid.
                TDMA_stretch_beta: 2.0
                # NonEqTBL only - floor (> 0) on the previous-step friction
                # velocity used to compute the LaRTE adaptive EMA coefficient
                # on the streamwise pressure gradient,
                # alpha_node = u_friction_lagged / (dist_ref + u_friction_lagged).
                # Keeps alpha_node well-defined at start-up and at
                # near-stagnation nodes. Default 1e-4 (lattice); a tighter
                # case-specific guideline is one tenth of the log-law
                # equilibrium friction velocity at the freestream scale.
                NeqWM_u_friction_floor: 1.0e-4
                # NonEqTBL only - static multiplier (>= 0) on the LaRTE-
                # filtered pressure gradient before it enters the TBL TDMA.
                # Physically motivated value: m = Ma_LBM = U_LBM * sqrt(3),
                # which rescales the raw LBM `dp/ds` magnitude back to the
                # incompressible-equivalent expectation - see
                # `theory/wall_model/neq_pres_filter`. Set to 0.0 to zero the
                # streamwise pressure-gradient source seen by the TDMA
                # (NonEqTBL then behaves like the equilibrium TBL) while
                # pres_grad and pres_grad_filt keep updating - a passive
                # raw-signal sensor for spectral diagnostics.
                NeqWM_pres_grad_mult: 1.0
              # Optional flux-driven Monin-Obukhov stratified-surface-layer
              # stability correction. Omit it (the default) for the neutral
              # log-law path. The buoyancy strength B = beta*|gravity| is read
              # from the coupling scalar's `buoyancy` block, so `coupling_scalar`
              # must name a buoyant `scalar_transports` entry.
              mo_stability:
                # Scalar the buoyancy strength B = beta*|gravity| is read from
                # (must declare `buoyancy`).
                coupling_scalar: temperature
                # Prescribed surface kinematic heat flux <w'phi'>_0 (lattice
                # units), the MO driver. Positive -> unstable (fuller profile);
                # negative -> stable (more sheared); zero -> neutral.
                surface_heat_flux: 1.0e-6
                # Optional startup ramp: linearly ease the surface flux in from
                # zero (neutral) over this many coarsest-level steps. 0 (default)
                # applies it at full strength from step 0.
                ramp_steps: 0
          terrain_cfg:
            n_iterations: 3
            forces_factor: 1
            wall_model:
              name: EqLog
              dist_ref: 2.5
              dist_shell: 0.5
              # Parameters required by EqLog
              params:
                z0: 0.0001
          constant_force_cfg:
            n_iterations: 1
            forces_factor: 0.12  
            debug_cte_force:
              # Constant velocity to apply in body normal
              constant_velocity: 1e-2
              # Tangential direction force factor. 
              # 0 to not change tangential velocity, 1 to make tangential velocity zero
              tangential_force_factor: 0
              # Use body normals to make force direction
              use_normal: false
              # Direction to apply force, in case normal is not used
              # Defaults to null
              direction_apply: [0, 1, -1]

Note

The configurations for these model requires an understanding of how the solver works under the hood. So be careful when changing these values from the guideline ones.

An IBM wall model can optionally carry a mo_stability block, selecting a flux-driven Monin-Obukhov stratified-surface-layer stability correction on top of the neutral log-law. The correction is flux-driven: surface_heat_flux prescribes the surface kinematic heat flux <w'phi'>_0 (lattice units) directly, which forms the Obukhov length L. A positive flux (upward, warming) gives an unstable L < 0 and a fuller, less-sheared near-wall profile; a negative flux gives a stable L > 0 and a more-sheared profile; zero recovers the neutral log-law. coupling_scalar names the models.scalar_transports entry the buoyancy strength B = beta*|gravity| is read from (that scalar must declare a buoyancy block), so the wall model and the fluid Boussinesq force cannot disagree on beta*g; the scalar field itself is not sampled by the momentum correction. The optional ramp_steps linearly eases surface_heat_flux in from zero over the first ramp_steps coarsest-level steps (analogous to the IBM forces_accomodate_time), so the wall model starts neutral (L -> inf) and ramps into the stratified correction, removing the step-0 surface-flux shock in stratified starts; 0 (default) applies the full flux from step 0. Omitting mo_stability (the default) keeps the neutral log-law path unchanged.

Important

The IBM EqLog wall model applies the stability-corrected log law, with u* and L converged by a short Picard iteration. The empirical stability-function coefficients are hard-coded to gamma = 16, beta_s = 5 (a validation case should cite the coefficients it adopts).

Scalar transport

Nassu can solve the advection-diffusion equation (ADE) for one or more passive scalars on top of the fluid simulation. Typical applications include temperature in atmospheric and urban-scale flows, pollutant or tracer concentration, and any quantity transported by the resolved velocity field with linear diffusion.

Each entry under models.scalar_transports adds a separate scalar LBM coupled to the fluid: an independent set of populations \(g_i\) with its own velocity set, evolving under streaming and a recursive regularised collision and reading the fluid velocity at every step. The dictionary key (e.g. temperature) becomes the symbol-name scope for the scalar’s macroscopics, so two scalars produce disjoint variables (temperature_phi, pollutant_phi, …) and never collide in the generated CUDA.

simulations:
  - name: example
    models:
      scalar_transports:
        # Each named entry adds one scalar advection-diffusion LBM.
        temperature:
          # Lattice for the scalar populations g_i.
          # Options: D3Q7 (production default in 3-D), D2Q5 (2-D debug
          # analogue), D3Q15, D3Q19, D3Q27 (higher cost, better isotropy
          # at high Peclet number).
          velocity_set: D3Q7
          # Scalar collision operator. Only RRBGK is supported.
          collision_operator: RRBGK
          # When to advance the scalar. Defaults to (0, 0) which means
          # the scalar is not stepped; set start_step / end_step to enable.
          interval: { start_step: 0, end_step: 0 }
          adv_diff_equation:
            # Molecular diffusivity D in lattice units.
            # When LES is active, the scalar diffusivity is augmented by
            # nu_SGS / Sc_t (turbulent Schmidt number).
            D: 1.0e-5
            # Source / sink term S(x, y, z, t, phi). "0" for a passive
            # scalar with no source.
            S: "0"
          # Turbulent Schmidt (or Prandtl) number. With LES active, the
          # node-local eddy diffusivity nu_SGS / Sc_t adds to D in the
          # collision and in the multiblock level transfers. Runtime
          # tunable (no kernel regeneration); ignored without LES.
          Sc_t: 0.7
          # Optional Boussinesq buoyancy: the scalar exerts the body
          # force F_buoy = -rho0 * beta * (phi - phi_ref) * gravity on
          # the fluid through the Guo source term. `gravity` points
          # along gravity in lattice units (e.g. [0, 0, -g_lbm] with z
          # up), so a parcel with phi > phi_ref rises. Compile-time
          # parameters; at most one scalar per simulation may declare
          # buoyancy. Storage cost is one float per node: only the
          # signed magnitude along gravity is stored, and the force is
          # projected onto the (compile-time) gravity direction at use.
          buoyancy:
            beta: 2.0e-3
            phi_ref: 0.0
            gravity: [0.0, 0.0, -1.0e-5]
            rho0: 1.0
            # Optional startup ramp: linearly ease the buoyancy force in
            # from zero over this many coarsest-level steps (analogous to
            # the IBM `forces_accomodate_time`), removing the step-0
            # buoyancy shock. 0 (default) applies it at full strength.
            ramp_steps: 0
          # Initial field expression; evaluated once at setup. May
          # reference (x, y, z) in lattice coordinates.
          initial_field: "0"
          # Per-scalar boundary conditions. Aliases live in
          # `AllScalarBCSchemesAliases`:
          # ScalarRegularizedDirichlet (prescribed phi),
          # ScalarRegularizedNeumann (prescribed flux, J_w=0 is
          # the zero-flux adiabatic wall),
          # ScalarRegularizedRobin (convective heat loss, h and
          # phi_inf), ScalarUniformInlet (equilibrium from prescribed
          # phi; the advection velocity comes from the fluid).
          # Periodicity is inherited from `models.BC.periodic_dims`.
          # `pos`, `system` and `wall_normal` follow the fluid `BC_map`
          # grammar: a cardinal alias or a boolean predicate over x, y, z,
          # optionally written in a `domain.systems` frame, with the
          # outward cardinal `wall_normal` stated explicitly.
          BC:
            BC_map:
              - pos: S
                BC: ScalarRegularizedNeumann
                J_w: 0.0
                wall_normal: S
              - pos: N
                BC: ScalarRegularizedNeumann
                wall_normal: N
                params:
                  J_w: 0.0
              - pos: E
                BC: ScalarRegularizedRobin
                wall_normal: E
                params:
                  h: 0.0058
                  phi_inf: 0.0
              # Emitting footprint on the ground: an interior surface
              # selected by a predicate instead of a whole domain face.
              - pos: "(z <= 0) & (x >= 120) & (x <= 180) & (y >= 40) & (y <= 90)"
                BC: ScalarRegularizedDirichlet
                wall_normal: B
                order: 1
                params:
                  phi_w: 1.0
          # Volumetric scalar sources (region emission): add a constant
          # rate to the scalar collision on every node inside a region,
          # complementing the `adv_diff_equation.S` expression. The region
          # is placed by the same `pos` a surface BC takes: a cardinal
          # alias, a predicate over x, y, z, or a typed shape.
          # The `rate` is compile-time; all regions of one scalar must
          # share the same `rate`, and the region geometry is shared across
          # scalars (each scalar bakes its own rate).
          source_regions:
            - pos: "(z >= 30) & (z <= 38)"
              rate: 1.0e-3

The minimal velocity sets D2Q5 and D3Q7 (rest plus axis neighbours) carry only the rank-1 Hermite tensor, which is exactly what the first-order scalar equilibrium needs. Larger sets (D3Q15, D3Q19, D3Q27, also D2Q9 in 2-D) are accepted for advection-dominated regimes where directional bias on D3Q7 becomes visible, at the cost of more populations per node.

Multiple scalars can coexist with different velocity sets - for instance temperature on D3Q7 and a high-Peclet pollutant on D3Q27 share the same fluid solver:

        scalar_transports:
          temperature:
            velocity_set: D3Q7
            # ... (other fields)
          pollutant:
            velocity_set: D3Q27
            # ... (other fields)

A scalar can also emit a volumetric source over a region of the domain through source_regions: each entry adds a constant rate to the scalar collision on every node inside its predicate region (a pollutant-emitting region), on top of the global adv_diff_equation.S expression and any surface emission. The region is placed by the same pos the surface BCs take, in any of its spellings - a cardinal alias, a boolean predicate or a typed shape (see Placing a region with a typed shape) - and the rate is baked into the kernel: all regions of one scalar must share a single rate, and the region geometry is shared across scalars so several species can emit in the same region with their own rates.

The prescribed value of a ScalarUniformInlet (phi_inlet) or a ScalarRegularizedDirichlet (phi_w) accepts either a constant float or an (x, y, z) equation string for a height-varying (stratified) inlet, e.g. phi_inlet: "0.05 * (z / 100.0) ** 0.25".

Important

For the equation (z-varying) form of phi_inlet / phi_w, a non-constant expression is baked per-node into the scalar BC kernel and evaluated at each node’s physical coordinates. A plain float is baked as a constant.

Energy field (DDF backend)

An entry can instead carry kind: energy (the default is kind: scalar) to transport energy with the double-distribution-function (DDF) energy backend. The energy field is a transported field on the same advection-diffusion path as a scalar: it streams a second distribution \(g_i\) whose conserved zeroth moment is the volumetric energy density \(E = \rho h\), advected and diffused by the fluid velocity. The two scalar diffusion knobs take their energy meaning - alpha is the thermal diffusivity \(\mu/(\rho_0\,\mathrm{Pr})\) (the slot the scalar uses for D) and Pr_t is the turbulent Prandtl number (the slot the scalar uses for Sc_t). An energy field takes no adv_diff_equation, buoyancy or source_regions: its source is the intrinsic deviation-cancelling correction, not a user term.

        scalar_transports:
          energy:
            kind: energy
            # The energy equilibrium is second order, so the set must be
            # rank-4 isotropic: D2Q9 / D3Q15 / D3Q19 / D3Q27. The minimal
            # D2Q5 / D3Q7 sets are rejected for energy.
            velocity_set: D3Q27
            collision_operator: RRBGK
            interval: { start_step: 0, end_step: 0 }
            # Thermal diffusivity mu/(rho_0 Pr) (constant).
            alpha: 1.0e-3
            # Turbulent Prandtl number (LES coupling), runtime tunable.
            Pr_t: 0.7
            # Initial energy density field, may reference (x, y, z).
            initial_field: "exp(-((x-8)*(x-8) + (y-8)*(y-8) + (z-8)*(z-8)) / 8.0)"

For the underlying numerics (DDF formulation, RRBGK collision, LES coupling, multiblock rescaling, boundary conditions), see the scalar transport theory pages. The DDF energy closure (conserved variable, advection-diffusion equilibrium, deviation correction, variable-density conduction) is documented in the DDF energy theory page.

Weakly-compressible thermal coupling

The optional models.energy block activates the weakly-compressible thermal coupling, Tier 2 of the thermal modelling hierarchy. The shared finite-difference energy field feeds a temperature deviation \(\theta = T/T_\mathrm{ref} - 1\) into the fluid equilibrium, promoting the lattice equation of state to \(p = \rho\,c_s^2\,(1 + \theta)\). The density stays \(\rho = \sum_i f_i\) (read from streaming), so this is the regime for moderate density variation at \(\mathrm{Ma} < 0.1\) - natural convection rather than fire-scale plumes. For the large density ratios of strongly buoyant flows, use the variable-density models.low_mach closure below.

This is a parallel, opt-in path. The block requires LBM.thermal_model (the existing thermal switch), is supported on D3Q27 only, and is mutually exclusive with models.low_mach. With it unset the legacy theta = 0 behaviour stays bit-identical.

models:
  LBM:
    thermal_model: true   # activates the Tier 2 coupling; required by models.energy
  energy:
    T_ref: 1.0            # reference temperature: theta = T / T_ref - 1
    Pr: 0.7               # molecular Prandtl number: lambda = cp mu(T) / Pr
    cp: 1.0               # specific heat (caloric relation T = T_ref + h/cp)
    mu_ref: 0.05          # reference dynamic viscosity mu_0; omit -> nu_0 = (tau - 1/2)/3
    mu_exponent: 0.0      # power-law exponent n of mu(T) = mu_0 (T/T_ref)^n
    field:                # the SAME enthalpy / temperature field the low-Mach closure uses
      initial_field: "1.0"          # initial T(x, y, z) in lattice units
      wall_bcs:                     # temperature wall BCs (Dirichlet / Neumann / Robin)
        - pos: "W"                  # hot wall
          BC: TempDirichlet
          T_w: 2.0
        - pos: "E"                  # cold wall
          BC: TempDirichlet
          T_w: 1.0

The field block is the same EnergyFieldConfig the low-Mach closure uses (see the energy description below): an initial_field temperature expression, optional volumetric source_regions, and per-face temperature wall_bcs. There is one temperature producer shared by both tiers; the tier only sets how the fluid responds to it (a deviation in the equilibrium here, the EOS density for the low-Mach closure). When field is omitted, a uniform T_ref start with no heat source is assumed.

Buoyancy in this regime is supplied separately - either a uniform body force (LBM.F) or the Boussinesq force of a scalar buoyancy coupling. That force and the equilibrium \(\rho\theta\) term model distinct effects (the buoyant drive versus the thermal-expansion pressure) and are additive; do not represent the same buoyancy through both.

For the equilibrium / non-equilibrium derivation, the \(p = \rho c_s^2 (1+\theta)\) readout, and the validity regime, see the thermal modelling-hierarchy theory page.

Variable-density low-Mach closure

The optional models.low_mach block activates the variable-density low-Mach thermal closure of Taha et al. (2024), Tier 3 of the thermal modelling hierarchy. It is the regime for strongly buoyant, large-density-ratio flows (fire-scale plumes) that the Boussinesq and weakly-compressible routes cannot represent. The density is slaved to the ideal-gas equation of state \(\rho = P/(r\,T)\) rather than read from the population sum, and buoyancy is the exact \((\rho - \rho_\infty)\,\mathbf{g}\) term, of which the Boussinesq scalar force is the linearised small-\(\Delta T\) limit.

This is a parallel, opt-in path. With the block unset the solver keeps \(\rho = \sum_i f_i\) and the isothermal / Boussinesq route stays bit-identical. It is supported on D3Q27 only and is mutually exclusive with the legacy LBM.thermal_model flag.

models:
  low_mach:
    r: 1.0            # specific gas constant (lattice units): rho = P / (r T)
    P_thermo: 1.0     # uniform thermodynamic pressure; open domain keeps dP/dt = 0
    T_ref: 1.0        # reference / ambient temperature (lattice units)
    rho_inf: 1.0      # ambient density in (rho - rho_inf) * gravity; omit for P/(r T_ref)
    Pr: 0.7           # molecular Prandtl number: lambda = cp mu(T) / Pr
    cp: 1.0           # specific heat (caloric relation T = T_ref + h/cp)
    mu_ref: 0.05      # reference dynamic viscosity mu_0; omit -> nu_0 = (tau - 1/2)/3
    mu_exponent: 0.0  # power-law exponent n of mu(T) = mu_0 (T/T_ref)^n (0 = constant)
    gravity: [0.0, 0.0, -1.0e-5]   # gravity vector (lattice units), length = dimension
    domain_closure: open           # 'open' pins P to ambient, 'closed' tracks pressurisation
    energy:                        # finite-difference enthalpy / temperature field
      initial_field: "1.0"         # initial T(x, y, z) in lattice units (uniform T_ref here)
      source_regions:              # volumetric heat-release sources Q (a burner)
        - pos: "(x >= 6) & (x <= 10) & (y >= 6) & (y <= 10)"
          rate: 1.0e-3             # one heat-release rate Q baked in; all regions share it
      wall_bcs:                    # temperature wall BCs on the domain faces
        - pos: "W"                 # x-min wall held at a fixed temperature
          BC: TempDirichlet
          T_w: 2.0
        - pos: "E"                 # x-max convective (Robin) wall
          BC: TempRobin
          h: 0.1
          T_inf: 1.0

The energy block transports the specific enthalpy \(h\) each fluid step by an explicit finite-difference solver of the low-Mach energy equation (Taha et al. 2024),

\[ \rho\,\frac{\mathrm{D}h}{\mathrm{D}t} = \partial_\alpha\!\left(\lambda\,\partial_\alpha T\right) + Q, \qquad \lambda = \frac{c_p\,\mu(T)}{\mathrm{Pr}}, \]

and derives the temperature through the constant-\(c_p\) caloric relation \(T = T_\mathrm{ref} + h/c_p\) (datum \(h = 0\) at \(T_\mathrm{ref}\)), which closes the EOS density read by the fluid collision. Conduction uses a per-node effective diffusivity \(\alpha_\mathrm{eff} = \mu(T)/(\rho\,\mathrm{Pr})\) built from the local EOS density and the viscosity law \(\mu(T) = \mu_\mathrm{ref}\,(T/T_\mathrm{ref})^{n}\) (mu_ref, mu_exponent; mu_ref defaults to \(\nu_0 = (\tau - 1/2)/3\)), discretised in the conservative variable-\(\lambda\) face-flux form; advection \(u_\alpha\,\partial_\alpha h\) uses first-order upwinding, in lattice units. See the energy-equation theory page. The energy block is optional: when omitted, a uniform T_ref initial field with no heat source is assumed.

The source_regions.rate is a volumetric heat-release rate \(Q\) (energy per unit volume per unit time, lattice units): it enters \(\rho\,\mathrm{D}h/\mathrm{D}t\), so the update adds \(Q/\rho\). The kernel bakes in one rate, so every region must share the same rate.

Temperature wall boundary conditions are declared per domain face in wall_bcs. Each entry pins one cardinal face (W/E/S/N/B/F, the same convention as the fluid BCs: W = x-min, E = x-max, S = y-min, N = y-max, B = z-min, F = z-max) with one of three aliases:

• TempDirichlet (T_w): fixed wall temperature.

• TempNeumann (J_w, default 0.0): prescribed wall flux \(J_w = -\alpha\,\partial T/\partial m\) along the inward normal \(m\) (positive injects heat). The default 0.0 is an adiabatic wall.

• TempRobin (h, T_inf): convective wall \(J_w = h\,(T_\infty - T_w)\), with h a level-0 lattice velocity and T_inf the ambient temperature.

The wall sits at the half-cell (\(\Delta x/2\)) position of the halfway bounce-back fluid wall. At most one entry is allowed per face; a face with no entry keeps the adiabatic (zero-gradient) fallback, and a face on an axis declared periodic in models.BC.periodic_dims may not carry a wall BC.

Note

The thermal wall BCs are exact for constant viscosity (mu_exponent: 0, the default). With a temperature-dependent viscosity (mu_exponent != 0), the conservative face-flux stencil evaluates the wall conduction with the face-averaged diffusivity (node and extrapolated ghost), not the wall-temperature value, so the realized TempDirichlet wall conduction and the prescribed TempNeumann / TempRobin flux carry a first-order error in the one-cell conductivity variation at the wall. Refine at the wall if a variable-viscosity thermal BC needs higher fidelity.

The domain_closure selector sets how the spatially-uniform thermodynamic pressure \(P\) is closed, and it must match the fluid boundary conditions:

• open (default) pins \(P\) to ambient (\(\mathrm{d}P/\mathrm{d}t = 0\)). This is valid when the domain vents to the environment through an inlet or outlet. A domain sealed by wall BCs (every models.BC.BC_map entry a wall, with no inlet or outlet) that is also heated - a volumetric source_regions rate, or a non-adiabatic wall_bcs (TempDirichlet, TempRobin, or TempNeumann with \(J_w \neq 0\)) - is rejected under open, because holding \(P\) fixed would model the wrong thermodynamics (the trapped gas should pressurise as it heats). Use closed, or add an outlet.

• closed evolves \(P(t)\) from total-mass conservation for a sealed, heated enclosure. With \(\rho = P/(rT)\) and a fixed total mass \(M = \int \rho\,\mathrm{d}V\), the pressure follows \(P(t) = M r / \int (1/T)\,\mathrm{d}V\). The inverse-temperature integral is reduced over the leaf fluid nodes once per coarse step and the new \(P\) written back into the equation of state. closed assumes a genuinely sealed domain, so it is rejected if BC_map declares an inlet or outlet through which mass would leave.

Note

A fully periodic domain (no BC_map) is mass-conserving and is the homogeneous / slab-heating idiom: it stays on the open path by default, and may opt into closed to track the pressurisation.

The conserved energy variable (the energy density \(\rho h\)) is transported on a dedicated double-distribution lattice (the DDF energy backend), conservative by construction. Temperature is a derived diagnostic, never separately solved. See the double-distribution energy backend theory page.

For the governing equations, the thermodynamic / hydrodynamic pressure split, and the zeroth-moment reinterpretation, see the thermal modelling-hierarchy theory page.

Volumetric regions

A surface boundary condition acts on the domain faces; a volumetric region acts on every fluid node inside a region of the domain. Volumetric regions are configured under models.volumetric_regions and apply a porous-medium momentum sink to the fluid through the Guo body force, combining two contributions on region nodes:

• porous_alpha - a linear (Darcy) sink F[a] = -porous_alpha * u[a], damping momentum proportionally to the local velocity (e.g. an outlet sponge that absorbs pressure waves and near-outlet overshoots).

• porous_beta - a quadratic (Forchheimer / canopy) drag F[a] = -porous_beta * |u| * u[a], modelling pressure-decay / canopy drag whose resistance grows with the square of the speed (e.g. a forest canopy momentum sink).

simulations:
  - name: example
    models:
      volumetric_regions:
        # Linear Darcy outlet sponge (beta omitted -> 0).
        - pos: "(x >= 480) & (x <= 500)"
          porous_alpha: 0.2
          porous_beta: 0.05
        # Quadratic canopy drag region (shares the single alpha/beta).
        - pos: "(z >= 0) & (z <= 20) & (x >= 100) & (x <= 200)"
          porous_alpha: 0.2
          porous_beta: 0.05

Each region is placed by the same pos the surface BCs take, so a cuboid is (x >= xa) & (x <= xb) & ... as a predicate, a ground layer is B as a cardinal alias, and either can be written as a typed shape instead (see Placing a region with a typed shape). As on a surface BC, system writes the position in a non-lattice domain.systems frame. Both coefficients are >= 0 level-0 lattice values baked into the generated kernel and rescaled by 1/2**lvl per refinement level. All configured regions must share the same porous_alpha and the same porous_beta; porous_beta defaults to 0 (linear-only).

For the physics and the analytic validation anchors, see the porous-media theory page.

Boundary conditions

For CFD simulations, boundary conditions (BC) deserve a very very special topic. They can be the difference between a great result and a simulation diverging.

Here is how to configure the BCs for the simulations in Nassu

simulations:
  - name: example
    models:
      BC:
        # Periodic dimensions in the domain (x, y, z)
        periodic_dims: [false, false, false]
        # Wall model configurations
        WM_cfg:
          # Max error for TDMA
          TDMA_max_error: 5e-06
          # Max iterations for TDMA 
          TDMA_max_iters: 50
          # Number of divisions in TDMA
          # Must be an odd number
          TDMA_max_div: 25
        # Map of boundary conditions to use
        BC_map:
          # The boundary conditions are added from first to last
          # So if there is one BC that conflicts nodes and order with order,
          # The last one will be the one that stands
          
          - # Boundary condition name
            # Available ones are: 
            #  RegularizedHWBB
            #  RegularizedVelocityWall (alias RegularizedVelocityBB)
            #  UniformFlow
            #  Neumann
            #  RegularizedNeumannSlip
            #  RegularizedNeumannOutlet
            BC: RegularizedNeumannOutlet
            # Order in which apply the BC
            # 0 are the first applied, then 1, then 2 and so on
            order: 2
            # Position of the BC (x, y, z)
            # Options are: (N is the domain limit)
            # E: at (N, _, _)
            # W: at (0, _, _)
            # N: at (_, N, _)
            # S: at (_, 0, _)
            # F: at (_, _, N)
            # B: at (_, _, 0)
            # It's possible to combine these positions, as in NW: at (0, N, _)
            pos: E
            # Normal position to consider in BC
            # It points to outside the domain, so for a wall at the top (_, N, _)
            # you would use N
            # It's also possible to combine these directions, as in the position
            # Some BCs, such as UniformFlow, don't require a wall normal
            wall_normal: E
            # BC-specific parameters live under a nested `params:` field,
            # validated at parse time against the selected BC. The entry is
            # extra="forbid", so a misspelled key is rejected on load. BCs that
            # take no parameters omit `params:` entirely.
            params:
              rho: 1.0
          # Others BCs to apply
          - BC: RegularizedNeumannSlip
            order: 1
            pos: F
            wall_normal: F
          - BC: RegularizedHWBB
            order: 1
            pos: B
            wall_normal: B
          - BC: RegularizedNeumannSlip
            order: 0
            pos: N
            wall_normal: N
          - BC: RegularizedNeumannSlip
            order: 0
            pos: S
            wall_normal: S
          - BC: Neumann
            order: 0
            # It's also possible to combine the positions and normals to
            # combine the limitations of each name
            pos: NF
            wall_normal: N
          - BC: Neumann
            order: 0
            pos: SF
            wall_normal: S
        # Normalises density at a domain face after each step, preventing pressure drift.
        # List of normalization entries; empty by default (no normalization).
        rho_normalization:
          - pos: W   # Face to average rho on. Options: N, S, W, E, F, B
            rho_norm: 1.0  # Target density value. Defaults to 1

For the boundary conditions, check the documentation on each one for what they apply and how.

Predicate positions and coordinate systems

Besides a cardinal alias (N, S, … and their combinations), a BC pos accepts a boolean predicate over the node coordinates x, y, z, selecting every node that satisfies it (for example the interior baffle (x >= 14) & (x <= 17) & (8 <= y) & (y <= 24)). The grammar covers + - * / and integer powers, the comparisons < <= > >= == != and boolean and/or/not (&, |, ~); a non-cardinal predicate region also needs a wall_normal.

A predicate is interpreted in lattice coordinates by default. Set system to the name of a domain.systems frame to write the condition in that frame instead: each lattice node is mapped into the system and the predicate checked there, with the transform (translation, rotation and scaling) baked into the predicate at config load. This lets a region be placed in a physical or building-aligned frame rather than hand-converted to lattice units. system is only valid on predicate and shape positions, not on cardinal aliases. The same system field is available on models.volumetric_regions and on each scalar’s source_regions.

Placing a region with a typed shape

pos takes a kind-tagged mapping as readily as it takes a string, on a BC, a volumetric region and a scalar source region alike: a typed shape (point, points, line, plane, box, predicate) carrying the same field names an export probe uses, lowered to the equivalent node predicate at config load. The shape carries its own transformation (the entry’s system names the frame), and takes no dist (a spacing belongs to a probe, which generates sample points; a placement selects the nodes the lattice already has). See the boundary-condition guide for the field table and examples.

      BC:
        BC_map:
          - pos: {kind: plane, axis: z, axis_pos: 0.0, min: [10, 10], max: [50, 50]}
            BC: RegularizedHWBB
            wall_normal: B
            order: 1

Wall normals under a rotated system: a single-surface predicate (exactly one comparison, e.g. a half-space or a sphere) derives its wall normal analytically from the predicate, so a rotation is carried into the normal automatically and correctly. A composed region (a cuboid, any predicate with several comparisons) has no single analytic normal and falls back to the BC’s explicit cardinal wall_normal, which is interpreted in the lattice frame and is not rotated; under a rotating system that cardinal normal will not align with the rotated faces (the config load warns). Translate/scale-only systems preserve the axis directions, so a cardinal wall_normal stays valid. For a rotated wall whose BC needs an accurate normal, prefer a single-surface predicate, or supply a wall_normal already expressed in the lattice frame.

SEM

One special and important BC is the SEM (Synthetic eddy method) that applies a velocity profile with synthetic eddies to produce turbulence in the inlet. It enforces this condition on all nodes at the start of the domain, so at x=0. The height is obligatory in z direction and the “width” in y direction.

simulations:
  - name: example
    models:
      # Synthetic eddy method configuration
      SEM:
          # Specifications for SEM eddies generation
          eddies:
            # Lengthscale to use, may specify different values for each dimension
            lengthscale: {x: 14, y: 14, z: 14}
            # Volume density of eddies (n_eddies=eddies_vol_density*(SEM_volume/eddy_volume))
            eddies_vol_density: 15
            # Seed for random numbers to generate the eddies
            seed_rand: 0
            # All eddies are generated in this limit in yz domain (x is 2*lenghtscale.x)
            domain_limits_yz:
              start: [16, 0]
              end: [48, 96]
          profile:
            # Profile with velocity and R (rate od strain) for SEM
            csv_profile_data: "fixture/SEM/example/real_profile.csv"
            # Sums in z profile
            z_offset: 0
            # K multiplication factor
            K: 1
            # Multiplier for profile length.
            # Scales Z (z*height_mul) for csv_profile_data
            length_mul: 1
            # Multiplier for velocity values. 
            # Scales velocity values (ux*vel_height) for csv_profile_data
            # Scales Reynolds stress tensor (R*(height_mul^2)) for csv_profile_data
            vel_mul: 1

Note

For more informations on SEM, check its specific documentation.

Initialization

Three initialization modes are available: equation-based (default), file-based, and SEM field.

Equation-based initialization (default)

All macroscopic fields are specified as SymPy-compatible math expression strings in terms of the physical node coordinates x, y, z (lattice units at level 0). Constant values are a degenerate case (e.g. rho: "1").

Allowed functions: sin, cos, tan, asin, acos, atan, atan2, log, log2, log10, exp, sqrt, Abs, ceil, floor, Min/min, Max/max. Allowed constants: pi, e.

simulations:
  - name: example
    models:
      initialization:
        equations:
          # Constant density
          rho: "1"
          # ABL power-law profile (x is streamwise, z is height)
          ux: "0.05 * (z / 100.0) ** 0.25"
          uy: "0"
          uz: "0"

The scalar fields rho and u are also accepted and are automatically converted to equations:

      initialization:
        rho: 1
        u: { x: 0.05, y: 0, z: 0 }

File-based initialization

Load macroscopics from an .xdmf snapshot. The volume must contain rho and u fields covering at least the domain size. S is generated from finite differences of u if not present.

      initialization:
        macrs_filename: ./path/to/macrs.xdmf

Inlet field initialization

Prefills every fluid node with the mean velocity profile ux(z) of the active inlet turbulence method, broadcast homogeneously over the whole domain (instead of a constant). Requires an inlet turbulence method - models.BC.SEM. Eliminates the long spin-up period needed when starting from a uniform field.

uy, uz are zero-initialized. The rate-of-strain S is computed from finite differences of the ux(z) profile. theta and omega_LES are zeroed. rho is taken from equations.rho (must be a constant, default "1").

      initialization:
        inlet_field: true
        equations:
          rho: "1"

The name sem_field is accepted as an alias.

Important

For checkpoint details, check the data section.

Devices and precision to use

Nassu is configured to run in GPU devices, using NVIDIA CUDA library. Some options are available for the setup of the devices.

It’s also required to specify the floating point precision to use for the simulation.

simulations:
  - name: example
    models:
      # Engine options
      engine:
        # Devices numbers to use
        # Defaults to none, devices are chosen automatically in ascending order
        devices_numbers: [0]
        # Number of devices to use
        # Defaults to one (no support for multiple devices)
        n_devices: 1
        # Engine to use
        # Options are: CUDA
        name: CUDA
      # Numerical precision options
      precision:
        # The default precision, must be specified
        # Options are: single | double
        default: single
        # All options below are: single | double | default
        # Defaults to: default

        # Precision to perform calculations
        calculations: default
        # Precision to store macroscopics fields
        macroscopics: default
        # Precision to use for populations in shared memory
        populations: default

Note

The precision used for the populations impacts the block size that can be used, because it’s allocated in the GPU shared memory, which is very limited.

Important

Multiple devices are not supported.