statistics_target.h File Reference

Spatial target resolution for the field-statistics pipeline. More...

#include "variables.h"
#include "field_catalog.h"

Include dependency graph for statistics_target.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Data Structures
struct	SpatialTargetPlan
	Resolved iteration domain for one field on one block. More...

Macros
#define	PICURV_STATISTICS_FLUID_THRESHOLD   0.1
	Threshold below which a cell counts as fluid for the default mask.
#define	PICURV_SPATIAL_AVERAGE_MAX_BUFFER   4000000
	Largest reduction buffer a directional average will allocate, in entries.

Enumerations
enum	PicurvTargetKind { PICURV_TARGET_POINTWISE = 0 }
	Spatial mapping kind. More...
enum	PicurvStatisticsMask { PICURV_STATISTICS_MASK_FLUID = 0 }
	Point-eligibility mask. More...

Functions
PetscErrorCode	PicurvLayoutDimensionIsNodeLike (FieldLayout layout, PetscInt dim, PetscBool *node_like)
	Reports whether a layout is node-like in one dimension.
PetscErrorCode	SpatialTargetPlanCreate (UserCtx *user, FieldId field_id, PicurvStatisticsMask mask, SpatialTargetPlan *plan)
	Resolves the iteration domain for one field on one block.
PetscErrorCode	SpatialTargetPlanLocalPointCount (const SpatialTargetPlan *plan, PetscInt *count)
	Counts the points this rank contributes.
PetscErrorCode	SpatialTargetPlanGlobalPointCount (const SpatialTargetPlan *plan, MPI_Comm comm, PetscInt *count)
	Counts the points contributed across a communicator.
PetscBool	SpatialTargetPlanMaskAllows (const SpatialTargetPlan *plan, PetscReal nvert_value)
	Reports whether a point passes the plan's mask.
PetscErrorCode	PicurvSpatialRatioAverage (UserCtx *user, const SpatialTargetPlan *plan, Vec numerator, Vec denominator, Vec inclusion, const PetscBool average_direction[3], MPI_Comm comm, Vec ratio, PetscReal *scalar)
	Averages two fields over a target domain and divides the results.

Detailed Description

Spatial target resolution for the field-statistics pipeline.

Implements the pointwise identity mapping described in 5. Spatial Target, Layout, and Mask. The abstraction exists with only that one mapping so later spatial bins, profiles, regions, and surfaces extend it rather than retrofit it; see Field Statistics Planned Extensions.

The central responsibility is producing an iteration domain that excludes two distinct categories which Field Identity and Layout Catalog section 4 warns must not be conflated:

1. PETSc/MPI halo entries created by DMDA decomposition, and
2. solver-layout boundary, dummy, or duplicate-periodic indices, which exist even on a single MPI rank.

Nothing here resolves field identity or layout metadata itself; that comes from the typed catalog through FieldGetDescriptor.

Definition in file statistics_target.h.

---

Data Structure Documentation

SpatialTargetPlan

struct SpatialTargetPlan

Resolved iteration domain for one field on one block.

Bounds are half-open [start, end) in DMDA index order (i, j, k) and are already intersected with this rank's owned range, so iterating them touches neither halo storage nor layout boundary/dummy/duplicate-periodic indices. An empty domain on a rank is represented by end[d] <= start[d].

Definition at line 49 of file statistics_target.h.

Collaboration diagram for SpatialTargetPlan:

Data Fields
PicurvTargetKind	kind	Spatial mapping; always pointwise.
PicurvStatisticsMask	mask	Point-eligibility mask.
const FieldDescriptor *	descriptor	Catalog metadata for the targeted field.
PetscInt	start[3]	Inclusive start per dimension (i, j, k).
PetscInt	end[3]	Exclusive end per dimension (i, j, k).

Macro Definition Documentation

PICURV_STATISTICS_FLUID_THRESHOLD

#define PICURV_STATISTICS_FLUID_THRESHOLD   0.1

Threshold below which a cell counts as fluid for the default mask.

Definition at line 29 of file statistics_target.h.

PICURV_SPATIAL_AVERAGE_MAX_BUFFER

#define PICURV_SPATIAL_AVERAGE_MAX_BUFFER   4000000

Largest reduction buffer a directional average will allocate, in entries.

Definition at line 127 of file statistics_target.h.

Enumeration Type Documentation

PicurvTargetKind

enum PicurvTargetKind

Spatial mapping kind.

