Procedures

Where the ABL wall model samples the velocity, dsampling grid spacings above the floor (vertex 1). As in Incompact3d's wall_sgs_noslip, the sample is interpolated linearly between the vertices below and above that height, so a non-integer dsampling is honoured rather than rounded to a grid plane.

The y-vertex above the ABL wall-model sample plane, checked against the interpolation weight. With a zero weight the upper plane is not read, so the sample may then sit on the last vertex.

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.

a, b, c += src(4:6, s)*exp(-dist^2/eps^2) within r_cut of the point src(1:3, s), over the global index box box(1:3, s)..box(4:6, s); see source_index_boxes. Box indices are wrapped across periodic boundaries and skipped where this rank does not own them. Launch with blocks (n_src, box extent in z) and threads striding y within the box; sources overlap, so the adds are atomic.

a, b, c += src(4:6, s)*exp(-dist^2/eps^2) within r_cut of the point src(1:3, s), over the global index box box(1:3, s)..box(4:6, s); see source_index_boxes. Box indices are wrapped across periodic boundaries and skipped where this rank does not own them. One thread per (source, z, y) of the box, each sweeping x; sources overlap, so the adds are atomic. grid holds the local VERT dims, global VERT dims, rank offsets and periodic flags (1/0), three entries each, then the padded y-block count.

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

Advance the rotor by one time step: update the angular velocity (and pitch) from the operating mode, then rotate the blades. For the controlled modes u_upstream must already hold the probe velocity.

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

Build the flow case with the given name.

(Re)allocate the station and element arrays of line.

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

Sets up the state every case shares. The backend's init calls it once it has set bc_code, which the spectral extents it passes in already depend on. The supported BC codes are the cases of the backend factories, create_omp_poisson_fft and create_cuda_poisson_fft.

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.

True when another process changes its card memory by more than 128 MiB or an NVML sample fails. Log the skip line once when emit_line is true

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.

Abort the run if an MPI call did not succeed.

Checks the fields passed to the point sampling and source operations are DIR_X VERT fields.

The FITS, BORDERLINE or DOES_NOT_FIT verdict of every estimate.

The BC parse and four-way split of classify_boundaries in m_mesh.

Parse the BC strings of the three directions. A direction is periodic only if both ends are 'periodic'. bc_code is the FFT_BC_* class of the resulting periodicity, or FFT_BC_UNSUPPORTED. Shared by mesh_init and callers with only BC strings (no mesh_t, e.g. static memory estimation).

Compute derived fields from the velocity gradient tensor.

Blade-element loads of every element from the sampled fluid velocity, the body velocity and the static polars.

Convective outflow number and inlet/outlet flow-rate imbalance for a case with an x inflow and a convective x outflow (the X_FACE form of field_set_face_from_field). Three slice reductions on the device and two batched MPI_Allreduces; no D2H copy of the full field. Assumes x is not decomposed, as the FFT Poisson solver requires. Without flow_rate_diff only the outflow number is computed.

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.

Tower loads: constant drag and a lift oscillating at the shedding frequency St*U/D, perturbed by noise (uniform in [0, 1) per element).

Loads on the blades (with end effects) and tower from the element velocities, then the rotor performance.

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

Count the leading steps of a file whose scalar variable_name is no greater than limit. Returns 0 if the variable is not in the file.

No file can be read without ADIOS2, so there are no steps to keep

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

Allocates the solver type of the mesh's boundary condition case and initialises it.

Helper subroutine to create FFT plan with automatic cuFFTMp fallback

Creates the forward and the backward FFT plan the way the solver does. A backward plan that forces the fallback to plain cuFFT makes the forward plan, built on cuFFTMp, get rebuilt on plain cuFFT, so both plans end on the same library. Shared by the solver init and the x3d2-memcheck probe. The optional arguments are the ones of create_fft_plan. label is appended to the plan names in the messages. When stat is present and not zero, the plans made so far are destroyed and the other plan is left undefined.

