Procedures

Report a fatal error and stop every rank of the run.

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

Adds a uniform constant c to a single x-normal plane (pencil index i_plane) of a DIR_X VERT field, in place.

Allocates the backend's own tdsops type and fills in the coefficients of the requested operation.

Pre-correction (called per substep after the integrator step): enforce the inflow Dirichlet and convective outflow BCs, then apply an exact uniform mass correction to the outlet, following Incompact3d's pre_correc.

Add the configured explicit SGS closure to the momentum RHS.

Apply the explicit low-pass filter to the velocity, as Incompact3d does for wall-modelled ABL runs (ifilter=1, C_filter).

Copies the first and last four points of every line into the start and end send buffers. n is the line length, n_u the padded extent.

Abort the run if a CUDA runtime call did not return cudaSuccess.

Abort the run if a cuFFT call did not return CUFFT_SUCCESS (0). cuFFT has no error-string function, so only the code is reported.

Checks the fields passed to the labelled field operations share a valid data location and the DIR_X layout.

Checks the fields passed to the labelled field operations share a valid data location and the DIR_X layout.

Abort the run if an MPI call did not succeed.

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.

Computes the Q-criterion from velocity gradients

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

Forms nut(scale_agradient_a + scale_b*gradient_b).

Computes nut=l_s^2sqrt(2S_ij*S_ij) from velocity gradients.

Computes the vorticity magnitude from velocity gradients

Copies cell_dims information to vert_dims taking periodicity into account

Copy a regular 3D array in host memory into the specialist data structure field that lives on device or host

Copy the specialist data structure from device or host back to a regular 3D data array in host memory.

Periodic wrap-around on a single rank: copies the first n_data elements of a halo buffer, the groups of the direction being exchanged. The buffers are sized for the direction with the most groups, so a whole-array assignment would copy that maximum on every exchange (16 times the data in use on an elongated mesh).

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

Boundary Conditions hook (called per substep before transeq). Builds the inlet Dirichlet profile only; outflow parameters are computed in apply_BC from the post-integration velocity. u = bc_start_u + inlet_noise(1) * um * (2r - 1) v = bc_start_v + inlet_noise(2) * um * (2r - 1) w = bc_start_w + inlet_noise(3) * um * (2r - 1)

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

Full non-periodic pentadiagonal solve on a single subdomain.

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

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,...).

Cyclic pentadiagonal solve via Sherman-Morrison-Woodbury (rank 4): A_cyc = A_np + UC with U = [e_1|e_2|e_{n-1}|e_n]. z holds A_np^{-1}U (see omptgt_tdsops_t%penta_z_dev); each line forms M = I + CZ, solves Mc = Cy and applies y -= Zc.

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

First phase of the distributed solve: builds the RHS and runs the forward and backward hybrid sweeps, leaving the 2x2 boundary values in the send buffers.

Second phase of the distributed solve: resolves the 2x2 boundary systems with the neighbours' values and substitutes back.

Non-periodic Thomas solve along each line of a single subdomain.

Periodic Thomas solve along each line of a single subdomain.

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.

Allocates a device array of the given extents.

Allocates a device array of a host array's shape and copies it in.

Copies a device array into a host array of the same size.

Frees a device array allocated with dev_alloc.

Copies a host array into a device array of the same size.

Waits until the device work queued so far has completed.

Unnormalised smearing kernel of disc at coords: Gaussian in the rotor-normal direction (delta_N), super-Gaussian (^8) in the radial direction (delta_R) to give a flat-topped disc.

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 paired split can read its partner mode without going looking for it in another rank's memory.

Single-subdomain non-periodic pentadiagonal compact solve.

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 solve via Sherman-Morrison-Woodbury. Halos u_recv_s/e must hold the periodic extension.

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). 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). On a single rank X and Y are both whole, so 7 kernel launches avoid cross-block race conditions. On multiple ranks cuFFTMp splits X across ranks, so the X pairing reaches its partner through the dim2 mirror exchange (see exchange_dim2_mirror) instead, adding two more launches.

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.