Only the pointwise identity is implemented.

Enumerator
PICURV_TARGET_POINTWISE

Definition at line 32 of file statistics_target.h.

             {
    PICURV_TARGET_POINTWISE = 0
} PicurvTargetKind;

PicurvStatisticsMask

enum PicurvStatisticsMask

Point-eligibility mask.

Only the fluid mask is implemented.

Enumerator
PICURV_STATISTICS_MASK_FLUID

Definition at line 37 of file statistics_target.h.

             {
    PICURV_STATISTICS_MASK_FLUID = 0
} PicurvStatisticsMask;

Function Documentation

PicurvLayoutDimensionIsNodeLike()

PetscErrorCode PicurvLayoutDimensionIsNodeLike	(	FieldLayout	layout,
		PetscInt	dim,
		PetscBool *	node_like
	)

Reports whether a layout is node-like in one dimension.

Node-like dimensions carry one more physical entry than cell-like dimensions, because they sit on grid lines rather than between them. I_FACE is node-like in i and cell-like in j and k, and correspondingly for the other face families.

Parameters

[in]	layout	Catalog layout.
[in]	dim	Dimension index: 0 for i, 1 for j, 2 for k.
[out]	node_like	Resolved classification.

Returns

Zero on success, or PETSC_ERR_ARG_OUTOFRANGE for an invalid dimension, or PETSC_ERR_ARG_WRONGSTATE for a component-staggered layout, whose components do not share one classification.

Reports whether a layout is node-like in one dimension.

See also

PicurvLayoutDimensionIsNodeLike()

Definition at line 16 of file statistics_target.c.