Allocates the solver type of the mesh's boundary condition case and initialises it.

Standard R2C/C2R plans on the grid permuted into the spectral layout, with cuFFTMp fallback. This serves the 000, 010 and 100 cases. fft_n and fft_n_loc are the plan dimensions from cuda_fft_plan_setup.

R2C with Z-transpose, on (nz, nx, ny). A single rank keeps the plain cuFFT path; on multiple ranks cuFFTMp distributes the transform, and the physical Z-slab is redistributed into its Y-slab input through an MPI_Alltoall (see redistribute_110_to_yslab).

The reason the solver cannot run the BC class bc_code on the cell dimensions cdims and the decomposition nproc_dir, blank when it can. The text is the message the solver init aborts with. Shared by the solver init and x3d2-memcheck.

The cuFFT plan dimensions and types of a BC class. The 100 case transposes X<->Y and runs the FFT on (ny, nx, nz), the 110 case transposes to (nz, nx, ny) with R2C along Z, and the 000 and 010 cases run it on (nx, ny, nz). fft_n holds the global plan dimensions and fft_n_loc the local extents of the first two. Shared by the solver init and the x3d2-memcheck probe.

The size n_spec and the offset n_sp_st of the local spectral pencil in the permuted state, for the BC class bc_code on the cell dimensions dims_glob and dims_loc, the decomposition nproc_dir and the rank position nrank_dir. A BC class the solver does not support gives zeros. Shared by the solver init and the x3d2-memcheck estimate.

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 non-periodic pentadiagonal solve on a single subdomain.

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

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

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.

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

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

Reorder hook of solve_poisson, overridden by every case that has a non-periodic direction.

Gathers the non-periodic x line into a periodic one so that a plain FFT can stand in for the cosine transform.

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.

The per-GPU estimate at num_gpus GPUs. The formula floor stands unless the full-grid probe gave a result and the floor is not over the limit. workspace_gib, overhead_gib and io_gib (optional) split the result. used_gib and npeak (optional) are the num_gpus=1 measurement of the real build and its peak_fields: the estimate is then that measurement scaled to num_gpus, and workspace_gib holds the fields and halo alone.

Interpolate the static polars at Reynolds number re and angle of attack alpha [deg]. Port of CACTUS' EvalStaticCoeff.

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.

The FFT_BC_* class of a poisson_bc_code case, or FFT_BC_UNSUPPORTED. The CUDA solver passes it to the helpers it shares with x3d2-memcheck, so it classifies through poisson_bc_code alone. x3d2-memcheck has no mesh and takes the class from m_mesh's classify_boundaries; test_poisson_bc_code checks that the two agree.

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.

Creates and destroys the forward and backward FFT plans the way the solver does, for the cell dimensions cdims, and for cuFFTMp also the xtdesc buffer, each bracketed with cudaMemGetInfo. worksize_bytes is the plan workspace (used only on the plain cuFFT path, in cuFFTMp it is part of the heap), heap_bytes the memory the plans took, xtdesc_bytes the xtdesc buffer. The 110 BC never uses cuFFTMp. stat is the first non-zero cuFFT/CUDA code, or 0, and the other outputs are then not valid. A plan that aborts inside the library never returns. With query_only set, a cuFFTMp plan is made with auto-allocation off. heap_bytes is then the plan memory without the workspace, worksize_bytes the workspace of both plans (not allocated), xtdesc_bytes is zero and no xtdesc buffer is allocated. A plain cuFFT plan is not changed by it.

Add the constant c to every point of the x plane i_plane of f.

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

Adds isotropic Gaussian sources to the DIR_X VERT fields a, b and c: at every grid point x within distance r_cut of points(1:3, s), a += coeffs(1, s)*exp(-|x - points(:, s)|^2/eps^2), and likewise b and c with coeffs(2, s) and coeffs(3, s). The kernel wraps across periodic boundaries, consistent with field_sample_points, and is cut at non-periodic ones. Needs a uniform grid.

