Procedures

Incremental running-mean update: $\bar{x}_n = \bar{x}_{n-1} + (x_n - \bar{x}_{n-1})/n$ where stat_inc = 1/n.

Add the configured explicit SGS closure to the momentum RHS.

Compute derived fields from the velocity gradient tensor.

Interpolates the pressure field from CELL (DIR_Z) to VERT (DIR_X) for snapshot output and rescales from pseudo-pressure to physical (kinematic) pressure.

From the global rank maps, fills in the rank position as well as the previous and next rank in the par structure

Copies cell_dims information to vert_dims taking periodicity into account

Copies vert_dims information to cell_dims taking periodicity into account

Allocate memory for a new block and return a pointer to a new m_allocator object.

Constructor function for the cuda_tdsops_t class. See tdsops_t for details.

Wrapper for curl

Curl of a vector field (u, v, w).

Performs 2D mesh decomposition using 2decomp&fft

Full (forward + backward) non-periodic pentadiagonal Thomas solve.

Full (forward + backward) non-periodic pentadiagonal Thomas solve on CPU. Mirrors the CUDA der_penta_full kernel; operates on 2D (SZ, n) arrays.

Periodic pentadiagonal Thomas solve via Sherman-Morrison-Woodbury rank-4.

Cyclic pentadiagonal compact solve via Sherman-Morrison-Woodbury rank-4 correction. A_cyc = A_np + U*W^T where U = [e_1,e_2,e_{n-1},e_n] and W^T encodes the four corner coupling pairs; A_np is the non-periodic Thomas factorisation in (ffr,...).

Go through the block list from head to tail, deallocating each memory block in turn. Deallocation of a m_allocator object automatically deallocates its internal allocatable data array.

Divergence of a vector field (u, v, w).

Wrapper for divergence_v2p

Gathers the non-periodic x line into a periodic one so that a plain FFT can stand in for the cosine transform. Unlike the y version below this handles an odd nx, which the 100 case requires.

Combined X and Y periodicity enforcement (interleave shuffle).

Gathers the non-periodic x and, when it is not split across ranks, y lines into periodic ones so that a plain FFT can stand in for the cosine transform in both directions.

Fills c_mirror_dev and c_plane_recv_dev with the dim2 mirror of the local spectral slab, so that the 100 paired split can read its partner mode without going looking for it in another rank's memory.

Single-GPU non-periodic pentadiagonal compact-FD solve.

Single-pass non-periodic pentadiagonal compact-FD solve on CPU. No MPI exchange is needed; each group is solved independently.

Cyclic pentadiagonal compact-FD solve via Sherman-Morrison-Woodbury. Halos u_recv_s/e must be set to periodic-extension values by the caller.

Cyclic pentadiagonal compact-FD solve via Sherman-Morrison-Woodbury. Halos u_recv_s/e must be set to periodic-extension values by the caller.

Pack the first dims Cartesian points of field f into a contiguous device buffer, in (i, j, k) order without padding.

DIR_X fields are reordered and unpadded by a single kernel. DIR_C fields are unpadded directly, and DIR_Y/DIR_Z fields are reordered to DIR_C first. Precision conversion happens in the same kernel.

Backward FFT for the non-periodic x case, undoing the x-y transpose applied by fft_forward_100_omp. Only the unpadded extents are written, which is all undo_periodicity_x_omp goes on to read.

Backward FFT for 110 case: C2R then transpose (nz,nx,ny)->(nx,ny,nz)

Backward FFT for the non-periodic x and y case, undoing the transpose applied by fft_forward_110_omp. Only the unpadded extents are written, which is all undo_periodicity_xy_omp goes on to read.

Forward FFT for non-periodic-X case We transpose X<->Y so that the non-periodic direction becomes the "Y" direction in the transposed space, then use the same FFT approach as 010

Forward FFT for the non-periodic x case. Transposes x and y so that the non-periodic direction sits in dim2, where it keeps its full mode range, then transforms the (ny, nx, nz) buffer with the swapped plan set up in init.

Forward FFT for 110 case: transpose (nx,ny,nz)->(nz,nx,ny) then R2C

Forward FFT for the non-periodic x and y case. Transposes the physical field (nx, ny, nz) to (nz, nx, ny) so that the r2c transform runs along the periodic z direction, then transforms with the plan set up in init.

if stretching in y is 'bottom'

Post-process div U* in spectral space for the non-periodic x case.

Spectral post-processing for 110 case (R2C Z-transpose) 7 separate kernel launches to avoid cross-block race conditions. Spectral array is (nz/2+1, nx, ny) = (nx_spec, ny_spec, nz_spec). Thread i->X (ny_spec), blockIdx%y k->Y (nz_spec), serial j->Z (nx_spec).