Adds a uniform constant c to one x-normal face plane of a DIR_X field, in place. Used to enforce exact outlet mass conservation: after the convective outflow update, the outlet-plane mean is shifted uniformly so that it matches the prescribed inlet mean (mirrors Incompact3d's bxxn = bxxn - ut + ut1). at_end = .true. -> the x = nx plane (outlet) at_end = .false. -> the x = 1 plane (inlet)

field_add_const_x_face_base

field_add_const_x_face_base

Add the face planes of g onto those of f, on Dirichlet faces only. Both fields are DIR_X; X_FACE and Y_FACE are supported.

m_base_backend

m_base_backend

f += g on the x faces flagged by set_start/set_end.

f += g on the x-faces (pencil index 1 and nx) flagged by set_start/set_end. Launch as field_set_x_face_from_field.

f += g on the y faces flagged by set_start/set_end.

f += g on the y-faces (j = 1 and ny) flagged by set_start/set_end. Launch as field_set_y_face.

Adds coeffs(1, l)g, coeffs(2, l)g and coeffs(3, l)*g to a, b and c at every point of the DIR_X fields labelled l in 1..size(coeffs, 2). Points with any other label are left unchanged.

field_labelled_axpy_base

field_labelled_axpy_base

Computes this rank's sums of ga, gb and g*c separately for each label in a single pass over DIR_X fields. label holds integer values stored as reals; points labelled l in 1..size(sums, 2) add to sums(:, l) and all other points are skipped. No MPI reduction is done, so the caller can reduce all labels together in one MPI call.

field_labelled_sums_base

field_labelled_sums_base

Obtains maximum and mean values in a field

field_max_mean_base

field_max_mean_base

Maximum and sum of |f|.

First pass of field_plane_sums: partial(j, z) is the sum over x of y index j at z, summed in a fixed order.

First pass of field_plane_sums: partial(j, z) is the sum over x of y index j at z. Launch with threads = dim3(SZ, 1, 1) and blocks = dim3(n_y_blocks, nz, 1).

Sum a DIR_X field over x and z on every local y-plane: sums(j) is the sum over the plane of y index j. The order of summation is fixed, so the result is reproducible.

m_base_backend

m_base_backend

Second pass of field_plane_sums: sums(j) is the sum over z of partial(j, z), in a fixed order.

Second pass of field_plane_sums: sums(j) = sum over z of partial(j, z), in a fixed order. Launch with one thread per y index.

Scales a field by a

Drag law tau = drag_coeffu_c|u_h| from each column's own sampled velocity, written into one stress plane (component 1: x, 3: z).

Drag law tau = drag_coeffu_c|u_h| from each column's own sampled velocity, written into one stress plane (component 1: x, 3: z).

Replace one SGS stress plane using the local sampled velocity.

m_base_backend

m_base_backend

A field is a subdomain with a rectangular cuboid shape. It has 6 faces, and these faces are either a subdomain boundary or a global domain boundary based on the location of the subdomain. This subroutine allows us to set any of these faces to a value, 'c_start' and 'c_end' for faces at opposite sides. 'face' is one of X_FACE, Y_FACE, Z_FACE from common.f90 Optionally, bc_start/bc_end select the BC type (default BC_DIRICHLET).

As field_set_face but with a spatially-varying inlet face field instead of a scalar c_start.

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).

field_set_face_base

Sets the x faces: a Dirichlet value c_start at the inlet (pencil index 1) and the convective outflow du/dt + Uc*du/dx = 0, Uc = c_end, at the outlet (pencil index nx). BC_NEUMANN (1) leaves a face untouched.

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

As field_set_x_face, with the inlet taken from the inlet plane of f_start rather than a constant.

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

Overwrites the x = 1 plane (pencil index 1) of a DIR_X field with the host array plane(ny, nz), in local VERT dims. Only the plane crosses to the device, so a time-varying inlet can be refreshed without uploading and reordering a full 3D field.

field_set_x_inlet_plane_base

field_set_x_inlet_plane_base

Sets the bottom (y = 1) and top (y = ny) faces where flagged, so the two y faces can differ (ABL: no-slip floor, free-slip lid).

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_start/set_end: whether that face carries a prescribed value at all, so the two y-faces can differ (ABL: no-slip floor, free-slip lid).

Copies the bottom and top y planes of f_start onto those of f.

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.

Sets the y plane plane to c.

Set one interior y-plane to a constant. The caller resolves the plane into its position within a y-block and the group offset for that block, so this kernel stays free of layout arithmetic.

Set one interior y-plane of a DIR_X field to a constant.

m_base_backend

m_base_backend

Shifts a field by a

Reduces field to a scalar by integrating over the volume

volume integral of a field

volume integral of a field