field_add_gaussian_sources_base

field_add_gaussian_sources_base

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

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

f(i_start + i - 1, j, k) = lambda(i)band(i, j, k) + (1 - lambda(i))f(i_start + i - 1, j, k) for a DIR_X vertex field; one thread per i, one block per (j, k)

f(i_start + i - 1, j, k) = lambda(i)band(i, j, k) + (1 - lambda(i))f(i_start + i - 1, j, k) for a DIR_X vertex field

band(i, j, k) = f(i_start + i - 1, j, k) for a DIR_X vertex field; one thread per i, one block per (j, k)

band(i, j, k) = f(i_start + i - 1, j, k) for a DIR_X vertex field

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. Launch with threads = dim3(SZ, 1, 1) and blocks = dim3(n_y_blocks, nz, 1).

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.

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) = sum over z of partial(j, z), in a fixed order. Launch with one thread per y index.

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

Computes this rank's part of the trilinear interpolation of the DIR_X VERT fields a, b and c at points(1:3, p): vals(1:3, p) sums the interpolation weight times a, b and c over the corners of the cell holding point p that this rank owns. Every grid point is owned by one rank, so an MPI sum over ranks gives the full interpolation without a halo exchange. Corners beyond a periodic boundary wrap around; corners outside a non-periodic boundary are skipped. Needs a uniform grid. No MPI reduction is done, so the caller can reduce once for all points.

field_sample_points_base

field_sample_points_base

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). The velocity is interpolated between the sample plane and the one above it.

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). The velocity is interpolated between sample_plane and upper_plane.

Replace one SGS stress plane using the local sampled velocity. The velocity is interpolated linearly between y-vertices sample_plane and sample_plane + 1, with weight sample_weight on the upper one, as Incompact3d samples at a non-integer dsampling.

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

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

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 with a spatially-varying inlet.

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

Write two (ny, nz) planes into the x = 1 and x = nx planes of a DIR_X field. p_end is read only when apply_end = 1. Launch convention matches field_set_x_face_from_field.

Writes two (ny, nz) planes into the x = 1 and x = nx planes of f. p_end is read only when apply_end is true.

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. A thin wrapper on set_field_face_data with face = X_FACE.

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

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

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

Write two (nx, nz) planes into the y = 1 and y = ny planes of a DIR_X field. p_end is read only when apply_end = 1. Launch convention matches field_set_y_face_from_field.

Writes two (nx, nz) planes into the y = 1 and y = ny planes of f. p_end is read only when apply_end is true.

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.

Sets the y plane plane to c.

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

m_base_backend

m_base_backend

Shifts a field by a

Blend an x band of a DIR_X vertex field towards a spanwise-shifted copy of an upstream band (the ABL shifted periodic fringe): f(i + i_offset, j, k) = lambdaf(i, j, k - nshift) + (1 - lambda)f(i + i_offset, j, k) for i = i_start, ..., i_start + size(lambda) - 1, with k circular over the global z planes. x and y must not be decomposed; z_plane_rank(kg) is the rank holding global z plane kg.

m_base_backend

Reduces field to a scalar by integrating over the volume

volume integral of a field

volume integral of a field

The fields and halo term at num_gpus GPUs for npeak fields.

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.

Release a turbine model set up by init_turbine_model; a no-op when the case has none

The per-iteration work after the substeps: sets current_iter and computes the pressure and the derived fields. Shared by run() and x3d2-memcheck.

Source index of entry i of a non-periodic line of n points gathered into a periodic one: the odd entries ascending, then the even entries descending, so that a plain FFT can stand in for the cosine transform. For an odd n the last entry lands at the centre, n/2 + 1. The single definition of the fold for the periodicity kernels: the enforce kernels gather through it and the undo kernels scatter through it, which keeps their reads contiguous. The OpenMP backend has the same function.

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