Post-process div U* in spectral space for the non-periodic x and y case. After the (nx, ny, nz) -> (nz, nx, ny) transpose applied by fft_forward_110_omp, dim1 holds the periodic z r2c modes, dim2 holds the non-periodic x modes and dim3 holds the non-periodic y modes.

m_base_backend

m_base_backend

Set a face of f using values supplied by another field f_start. Both fields must be DIR_X VERT pencil-layout. X_FACE uses the inlet plane (pencil index 1) of f_start plus a convective outflow update at the right face. Y_FACE copies the bottom and top y-pencil planes of f_start onto the corresponding planes of f (Dirichlet only).

m_base_backend

Set domain X_FACE boundary values. c_start: Dirichlet value applied at the left face (i = 1) c_end: convective velocity Uc = uxmax * gdt / dx, used as multiplier in the outflow scheme du/dt + Uc*du/dx = 0

Set domain X_FACE boundary values with a spatially-varying inlet.

Set domain Y_FACE boundary values. c_start: Dirichlet value applied at the bottom face (j = 1) c_end: Dirichlet value applied at the top face (j = ny)

Set domain Y_FACE boundary values from another field. Both f and f_start are DIR_X VERT pencil-layout. Only the bottom (j = 1) and top (j = ny) y-pencil planes of f_start are read.

volume integral of a field

volume integral of a field

Volume-equivalent grid-filter width, Delta=(dxdydz)^(1/3).

Release resources owned by the solver and its runtime models.

Return a pointer to the first available memory block, i.e. the current head of the block list. If the list is empty, allocate a new block with create_block first.

Utility function that returns an array made of the id of the block currently in the block list. Return the array [0] if block list is empty.

Get the coordinates of a vertex with i, j, k local cartesian indices Avoid calling this in hot loops

Getter for local domain dimensions

Getter for domain dimensions

Extract data from field f optionally reordering into dir orientation. To output in same orientation as f, use call ...%get_field_data(data, f, f%dir)

Resolve a field name to the corresponding solver field pointer. Centralises the name-to-field mapping used by I/O managers.

Getter for local domain dimensions

Get application storage directional index from cartesian index

Get cartesian index from application storage directional one

Converts indices in between any two DIR_?

Getter for the main dimension a field oriented along dir with data on data_loc

Getter for the main dimension of field phi

Returns RDR_?2? value based on two direction inputs

Gradient of a scalar field 'p'.

Wrapper for gradient_p2v

Check whether a field name is present in the output_fields list.

if stretching in y is 'bottom'

Interpolation of a scalar field from cell centres to vertices.

Returns whether or not the current rank is the root rank

Laplacian of a scalar field 'u'.

Copy data between x3d2 padded arrays and cuFFTMp descriptors

Copy with transpose: src(nx, ny, nz) -> dst(ny, nx, nz) Used for 100 case forward FFT

Copy with transpose back: src(ny, nx, nz) -> dst(nx, ny, nz) Used for 100 case backward FFT

Completely initialise the mesh object. Upon initialisation the mesh object can be read-only and shouldn't be edited Takes as argument global information about the mesh like its length, number of cells and decomposition in each direction

The only rank is the root, so its buffer is already the broadcast value.

Reported as initialised whether or not MPI_Init has been called. There is no MPI environment to bring up here, so callers asking whether it is safe to proceed always have their answer, and a library handed a communicator never has to wait for one.

Reachable only where a call site has already branched on a rank count above one, as the exchanges in the CUDA Poisson solver do. It exists so that callers can name it in a use m_mpi, only: list that has to compile in both builds.

Elapsed wall-clock seconds, from an arbitrary origin as MPI defines it.

Whether ranges are emitted: X3D2_NVTX (default on), read on first use.

Obtains global coordinates for all the vertices and midpoints

Copy the unpadded (nx, ny, nz) part of a DIR_C field

Single precision output variant of pack_c_dp

Copies the first dim2 plane of the local spectral slab into a contiguous buffer so it can be handed to MPI.

Reorder a DIR_X field into an unpadded (nx, ny, nz) Cartesian buffer in one pass: reorder_x2c followed by dropping the padding.

Single precision output variant of pack_x2c_dp

Backward counterpart of pair_fw, symmetric in the same way.

Forward paired split of one dim2 mode against its mirror. The rule is symmetric, so the partner entry is pair_fw(partner, own, i2_rev).

Divide one entry by its spectral equivalence constant, with the eps_wave guard that keeps roundoff noise out of the pressure field.

Backward counterpart of postprocess_dim13_fw.

Forward postprocess of one spectral entry in dim3 then dim1.

LU preprocessing for non-periodic pentadiagonal Thomas algorithm.