Fill host_field with the right-hand side of the Poisson equation. The routine uses the host-associated mesh, dims and n_pi. Edit this routine to supply another right-hand side, then rebuild. The coordinates are the cell centres from mesh%get_coordinates(i, j, k, CELL). The right-hand side must have zero mean. In every supported combination of boundary conditions, the pressure has periodic or homogeneous Neumann conditions. The default is cos(2 pi x).

Filter one velocity component in all three directions, in place.

Filter a DIR_X field along one direction, in place.

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

Release resources owned by the solver and its runtime models.

Friction velocity from the magnitude of a mean wall-stress vector.

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

Address of the block's storage in a separate device allocation, null for a host block. A device field type that keeps its storage behind a raw device pointer overrides this, letting its backend reach the storage by dynamic dispatch rather than by a type query.

Reports the memory clock rate and bus width of the target device, used to work out a theoretical peak bandwidth.

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'

Allocates and initialises the backend's FFT based Poisson solver.

Interpolation of a scalar field from cell centres to vertices.

Returns whether or not the current rank is the root rank

a, b, c += coeffs(1:3, l)*g at the points labelled l.

a, b, c += coeffs(1:3, l)*g at the points labelled l; launch as vector_norm_squared.

s(1:3, l) += sums of ga, gb and g*c over the points labelled l, as cuda/kernels/fieldops.f90: each thread accumulates runs of equal labels along its pencil and flushes a run when the label changes, so the atomics scale with the label changes rather than the points. Atomics instead of reduction(+:s), which flang cannot do on the GPU for an array whose shape is only known at run time.

s(1:3, l) += sums of ga, gb and g*c over the points labelled l; launch as vector_norm_squared. Each thread accumulates runs of equal labels along its pencil and flushes a run when the label changes, so the atomics scale with the label changes rather than the points.

Laplacian of a scalar field 'u'.

Detail for developers, hidden unless the level is LOG_LEVEL_DEBUG.

An error the caller handles itself. Use abort_run to stop the run.

Progress and setup information.

Something is suspicious, but the run can carry on.

Tridiagonal low-pass filter, matching Incompact3d's filter.

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

Copies n reals between two device addresses.

Copies n reals from device address src to host address dst.

Copies n reals from host address src to device address dst.

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

Neutral rough-wall drag coefficient, (kappa/ln(h/z0))2, so that the wall stress is drag_coeffu|u| at the sampling height h.

Neutral rough-wall drag law. Returns the wall value of the SGS stress tau_xy (and tau_yz), signed so that a positive sample gives a positive stress and hence a momentum sink once differentiated.

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

Obtains global coordinates for all the vertices and midpoints

Initialises tdsops in place. See tdsops_t for the arguments.

Packs the physical Z-slab field block (nx, ny, nz_l) into the MPI_Alltoall send buffer, selecting for each destination rank s the Y range [sny_l+1, (s+1)ny_l] it owns in the Y-slab layout.

Inverse of unpack_110_yslab: packs the Y-slab real view into the MPI_Alltoall send buffer, selecting for each destination rank s the Z range it originally owned in the physical Z-slab layout.

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).

Read nx ny nz, then the optional bc_x bc_y bc_z, then the optional --repeat N. The bc values are periodic or dirichlet. The default is periodic. The default for --repeat is 1. Any other input prints one usage line and stops.

Solves A_np*z_k = e_{p_k} for p = 1, 2, n-1, n with the non-periodic pentadiagonal LU factors, as thread 1 of der_penta_periodic does in the CUDA backend.

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.

Incompact3d's pre_correc: on each Dirichlet face, add the previous projection's tangential correction to the value set by the case. The projection that follows then leaves the face at its prescribed value plus dt(grad p_old - grad p_new), instead of -dtgrad p_new. The kept correction is dt_last*grad p_old, so it is rescaled by dt_now/dt_last first, as Incompact3d does with gdt(itr). The timestep is time_intg_t%stage_dt, the time the substep has advanced from the divergence-free start of the step, which is what the projection scales its correction by. For AB this ratio is always 1; for RK it changes between stages (1/2, 3/4, 1 for RK3). Before the first projection (or after a restart, as in Incompact3d) there is no correction to add.

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. nx, ny are the local loop bounds (equal to nx_glob, ny_glob on a single rank); the normalisation divisor always needs the globals.

Step 4: Poisson solve — divide by waves nx is the local loop bound; the Nyquist test always needs nx_glob.

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

Multi-rank counterpart of process_spectral_110_x_pair_bw, see process_spectral_110_x_pair_fw_mirror.

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).

Multi-rank counterpart of process_spectral_110_x_pair_fw: the X paired split, but the partner mode lives on another rank, reached through the dim2 mirror exchange (see exchange_dim2_mirror). Unlike the single-rank kernel this only updates this rank's own modes; the partner rank updates its own copy the same way. At the global Nyquist mode (even grid size) the mirror partner is this rank's own pre-update value, so the formula below self-pairs it exactly as the single-rank kernel does.

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

Memory clock (kHz) and bus width (bits) of OpenMP device device. The vendor this build targets is tried first, then the other one, as the build only guesses the vendor for some compilers.

Read turbine disc parameters from a .ad file. Format: one header line followed by one row per disc of CoR(x) CoR(y) CoR(z) Yaw[deg] Tilt[deg] RotorDiam C_T alpha

Redistributes the 110 field from the physical Z-slab (nx, ny, nz_l) to the cuFFTMp Y-slab real view (2*(nz/2+1), nx, ny_l), ahead of the forward transform. See the comment above pack_110_xyz_to_yslab.

Inverse of redistribute_110_to_yslab, after the backward transform.

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.

Reorders u_i in layout dir_from into u_o in layout dir_to.

reorder subroutines are straightforward, they rearrange data into our specialist data structure so that regardless of the direction tridiagonal systems are solved efficiently and fast.

Views the device storage of u as an array of its current shape. The pointer holds a device address: pass it to a kernel, never index it on the host.

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

Calculates the scalar product of two input fields

Scalar product of two DIR_C fields, whose storage is the padded Cartesian box n_c with the physical points in 1:dims.

scalar_product_base

scalar_product_base

Deterministic seeding of the intrinsic RNG, distinct for every (rank, start iteration) pair. A fresh run starts at iter = 0.

Drop messages below level (one of the LOG_LEVEL_* parameters).

Copies plane(y, z) into the x = 1 plane (pencil index 1) of a DIR_X VERT field, where f(i, 1, b) holds y = (y_block - 1)SZ + i and z = b - (y_block - 1)nz. Launch as add_scalar_x_face.

Signed maximum and sum of the slice f(:, i_slice, :).

Reduces a single slice of f at index i_slice along f's DIR axis. Returns signed max (not abs) and signed sum. No division by count. Caller is responsible for MPI_Allreduce across ranks.

slice_max_sum_base

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.

slice_max_sum_base

|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.

u += u_d, where u is DIR_X and u_d the same field in layout dir_from.

sum_yintox combines a y directional field into the corresponding x directional field.

As sum_yintox, for a z directional field.

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.

Waits for outstanding device work.

Applies a tridiagonal operator to u, the exact algorithm used being decided at runtime. Backend implementations are responsible from directing calls to tds_solve to the correct algorithm.

Constructor function for the tdsops_t class.

As tds_solve, solving the tridiagonal system with the Thomas algorithm.

First phase for du = d(u), dud = d(u*v) and d2u = d2(u) in one sweep.

Second phase: substitutes back and combines r_du = -1/2(convdu + dud) + nu*d2u.

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.

As transeq_x, for a scalar species transported by uvw.

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

transeq equation obtains the derivatives direction by direction, and the exact algorithm used to obtain these derivatives are decided at runtime. Backend implementations are responsible from directing calls to transeq_x into the correct algorithm.

As transeq_x, for the y direction.

As transeq_x, for the z direction.

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.

Unpacks the MPI_Alltoall receive buffer into the padded real view (2(nz/2+1), nx, ny_l) of the cuFFTMp Y-slab descriptor: data received from source rank s reconstructs the global Z range [snz_l+1, (s+1)*nz_l].

Inverse of pack_110_xyz_to_yslab: unpacks into the physical Z-slab field block (nx, ny, nz_l). dst must be zeroed by the caller first, as fft_backward_110_cuda already does for the single rank path.

Validate relationships between ABL namelist values in one place.

adds two vectors together: y = ax + by

copy vectors: y = x

pointwise multiplication between two vectors: y(:) = y(:) * x(:)

vecmult_base

vecmult_base

Computes the global discrete squared norm of a three-component field.

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.

Warn on rank 0 if the right-hand side does not have zero mean. The solver needs a zero-mean right-hand side for these boundary conditions. It discards the mean component. The solve runs in both cases.

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)

Groups between successive z in the OMP DIR_X ordering, group = stride*(z - 1) + y_block. It is the padded y-block count (see get_index_dir in m_ordering), so it can exceed the number of y-blocks a field uses, e.g. a cell field on a non-periodic y.