Global number of cells per direction: a periodic direction has as many cells as vertices, any other has one less. Shared by mesh_init and callers with only the global vertex dims (no mesh_t).

Device memory the GPU-aware ADIOS2 write path (src/io/adios2/io.f90, write_device_field_adios2) holds on top of the solver's own peak: a single reused staging buffer (file_handle%device_staging), sized to one unpadded LOCAL field in whichever precision that write uses. The buffer is grown, never shrunk (see write_device_field_adios2), so its steady-state size is the LARGEST single write it is asked to hold: the build precision whenever a checkpoint is enabled (always native, and always at least as large as a single-precision snapshot), otherwise whatever precision the snapshot uses.

Which writes can take the device path: a snapshot only at unit stride (striding falls back to the host path in snapshot_manager.f90), a checkpoint never strides. Shared by gpu_io_staging_bytes and the x3d2-memcheck report.

The gpu_write_mode_* value of the environment variable X3D2_ADIOS2_GPU_WRITE_MODE. An unset or empty variable gives auto.

Gradient of a scalar field 'p'.

Wrapper for gradient_p2v

CUDA-backend halo-exchange buffer bytes. Mirrors the 24 allocate statements in src/backend/cuda/backend.f90:128-152: 12 arrays shaped (sz, n_halo, n_groups) for u/v/w send/recv start/end, plus 12 more shaped (sz, 1, n_groups) for du/dud/d2u send/recv start/end - all real(dp). n_groups is the max over X/Y/Z of ny_pnz_p/sz, nx_pnz_p/ sz, nx_pny_p/sz (src/allocator.f90:81-83, allocator_t% n_groups_dir/get_n_groups), where _p are the sz-padded dims (x/y padded via padded_dim, z untouched, same rule padded_cells uses).

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

The ibm_.bp mask filename for the periodic flags of the three directions, like "ibm_000.bp" or "ibm_110.bp".

File holding block block (1-based) of a plane database

Sets up the device state of the case ahead of the FFT plans. The 000 case needs nothing beyond the shared state.

if stretching in y is 'bottom'

The 010, 100 and 110 cases pair each dim2 spectral mode with its mirror, which cuFFTMp may place on another rank. Slab index r of P owns global dim2 modes [rm + 1, (r + 1)m] with m = ny_spec, and the mirror of that range is slab P - 1 - r shifted up by one plane, so the extra plane comes from slab P - r. Both partners are symmetric. A single rank holds every partner itself and allocates nothing. This is the only place c_mirror_dev is allocated, see the cuFFTMp check in init_cuda.

Start <line name>.load with its header, unless appending to an existing file on restart.

Allocates c_pair, the y-pencil of sp, and when with_waves is set also waves_pair, redistributing waves into it once.

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

Allocate the turbine model wind_turbine_nml selects, initialise it and register it for checkpointing

Interpolation of a scalar field from cell centres to vertices.

Returns whether or not the current rank is the root rank

True when allocate_flow_case has a branch for the name.

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

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

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.

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.

Laplacian of a scalar field 'u'.

Local one-based index of the zero-based global VERT index g, wrapped across periodic boundaries; owned is false off this rank or outside a non-periodic boundary.

Local one-based index of the zero-based global index g, wrapped when periodic; 0 when it is off this rank or outside the domain.

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.

Element centres, unit vectors, areas and chords from the stations.

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

100-case multi rank only mirror/exchange buffers, zero at num_gpus<=1 - matches the nproc>1 guard at src/backend/cuda/poisson_fft.f90:332-349: c_mirror_dev + c_slab_send_dev: (nx_spec, ny_spec, nz_spec) x2 c_plane_send_dev + c_plane_recv_dev: (nx_spec, nz_spec) x2 all complex(dp). n_spec is the spectral shape of the 100 case.

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 a line of the file opens the namelist group (&group, in any case), so an optional namelist that is present but unterminated can be told apart from an absent one.

Shuts MPI down, then stops.

Number of old-step blocks per variable that a time_intg_t keeps: k - 1 for ABk (nstep - 1) and k for RKk (nstage). Static, so that x3d2-memcheck can size the history blocks without allocating them.

The reason GPU count num_gpus is not a supported decomposition, blank when it is.

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.

Advance the controller to time given the rotor speed [rad/s]. The commanded generator torque is left in gen_trq and the commanded blade pitch angles in pit_com.

The flow_rate_diff for the X_FACE stamp of field_set_face_from_field that makes the flow rate out equal the flow rate in exactly, as the pressure projection sees them. The x staggered derivative sums to (u(nx) - u(1))/dx over the cells, so only the two face planes enter, each weighted over y by y_weights: the case passes the weights of the y interpolation (1/2 on a non-periodic boundary plane, 1 elsewhere), with 0 on any plane it overwrites with the same value on both faces after the stamp. The stamp is linear in u, u(nx) <- (1 - c)u(nx) + cu(nx - 1) + flow_rate_diff, so the sums of the planes before it give the outlet flow rate after it. Call it with u as it is about to be stamped, once bc_start_u_x holds this substep's inlet plane. Assumes x is not decomposed, as compute_outflow_params does.

Whether a snapshot output field is actually driven at runtime - a field only triggers its associated allocations if BOTH it is listed in output_fields AND snapshotting is enabled. Mirrors the pattern used for keep_pressure (src/case/base_case.f90:120-122), shared by run() and x3d2-memcheck.

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

Cell count of one padded field at dims (x, y padded to sz, z untouched). i8 so callers can multiply by field count/nbytes without overflow, as the allocator's own int32 ngrid does for very large grids.

Smallest multiple of sz that is >= n. Used for the x/y padding in allocator_init; z is never padded so is not covered here.

halo_bytes for a local block of local_dims cells: pads x and y to sz (z untouched) the way the allocator does, then sizes the halo-exchange buffers of those padded dims.

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.

Map the text of X3D2_ADIOS2_GPU_WRITE_MODE to a gpu_write_mode_* value. Case is ignored, an unknown or empty text gives auto.

Static estimate of the allocator's work-field high-water mark (allocator_t%next_id after one substep - see src/allocator.f90), evaluated purely from config flags, no build/run required.

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.

Pitch every station by pitch_angle [rad] about the line's root-to-tip direction, then rebuild the element geometry.

Names the boundary condition case of the FFT Poisson solver, one digit per direction in x, y, z order: 0 for periodic, 1 otherwise. The 000 case is periodic everywhere, 110 is non-periodic in x and y.

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

Whether the Poisson solver of this build runs on 2decomp: only the FFT solver of the OMP backend does.

Give every element the airfoil whose thickness ratio is closest to its own, clamped to the thinnest and thickest available.

Called after every completed time step, unlike postprocess which runs only every n_output steps. No-op unless a case overrides it.

Post the sends and receives of a circular z shift by nshift planes for the planes held by other ranks; band(:, :, k) and shifted(:, :, k) are this rank's local plane k. The tag is the receiving global plane.

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.

The NVSHMEM heap of the forward and backward cuFFTMp plans with the workspace allocated and without the xtdesc buffer: the plan memory plus the workspace rounded up to 1 GiB. This is the definition of the measured heap of the normal FFT workspace query.

Decides the FFT probe once per run and runs it. The small-grid probe gives the context and the heap. The full-grid probe follows when the mode allows and the card has room for the memory the plans need. The decision is in ctx%probe_plan and ctx%probe_reason (see ensure_probe in the device).

LU preprocessing for non-periodic pentadiagonal Thomas algorithm.

The verdict line of the requested nproc_dir, or why it is not supported. Sets final_verdict.

A separator line of 60 characters.

One table row. basis is the mark of the overhead (see basis_name).

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.

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.

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 x and y are cosine transformed, so neither has a Nyquist mode and the only singular entry is the mean, which the eps_wave test handles. The x Nyquist zeroing of the 010 case must not be copied here: in this layout x index nx/2 + 1 is the regular cosine mode N/2 (#465).

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. x and y are cosine transformed, so neither has a Nyquist mode and only the mean is singular (handled by eps_wave); the x Nyquist zeroing of the 010 case does not apply here (#465).

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

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.

Rotate point v by theta [rad] about axis through origin (right-hand rule). Equivalent to Incompact3d's quaternion QuatRot, written in Rodrigues' form; axis need not be normalised.

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

Read an .air polar file into airfoil.

Read one block on this rank's y-z subdomain. The arrays are allocated here as (ntimesteps, ny, nz) on the local VERT dims, and times as (ntimesteps). The rank layout may differ from the precursor's, but the global y-z vertex counts must match.

The namelist is optional: a missing inflow_nml keeps the defaults.

Read an actuator line geometry (.al) file: <label>: <length>, <label>: <sx> <sy> <sz>, <label>: <n>, then n rows of r/L c/L pitch[deg] t/c.

Read a string attribute of a file. Returns an empty string if the attribute is not in the file.

No file can be read without ADIOS2, so there is no attribute

Read a &TurbineSpecs namelist file and the airfoils it lists. u_inflow is the initial reference speed of the NREL controller (Incompact3d uses 0.5*(u1 + u2)).

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.

The reason the nproc_dir of the input file is not supported, blank when it is. The CUDA backend only supports nproc_dir = [1,1,num_gpus].

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.

Resolves whether snapshot and checkpoint writes stage through a device buffer, with the environment reader the solver uses for the write mode.

Rotate the stations of line by theta about axis through origin (chordwise vectors about axis through zero). Element geometry is not rebuilt; call make_line_geometry afterwards.

The smallest multiple of step_bytes that is at least n_bytes.

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

The real build: a case build and one substep on 1 GPU in this process, measured with cudaMemGetInfo. The estimate stands when the num_gpus=1 row is DOES_NOT_FIT or the inputs are not valid.

vals(1:3, p) = this rank's part of the trilinear interpolation of a, b and c at points(:, p), one thread per point; see field_sample_points. (dx, dy, dz) is the grid spacing, (nx, ny, nz) and (gx, gy, gz) the local and global VERT dims, (ox, oy, oz) the rank offsets and (px, py, pz) flag periodic directions.

vals(1:3, p) = this rank's part of the trilinear interpolation of the DIR_X VERT fields a, b and c at points(:, p), one thread per point; see field_sample_points. grid holds the local VERT dims, global VERT dims, rank offsets and periodic flags (1/0), three entries each, then the padded y-block count of the DIR_X layout.

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.

Allocates the time integrator field lists and points curr at the solver velocity and species fields. Shared by run() and x3d2-memcheck.

Write host 2D planes into the face planes of a DIR_X field and touch nothing else. For Y_FACE the planes are (nx, nz) and are written to y = 1 (data_start) and y = ny (data_end). For X_FACE the planes are (ny, nz) and are written to x = 1 (data_start) and x = nx (data_end). data_end is optional, when absent only the start face is written.

m_base_backend The host planes are copied to small device arrays and one kernel scatters them into the face planes of the DIR_X field.

m_base_backend

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

Build the blades (and tower) at the current azimuth and collective pitch. Unless restart, the angular velocity is set to uref*tsr/R.

Global z plane that global plane kg takes its value from in a circular shift by nshift planes

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

Gathers each non-periodic line into a periodic one, so that a plain FFT can stand in for the cosine transform, solves in spectral space, and scatters the result back. The 000 case has nothing to reorder and overrides this to transform f in place.

The 000 case has nothing to reorder, so it transforms f in place.

The 000 case has nothing to reorder, so it transforms f in place.

Global zero-based VERT index box [lo, hi] of the grid points within r_cut of each point, on a uniform grid. In a periodic direction the box is left unclipped, so it may run past either end of the domain; the projection wraps those indices onto the grid, which keeps the whole kernel. In a non-periodic direction it is clipped to the domain, so a kernel crossing an inflow, outflow or wall loses that part, as in Incompact3d. Stops the run if a box is wider than a periodic direction, as the wrapped kernel would then cover points twice.

The physical direction held in each dimension of the spectral array. The r2c transform runs along dim1, so a case with x non-periodic permutes a periodic direction into dim1 ahead of the FFT: 100 swaps x and y, 110 brings z to the front.

110-case bytes on top of the waves_dev slab (never runs at num_gpus>1 in this solver): a SECOND complex array of the same spectral shape n_spec, c_dev (src/backend/cuda/poisson_fft.f90:416-418), plus a real, globally-shaped transposed workspace r_dev_110 (:414), which has the same total element count as cdims, just permuted.

Bytes of ONE complex spectral array (waves_dev, and c_dev on the plain-cuFFT fallback path) of the shape n_spec. n_spec comes from cuda_fft_spectral_dims in src/backend/cuda/poisson_fft.f90, the routine the solver init calls, so the shape per BC class and per number of ranks is the solver's. Whether a decomposition is valid is checked by cuda_fft_decomp_reason, not here.

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

010-case, non-uniform y-stretching Poisson coefficient matrices, src/poisson_fft.f90:176-186 (only 010 sets stretched_y, and only at num_gpus=1 since 010 error-stops at nproc>1) and :320-323,424-428 (allocation shapes), copied to the device at src/backend/cuda/poisson_fft.f90:390-411: 'bottom': a_re, a_im, shape (nx_spec, ny_spec, nz_spec, 5), real(dp) otherwise ('centred', 'top-bottom'): a_odd_re/im, a_even_re/im, shape (nx_spec, ny_spec/2, nz_spec, 5), real(dp) doubled by the store_a_*_dev copies when lowmem_fft is .false. (the default).

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

One time-integrator sub-stage: case-specific BCs, transeq, forcings, time integration, case BC re-application, IBM body force, and pressure correction. Extracted from run()'s per-sub-stage loop body so the same allocation path can be driven directly (e.g. by the memory-estimate probe) without going through the full time loop, I/O cadence, or diagnostics.

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.

Bytes to GiB, the one conversion every printed figure goes through.

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.

Inverse of enforce_periodicity.

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.

Fill bc_start_*_x with the precursor plane for the iteration being computed. The plane is held for every substep of the iteration, as Incompact3d does, so it is uploaded once per iteration.

Probe point for the controlled modes, two radii along -rot_n from the rotor centre (as Incompact3d; for a counterclockwise rotor rot_n points upstream, so the probe lies downstream).

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

Local one-based index of the zero-based global VERT index g, wrapped when periodic; 0 when it is off this rank or outside the domain.

Append the element loads of line at iteration iter. points are the element positions the loads were computed at, before the rotor advanced.

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

Sum of the slice f(:, i_slice, :) with the weight w(j) on the points of stacked index j.

Weighted sum of the x slice f(i_slice, :, :) of a DIR_X field: the point with local y index j counts y_weights(j) times. Gives the flow rate through an x face as the discrete divergence sees it, whose y interpolation weights the boundary planes differently from the interior ones. Caller is responsible for MPI_Allreduce across ranks.

x_slice_weighted_sum_base

Weighted sum of the x slice f(:, i_slice, :) of a DIR_X field, w(j) on the points of y index j. One thread per (y, z) point; z is the fast-varying component of the group index.

x_slice_weighted_sum_base

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.