Post-processes the divergence of velocity in spectral space, including scaling w.r.t. grid size.

Post-process div U* in spectral space for all periodic BCs.

Post-processes the divergence of velocity in spectral space, including scaling w.r.t. grid size.

Post-process div U* in spectral space, for non-periodic BC in y-dir.

Post-processes the divergence of velocity in spectral space, including scaling w.r.t. grid size.

Post-processes the divergence of velocity in spectral space, including scaling w.r.t. grid size.

Solve the Poisson equation at cell centres with non-periodic BC along y

Post-process div U* in spectral space, for non-periodic BC in x-dir.

Stage 1 (forward): normalisation followed by the dim3 and dim1 postprocess. Every mode is handled on its own here.

Stage 3 (backward): the dim2 paired recombination, then the dim3 and dim1 postprocess.

Stage 2 (forward): the dim2 paired split, then the Poisson solve.

Post-process div U* in spectral space, for non-periodic BC in both x- and y-directions (R2C Z-transpose layout: dim1 = Z r2c modes, dim2 = X modes, dim3 = Y modes).

Backward third of the spectral post-process for the 110 case: the X paired even/odd recombine, then the Z periodic post-process backward.

Forward third of the spectral post-process for the 110 case (non-periodic x and y, periodic z, R2C Z-transpose layout): normalise, Z periodic post-process forward, then the X paired even/odd split.

Step 1 (forward): normalise + Z periodic post-process Z is dim1 (serial j loop), periodic R2C — no sign flip needed since j only goes to nz/2+1.

Step 4: Poisson solve — divide by waves

Middle third of the spectral post-process for the 110 case: the Y paired even/odd split, the Poisson divide, and the Y paired even/odd recombine.

Step 6 (backward): X paired even/odd recombine

Step 2 (forward): X paired even/odd split X is dim2 (thread i). Only i in [2, nx/2+1] executes. Writes to i and nx-i+2. Race-free (pair doesn't enter).

Step 5 (backward): Y paired even/odd recombine

Step 3 (forward): Y paired even/odd split Y is dim3 (blockIdx%y = k). Only k in [2, ny/2+1] executes. Writes to k and ny-k+2. Race-free (pair block doesn't enter).

Step 7 (backward): Z periodic undo

Release memory block pointed to by HANDLE to the block list. It is pushed to the front of the block list, in other words it is made the head block.

Runs the solver forwards in time from t=t_0 to t=T, performing postprocessing/IO and reporting diagnostics.

m_base_backend

m_base_backend

m_base_backend

Reduces a single slice f(i_slice, :, :) (in the kernel's packed-pencil indexing). Signed max/sum, no abs. One thread per (y, z) point.

m_base_backend

|S|=sqrt(2 S_ij S_ij), where S_ij=(du_i/dx_j+du_j/dx_i)/2.

Stretching necessitates a special operation in spectral space. The coefficients for the operation are stored in matrix form.

Whether export_field_to_device can pack field f without staging it through host memory. Backends without device memory keep this default and report no support.

Constructor function for the tdsops_t class.

Thomas algorithm implementation. So much more easier than the distributed algorithm. It is intended to work only on a single rank so there is no MPI communication.

Skew-symmetric form of convection-diffusion terms in the incompressible Navier-Stokes momentum equations, excluding pressure terms. Inputs from velocity grid and outputs to velocity grid.

low memory version of the transport equation, roughly %2 slower overall

Skew-symmetric form of convection-diffusion terms in the species equation. Inputs from velocity grid and outputs to velocity grid.

Compute the convection and diffusion for the given field in the given direction. Halo exchange for the given field is necessary When sync is true, halo exchange of momentum is necessary

Compute the convection and diffusion for the given field in the given direction. Halo exchange for the given field is necessary When sync is true, halo exchange of momentum is necessary

Used before R2C FFT for 110 case to put Z (periodic) in fast dim.

Transpose back: src(nz, nx, ny) -> dst(nx, ny, nz) Used after C2R FFT for 110 case.

Scatters the gathered x line back to its original ordering.

Combined X and Y periodicity undo (reverse interleave shuffle).

Scatters the gathered x and, when it is not split across ranks, y lines back to their original ordering.

m_base_backend

m_base_backend

Global sum of a2 + b2 + c**2, with one MPI reduction.

Global sum of a2 + b2 + c**2, with one MPI reduction.

Mason-Thomson blend used by Incompact3d's ABL Smagorinsky path.

Spectral equivalence constants

Computes and reports monitoring quantities, writing to both stdout and the monitoring file. Enstrophy: $\mathcal{E} = \frac{1}{2N} \sum |\nabla \times \mathbf{u}|^2$ Divergence: $\nabla \cdot \mathbf{u}$ max and mean (should be ~0)