{
    PetscFunctionBeginUser;
    PetscCheck(node_like != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Classification output is required.");
    PetscCheck(dim >= 0 && dim < 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Dimension must be 0, 1, or 2; got %" PetscInt_FMT ".", dim);
 
    switch (layout) {
    case FIELD_LAYOUT_NODE_CENTERED:
        *node_like = PETSC_TRUE;
        break;
    case FIELD_LAYOUT_CELL_CENTERED:
        *node_like = PETSC_FALSE;
        break;
    case FIELD_LAYOUT_I_FACE:
        *node_like = (PetscBool)(dim == 0);
        break;
    case FIELD_LAYOUT_J_FACE:
        *node_like = (PetscBool)(dim == 1);
        break;
    case FIELD_LAYOUT_K_FACE:
        *node_like = (PetscBool)(dim == 2);
        break;
    case FIELD_LAYOUT_COMPONENT_STAGGERED:
        /* x, y, and z live on I-, J-, and K-faces respectively, so no single
         * classification describes the packed vector. */
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                "Component-staggered layout has no single per-dimension classification.");
    default:
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
                "Unknown field layout %d.", (int)layout);
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SpatialTargetPlanCreate()

PetscErrorCode SpatialTargetPlanCreate	(	UserCtx *	user,
		FieldId	field_id,
		PicurvStatisticsMask	mask,
		SpatialTargetPlan *	plan
	)

Resolves the iteration domain for one field on one block.

Rejects component-staggered fields, whose x, y, and z components live on different face families and therefore cannot share a single pointwise domain. Configured statistics request only cell-centered fields, but the plan resolves node and face layouts correctly so later phases inherit a verified contract.

Parameters

[in]	user	Block context supplying the DMDA layout and periodicity.
[in]	field_id	Catalogued field to target.
[in]	mask	Point-eligibility mask.
[out]	plan	Resolved plan.

Returns

Zero on success, or a PETSc error for a null argument, an unknown field, or an unsupported layout.

Resolves the iteration domain for one field on one block.

See also

SpatialTargetPlanCreate()

Definition at line 79 of file statistics_target.c.

{
    const FieldDescriptor *descriptor = NULL;
    const DMDALocalInfo   *info = NULL;
    SimCtx                *simCtx = NULL;
    PetscInt owned_start[3];
    PetscInt owned_end[3];
    PetscInt global_size[3];
    PetscBool periodic[3];
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && plan != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Block context and plan output are required.");
    PetscCheck(mask == PICURV_STATISTICS_MASK_FLUID, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Only the fluid mask is implemented.");
    simCtx = user->simCtx;
    PetscCheck(simCtx != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "Block context must carry a simulation context.");
 
    PetscCall(FieldGetDescriptor(field_id, &descriptor));
    PetscCheck(descriptor->layout != FIELD_LAYOUT_COMPONENT_STAGGERED,
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "Field '%s' is component-staggered; its components live on different face "
               "families and cannot share one pointwise target domain.",
               descriptor->canonical_name);
 
    info = &user->info;
    owned_start[0] = info->xs; owned_end[0] = info->xs + info->xm; global_size[0] = info->mx;
    owned_start[1] = info->ys; owned_end[1] = info->ys + info->ym; global_size[1] = info->my;
    owned_start[2] = info->zs; owned_end[2] = info->zs + info->zm; global_size[2] = info->mz;
    /* The simCtx flags are authoritative here because they are what selects
     * DM_BOUNDARY_PERIODIC when the DMDA is built, and therefore what decides
     * whether the end planes are wrapped duplicates. */
    periodic[0] = (PetscBool)(simCtx->i_periodic != 0);
    periodic[1] = (PetscBool)(simCtx->j_periodic != 0);
    periodic[2] = (PetscBool)(simCtx->k_periodic != 0);
 
    plan->kind = PICURV_TARGET_POINTWISE;
    plan->mask = mask;
    plan->descriptor = descriptor;
 
    for (PetscInt dim = 0; dim < 3; ++dim) {
        PetscBool node_like = PETSC_FALSE;
        PetscInt  layout_lo = 0;
        PetscInt  layout_hi = 0;
 
        PetscCall(PicurvLayoutDimensionIsNodeLike(descriptor->layout, dim, &node_like));
        ResolveLayoutSpan(node_like, periodic[dim], global_size[dim], &layout_lo, &layout_hi);
 
        /* Intersect the layout-valid span with this rank's owned range: the
         * first exclusion removes solver-layout indices, the second removes
         * PETSc halo storage. */
        plan->start[dim] = PetscMax(owned_start[dim], layout_lo);
        plan->end[dim]   = PetscMin(owned_end[dim], layout_hi);
        if (plan->end[dim] < plan->start[dim]) plan->end[dim] = plan->start[dim];
    }
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SpatialTargetPlanLocalPointCount()

PetscErrorCode SpatialTargetPlanLocalPointCount	(	const SpatialTargetPlan *	plan,
		PetscInt *	count
	)

Counts the points this rank contributes.

Parameters

[in]	plan	Resolved plan.
[out]	count	Local point count; zero for an empty domain.

Returns

Zero on success, or PETSC_ERR_ARG_NULL for a null argument.

Counts the points this rank contributes.

See also

SpatialTargetPlanLocalPointCount()

Definition at line 145 of file statistics_target.c.

{
    PetscInt total = 1;
 
    PetscFunctionBeginUser;
    PetscCheck(plan != NULL && count != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Plan and count output are required.");
    for (PetscInt dim = 0; dim < 3; ++dim) {
        const PetscInt extent = plan->end[dim] - plan->start[dim];
        if (extent <= 0) { *count = 0; PetscFunctionReturn(0); }
        total *= extent;
    }
    *count = total;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SpatialTargetPlanGlobalPointCount()

PetscErrorCode SpatialTargetPlanGlobalPointCount	(	const SpatialTargetPlan *	plan,
		MPI_Comm	comm,
		PetscInt *	count
	)

Counts the points contributed across a communicator.

The result is decomposition independent: it depends only on the layout, global dimensions, and periodicity, never on how ranks divide the domain.

Parameters

[in]	plan	Resolved plan.
[in]	comm	Communicator to reduce over.
[out]	count	Global point count.

Returns

Zero on success, or PETSC_ERR_ARG_NULL for a null argument.

Counts the points contributed across a communicator.

See also

SpatialTargetPlanGlobalPointCount()

Definition at line 167 of file statistics_target.c.

{
    PetscInt local = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(plan != NULL && count != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Plan and count output are required.");
    PetscCall(SpatialTargetPlanLocalPointCount(plan, &local));
    PetscCallMPI(MPI_Allreduce(&local, count, 1, MPIU_INT, MPI_SUM, comm));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SpatialTargetPlanMaskAllows()

PetscBool SpatialTargetPlanMaskAllows	(	const SpatialTargetPlan *	plan,
		PetscReal	nvert_value
	)

Reports whether a point passes the plan's mask.

Because Nvert changes when immersed bodies move, the mask is treated as potentially moving: callers accumulate a per-point valid count and weight rather than assuming a point contributes to every accepted state.

Parameters

[in]	plan	Resolved plan.
[in]	nvert_value	Node-blanking value at the point.

Returns

PETSC_TRUE when the point is eligible.

Reports whether a point passes the plan's mask.

See also

SpatialTargetPlanMaskAllows()

Definition at line 185 of file statistics_target.c.

{
    if (plan == NULL) return PETSC_FALSE;
    return (PetscBool)(nvert_value < PICURV_STATISTICS_FLUID_THRESHOLD);
}

Here is the caller graph for this function:

PicurvSpatialRatioAverage()

PetscErrorCode PicurvSpatialRatioAverage	(	UserCtx *	user,
		const SpatialTargetPlan *	plan,
		Vec	numerator,
		Vec	denominator,
		Vec	inclusion,
		const PetscBool	average_direction[3],
		MPI_Comm	comm,
		Vec	ratio,
		PetscReal *	scalar
	)

Averages two fields over a target domain and divides the results.

Computes ratio = <numerator> / <denominator>, where each average is taken separately over the points the plan admits, and then divides. The order matters: the mean of the pointwise quotients is a different number from the quotient of the means, and several callers need the second one. The LES dynamic procedure is one, because Lilly's least-squares closure is defined that way.

average_direction selects which logical directions collapse into the average. Selecting none makes the operation pointwise, which lets a caller express a local and an averaged formulation through one code path rather than two. Selecting all three yields one number for the whole domain. Selecting a subset retains a profile along the directions left out — averaging over xi and zeta on a plane channel, for example, leaves a wall-normal profile.

Weighting is the caller's, not this function's: scale numerator and denominator by whatever weight the average should carry before calling. A caller wanting a volume-weighted average multiplies both by the cell volume; a caller wanting each admitted point to count once passes them unscaled.

The reduction runs over comm, so a result is independent of how the domain is decomposed across ranks. The plan already excludes PETSc halo storage and the solver's boundary, dummy, and duplicate-periodic indices, so no physical point is counted twice.

Parameters

[in]	user	Block context supplying the DMDA and the solid mask.
[in]	plan	Resolved iteration domain, from SpatialTargetPlanCreate().
[in]	numerator	Ghosted local field to average in the numerator.
[in]	denominator	Ghosted local field to average in the denominator, or NULL to divide by the number of admitted points.
[in]	inclusion	Optional second mask: points where this field is not positive are skipped. Use it to exclude points that hold a zero meaning "never measured" rather than a measurement. NULL applies the plan's mask alone.
[in]	average_direction	Directions to average over, in (xi, eta, zeta) order.
[in]	comm	Communicator to reduce over.
[out]	ratio	Ghosted local field receiving the quotient, or NULL when only the scalar is wanted. Zero wherever the averaged denominator underflows.
[out]	scalar	Receives the single averaged value, or NULL. Valid only when every direction is averaged, since otherwise there is no single value to report.

Returns

Zero on success, PETSC_ERR_ARG_OUTOFRANGE when the retained directions would need an unreasonably large reduction buffer, or PETSC_ERR_ARG_WRONGSTATE when a scalar is requested from a result that still varies in space.

Averages two fields over a target domain and divides the results.

See also

PicurvSpatialRatioAverage()

Definition at line 210 of file statistics_target.c.

{
    DM             da = NULL;
    DMDALocalInfo  info;
    PetscInt       retained_extent[3];
    PetscInt       buffer_size = 1;
    PetscInt       averaged_count = 0;
    PetscReal   ***num = NULL, ***den = NULL, ***inc = NULL, ***out = NULL, ***nvert = NULL;
    PetscReal     *num_sum = NULL, *den_sum = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && plan != NULL && numerator != NULL && average_direction != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context, plan, numerator, and direction selector are required.");
 
    da = user->da;
    PetscCall(DMDAGetLocalInfo(da, &info));
 
    retained_extent[0] = average_direction[0] ? 1 : info.mx;
    retained_extent[1] = average_direction[1] ? 1 : info.my;
    retained_extent[2] = average_direction[2] ? 1 : info.mz;
    for (PetscInt axis = 0; axis < 3; axis++) {
        if (average_direction[axis]) averaged_count++;
        buffer_size *= retained_extent[axis];
    }
    PetscCheck(scalar == NULL || averaged_count == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "A single averaged value exists only when every direction is averaged over; "
               "%" PetscInt_FMT " of 3 were selected.", averaged_count);
 
    PetscCall(DMDAVecGetArrayRead(da, numerator, &num));
    if (denominator) PetscCall(DMDAVecGetArrayRead(da, denominator, &den));
    if (inclusion) PetscCall(DMDAVecGetArrayRead(da, inclusion, &inc));
    if (ratio) PetscCall(DMDAVecGetArray(da, ratio, &out));
    PetscCall(DMDAVecGetArrayRead(da, user->lNvert, &nvert));
 
    if (averaged_count == 0) {
        /* Pointwise: the empty averaging set, which is a legitimate request rather
         * than a degenerate one. It is what a local formulation asks for. The masks
         * mean the same thing here as they do for an averaged set - a point they
         * exclude contributes nothing and receives nothing - so that a caller can
         * switch between the two without its exclusions quietly changing meaning. */
        for (PetscInt k = plan->start[2]; k < plan->end[2]; k++)
        for (PetscInt j = plan->start[1]; j < plan->end[1]; j++)
        for (PetscInt i = plan->start[0]; i < plan->end[0]; i++) {
            const PetscReal divisor = denominator ? den[k][j][i] : 1.0;
 
            if (!SpatialTargetPlanMaskAllows(plan, nvert[k][j][i]) ||
                (inc && inc[k][j][i] <= 0.0)) {
                out[k][j][i] = 0.0;
                continue;
            }
            out[k][j][i] = (PetscAbsReal(divisor) > 0.0) ? (num[k][j][i] / divisor) : 0.0;
        }
    } else {
        PetscCheck(buffer_size <= PICURV_SPATIAL_AVERAGE_MAX_BUFFER, PETSC_COMM_SELF,
                   PETSC_ERR_ARG_OUTOFRANGE,
                   "Averaging over %" PetscInt_FMT " direction(s) would need a reduction buffer "
                   "of %" PetscInt_FMT " entries, above the %d entry limit. Average over more "
                   "directions.", averaged_count, buffer_size, PICURV_SPATIAL_AVERAGE_MAX_BUFFER);
 
        PetscCall(PetscCalloc2(buffer_size, &num_sum, buffer_size, &den_sum));
 
        for (PetscInt k = plan->start[2]; k < plan->end[2]; k++)
        for (PetscInt j = plan->start[1]; j < plan->end[1]; j++)
        for (PetscInt i = plan->start[0]; i < plan->end[0]; i++) {
            PetscInt slot;
 
            if (!SpatialTargetPlanMaskAllows(plan, nvert[k][j][i])) continue;
            /* A point the caller marks unmeasured holds a zero that means absence, not
             * a measurement; counting it would scale the answer toward zero. */
            if (inc && inc[k][j][i] <= 0.0) continue;
 
            slot = SpatialAverageSlot(average_direction, retained_extent, i, j, k);
            num_sum[slot] += num[k][j][i];
            den_sum[slot] += denominator ? den[k][j][i] : 1.0;
        }
 
        PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, num_sum, (PetscMPIInt)buffer_size,
                                   MPIU_REAL, MPI_SUM, comm));
        PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, den_sum, (PetscMPIInt)buffer_size,
                                   MPIU_REAL, MPI_SUM, comm));
 
        if (out) {
            for (PetscInt k = plan->start[2]; k < plan->end[2]; k++)
            for (PetscInt j = plan->start[1]; j < plan->end[1]; j++)
            for (PetscInt i = plan->start[0]; i < plan->end[0]; i++) {
                const PetscInt slot = SpatialAverageSlot(average_direction, retained_extent, i, j, k);
 
                out[k][j][i] = (PetscAbsReal(den_sum[slot]) > 0.0)
                                   ? (num_sum[slot] / den_sum[slot]) : 0.0;
            }
        }
        if (scalar) {
            *scalar = (PetscAbsReal(den_sum[0]) > 0.0) ? (num_sum[0] / den_sum[0]) : 0.0;
        }
 
        PetscCall(PetscFree2(num_sum, den_sum));
    }
 
    PetscCall(DMDAVecRestoreArrayRead(da, numerator, &num));
    if (denominator) PetscCall(DMDAVecRestoreArrayRead(da, denominator, &den));
    if (inclusion) PetscCall(DMDAVecRestoreArrayRead(da, inclusion, &inc));
    if (ratio) PetscCall(DMDAVecRestoreArray(da, ratio, &out));
    PetscCall(DMDAVecRestoreArrayRead(da, user->lNvert, &nvert));
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: