statistics_accumulator.c File Reference

Per-window PETSc accumulator storage and pointwise application. More...

#include "statistics_accumulator.h"
#include "statistics_moments.h"
#include "statistics_target.h"
#include "field_catalog.h"
#include "io.h"
#include "logging.h"
#include <stddef.h>

Include dependency graph for statistics_accumulator.c:

Go to the source code of this file.

Macros
#define	STATISTICS_DERIVED_OUTPUT_LENGTH   256
	Longest output list a recipe may request.
#define	__FUNCT__   "PicurvWindowStorageCreate"
#define	__FUNCT__   "PicurvWindowDerive"
#define	__FUNCT__   "PicurvWindowSpatialMean"
#define	__FUNCT__   "PicurvWindowValidFractionRange"
#define	__FUNCT__   "PicurvWindowAccumulate"

Enumerations
enum	DerivedKind {   DERIVED_MEAN = 0 , DERIVED_REYNOLDS_STRESS , DERIVED_RMS , DERIVED_TKE ,   DERIVED_FLUX , DERIVED_KIND_COUNT }
	Output kinds a postprocessing recipe may request, in enumeration order. More...

Functions
static PetscInt	ProductDiagonalIndex (PetscInt component)
	Internal helper: the product index carrying one component's own variance.
static PetscErrorCode	ProductComponentLabel (PetscInt index, char *out, size_t size)
	Internal helper: writes the two-axis label of one product component.
PetscErrorCode	PicurvProductComponentCount (PetscInt dof, PetscInt *count)
	Implementation of PicurvProductComponentCount().
PetscErrorCode	PicurvCovarianceComponentCount (PetscInt dof_a, PetscInt dof_b, PetscInt *count)
	Implementation of PicurvCovarianceComponentCount().
PetscErrorCode	PicurvStatisticsComponentDM (UserCtx *user, PetscInt components, DM *dm)
	Implementation of PicurvStatisticsComponentDM().
PetscErrorCode	PicurvWindowStorageCreate (UserCtx *user, const PicurvWindowDefinition *definition, PicurvWindowStorage *storage)
	Implementation of PicurvWindowStorageCreate().
PetscErrorCode	PicurvWindowStorageDestroy (PicurvWindowStorage *storage)
	Implementation of PicurvWindowStorageDestroy().
static PetscErrorCode	FindFieldSlot (const PicurvWindowDefinition *definition, PetscInt field_id, PetscInt *slot)
	Internal helper: locates the storage slot holding one field's running mean.
PetscErrorCode	PicurvWindowStoragePayloadCount (const PicurvWindowStorage *storage, PetscInt *count)
	Implementation of PicurvWindowStoragePayloadCount().
PetscErrorCode	PicurvWindowStoragePayload (UserCtx *user, const PicurvWindowDefinition *definition, const PicurvWindowStorage *storage, PetscInt index, PicurvStatisticsPayload *payload)
	Implementation of PicurvWindowStoragePayload().
static PetscErrorCode	ParseDerivedKinds (const char *outputs, PetscBool wanted[DERIVED_KIND_COUNT])
	Internal helper: reports which output kinds a recipe requested.
static PetscErrorCode	DerivedKindExtent (const PicurvWindowDefinition *definition, const PicurvWindowStorage *storage, DerivedKind kind, PetscInt *count)
	Internal helper: counts the derived fields each kind contributes.
PetscErrorCode	PicurvWindowDerivedCount (const PicurvWindowDefinition *definition, const PicurvWindowStorage *storage, const char *outputs, PetscInt *count)
	Implementation of PicurvWindowDerivedCount().
static PetscErrorCode	SafeStandardDeviation (PetscReal variance, const char *label, PetscReal *result)
	Internal helper: takes a square root of a variance that may be barely negative.
static PetscErrorCode	ResolveDerivedIndex (const PicurvWindowDefinition *definition, const PicurvWindowStorage *storage, const PetscBool wanted[DERIVED_KIND_COUNT], PetscInt index, DerivedKind *kind, PetscInt *offset)
	Internal helper: resolves which kind and member one derived index selects.
PetscErrorCode	PicurvWindowDerive (UserCtx *user, const PicurvWindowDefinition *definition, const PicurvWindowStorage *storage, const char *outputs, PetscInt index, Vec scalar_target, Vec vector_target, PicurvDerivedField *field)
	Implementation of PicurvWindowDerive().
PetscErrorCode	PicurvWindowSpatialMean (UserCtx *user, const PicurvWindowDefinition *definition, const PicurvWindowStorage *storage, Vec field, PetscReal *mean)
	Implementation of PicurvWindowSpatialMean().
PetscErrorCode	PicurvWindowValidFractionRange (UserCtx *user, const PicurvWindowDefinition *definition, const PicurvWindowStorage *storage, PetscInt sample_count, PetscReal *minimum, PetscReal *maximum)
	Implementation of PicurvWindowValidFractionRange().
PetscErrorCode	PicurvWindowAccumulate (UserCtx *user, const PicurvWindowDefinition *definition, PicurvWindowStorage *storage, PetscReal weight)
	Implementation of PicurvWindowAccumulate().

Variables
static const PetscInt	kProductFirst [6] = {0, 0, 0, 1, 1, 2}
	Upper-triangular row-major component pairs for a three-vector self-product.
static const PetscInt	kProductSecond [6] = {0, 1, 2, 1, 2, 2}
static const char *const	kAxisName [3] = {"x", "y", "z"}
	Axis labels indexing the pair table above.
static const char *const	kDerivedKindName [DERIVED_KIND_COUNT]
	Recipe spellings of the output kinds.
static const PetscInt	kDerivedKindScaleExponent [DERIVED_KIND_COUNT]
	Power of the source field's reference scale each derived kind carries.

Detailed Description

Per-window PETSc accumulator storage and pointwise application.

Full API contract is documented with the declarations in include/statistics_accumulator.h.

Definition in file statistics_accumulator.c.

Macro Definition Documentation

STATISTICS_DERIVED_OUTPUT_LENGTH

#define STATISTICS_DERIVED_OUTPUT_LENGTH   256

Longest output list a recipe may request.

Definition at line 19 of file statistics_accumulator.c.

__FUNCT__ [1/5]

#define __FUNCT__   "PicurvWindowStorageCreate"

Definition at line 135 of file statistics_accumulator.c.

__FUNCT__ [2/5]

#define __FUNCT__   "PicurvWindowDerive"

Definition at line 135 of file statistics_accumulator.c.

__FUNCT__ [3/5]

#define __FUNCT__   "PicurvWindowSpatialMean"

Definition at line 135 of file statistics_accumulator.c.

__FUNCT__ [4/5]

#define __FUNCT__   "PicurvWindowValidFractionRange"

Definition at line 135 of file statistics_accumulator.c.

__FUNCT__ [5/5]

#define __FUNCT__   "PicurvWindowAccumulate"

Definition at line 135 of file statistics_accumulator.c.

Enumeration Type Documentation

DerivedKind

enum DerivedKind

Output kinds a postprocessing recipe may request, in enumeration order.

Enumerator
DERIVED_MEAN
DERIVED_REYNOLDS_STRESS
DERIVED_RMS
DERIVED_TKE
DERIVED_FLUX
DERIVED_KIND_COUNT

Definition at line 380 of file statistics_accumulator.c.

             {
    DERIVED_MEAN = 0,
    DERIVED_REYNOLDS_STRESS,
    DERIVED_RMS,
    DERIVED_TKE,
    DERIVED_FLUX,
    DERIVED_KIND_COUNT
} DerivedKind;

Function Documentation

ProductDiagonalIndex()

static PetscInt ProductDiagonalIndex ( PetscInt  component )

static

Internal helper: the product index carrying one component's own variance.

Local to this translation unit. Found by searching the pair table for the entry pairing a component with itself, so the diagonal cannot be stated separately from the order it belongs to.

Definition at line 56 of file statistics_accumulator.c.

{
    for (PetscInt c = 0; c < 6; ++c) {
        if (kProductFirst[c] == component && kProductSecond[c] == component) return c;
    }
    return -1;
}

Here is the caller graph for this function:

ProductComponentLabel()

static PetscErrorCode ProductComponentLabel ( PetscInt  index, char *  out, size_t  size  )

static

Internal helper: writes the two-axis label of one product component.

Local to this translation unit. Built from the pair table, so "xy" and the slot it names can never refer to different pairs.

Definition at line 69 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    PetscCheck(index >= 0 && index < 6, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Product component %" PetscInt_FMT " is outside the symmetric set.", index);
    PetscCall(PetscSNPrintf(out, size, "%s%s", kAxisName[kProductFirst[index]],
                            kAxisName[kProductSecond[index]]));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvProductComponentCount()

PetscErrorCode PicurvProductComponentCount	(	PetscInt	dof,
		PetscInt *	count
	)

Implementation of PicurvProductComponentCount().

Reports how many symmetric product components a field's second moment needs.

See also

PicurvProductComponentCount()

Definition at line 83 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    PetscCheck(count != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Count output is required.");
    PetscCheck(dof == 1 || dof == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Self-products are supported for scalar and three-vector fields, got dof %" PetscInt_FMT ".", dof);
    *count = (dof == 1) ? 1 : 6;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvCovarianceComponentCount()

PetscErrorCode PicurvCovarianceComponentCount	(	PetscInt	dof_a,
		PetscInt	dof_b,
		PetscInt *	count
	)

Implementation of PicurvCovarianceComponentCount().

Reports how many components a covariance between two fields needs.

See also

PicurvCovarianceComponentCount()

Definition at line 97 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    PetscCheck(count != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Count output is required.");
    /* Scalar-scalar and vector-scalar pairs only; a vector-vector cross product
     * would need a nine-component carrier that nothing allocates. */
    PetscCheck((dof_a == 1 && dof_b == 1) || (dof_a == 3 && dof_b == 1) || (dof_a == 1 && dof_b == 3),
               PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Covariance is supported for scalar-scalar and vector-scalar pairs, got dof %" PetscInt_FMT
               " and %" PetscInt_FMT ".", dof_a, dof_b);
    *count = (dof_a == 3 || dof_b == 3) ? 3 : 1;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvStatisticsComponentDM()

PetscErrorCode PicurvStatisticsComponentDM	(	UserCtx *	user,
		PetscInt	components,
		DM *	dm
	)

Implementation of PicurvStatisticsComponentDM().

Resolves the DM carrying a given number of accumulator components.

See also

PicurvStatisticsComponentDM()

Definition at line 115 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    PetscCheck(user != NULL && dm != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context and DM output are required.");
    switch (components) {
    case 1: *dm = user->da;   break;
    case 3: *dm = user->fda;  break;
    case 6: *dm = user->fda6; break;
    default:
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
                "No block DM carries %" PetscInt_FMT " accumulator components.", components);
    }
    PetscCheck(*dm != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
               "The block DM for %" PetscInt_FMT " accumulator components was never created.",
               components);
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvWindowStorageCreate()

PetscErrorCode PicurvWindowStorageCreate	(	UserCtx *	user,
		const PicurvWindowDefinition *	definition,
		PicurvWindowStorage *	storage
	)

Implementation of PicurvWindowStorageCreate().

Allocates the accumulator state one window owns on one block.

See also

PicurvWindowStorageCreate()

Definition at line 140 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && definition != NULL && storage != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Context, definition, and storage are required.");
    PetscCall(PetscMemzero(storage, sizeof(*storage)));
    storage->field_count = definition->field_count;
    storage->covariance_count = definition->covariance_count;
 
    /* Occupancy is per point and scalar regardless of what is accumulated. */
    PetscCall(DMCreateGlobalVector(user->da, &storage->count));
    PetscCall(VecSet(storage->count, 0.0));
    PetscCall(DMCreateGlobalVector(user->da, &storage->weight));
    PetscCall(VecSet(storage->weight, 0.0));
    PetscCall(DMCreateGlobalVector(user->da, &storage->weight_sq));
    PetscCall(VecSet(storage->weight_sq, 0.0));
 
    if (definition->field_count > 0) {
        PetscCall(PetscCalloc1((size_t)definition->field_count, &storage->mean));
        PetscCall(PetscCalloc1((size_t)definition->field_count, &storage->m2));
    }
    for (PetscInt field_index = 0; field_index < definition->field_count; ++field_index) {
        FieldView view;
        PetscInt  components = 0;
        DM        product_dm = NULL;
 
        PetscCall(FieldGetView(user, (FieldId)definition->fields[field_index].field_id, &view));
        /* The mean matches the source field's own layout, so it comes from that
         * field's DM and inherits its decomposition. */
        PetscCall(DMCreateGlobalVector(view.dm, &storage->mean[field_index]));
        PetscCall(VecSet(storage->mean[field_index], 0.0));
        if (!definition->fields[field_index].want_second) continue;
        PetscCall(PicurvProductComponentCount(view.descriptor->dof, &components));
        PetscCall(PicurvStatisticsComponentDM(user, components, &product_dm));
        PetscCall(DMCreateGlobalVector(product_dm, &storage->m2[field_index]));
        PetscCall(VecSet(storage->m2[field_index], 0.0));
    }
 
    if (definition->covariance_count > 0) {
        PetscCall(PetscCalloc1((size_t)definition->covariance_count, &storage->cm));
    }
    for (PetscInt pair_index = 0; pair_index < definition->covariance_count; ++pair_index) {
        FieldView view_a, view_b;
        PetscInt  components = 0;
        DM        pair_dm = NULL;
 
        PetscCall(FieldGetView(user, (FieldId)definition->covariances[pair_index].first, &view_a));
        PetscCall(FieldGetView(user, (FieldId)definition->covariances[pair_index].second, &view_b));
        PetscCall(PicurvCovarianceComponentCount(view_a.descriptor->dof, view_b.descriptor->dof, &components));
        PetscCall(PicurvStatisticsComponentDM(user, components, &pair_dm));
        PetscCall(DMCreateGlobalVector(pair_dm, &storage->cm[pair_index]));
        PetscCall(VecSet(storage->cm[pair_index], 0.0));
    }
 
    /* Report what this window costs and what it covers, once per block at setup.
     * The point count is a function of layout, dimensions, and periodicity alone, so
     * it tells an operator what "per point" will mean before any sample is taken. */
    {
        SpatialTargetPlan plan;
        PetscInt payloads = 0;
        PetscInt points = 0;
 
        PetscCall(PicurvWindowStoragePayloadCount(storage, &payloads));
        PetscCall(SpatialTargetPlanCreate(user, (FieldId)definition->fields[0].field_id,
                                          PICURV_STATISTICS_MASK_FLUID, &plan));
        PetscCall(SpatialTargetPlanGlobalPointCount(&plan, PETSC_COMM_WORLD, &points));
        LOG_ALLOW(LOCAL, LOG_DEBUG,
                  "Statistics window '%s' block %d: %d accumulator vector(s) over %d point(s).\n",
                  definition->name, (int)user->_this, (int)payloads, (int)points);
    }
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PicurvWindowStorageDestroy()

PetscErrorCode PicurvWindowStorageDestroy

(

PicurvWindowStorage *

storage

)

Implementation of PicurvWindowStorageDestroy().

Releases accumulator state previously created for one window.

See also

PicurvWindowStorageDestroy()

Definition at line 220 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    if (storage == NULL) PetscFunctionReturn(0);
    if (storage->count) PetscCall(VecDestroy(&storage->count));
    if (storage->weight) PetscCall(VecDestroy(&storage->weight));
    if (storage->weight_sq) PetscCall(VecDestroy(&storage->weight_sq));
    for (PetscInt field_index = 0; storage->mean && field_index < storage->field_count; ++field_index) {
        if (storage->mean[field_index]) PetscCall(VecDestroy(&storage->mean[field_index]));
    }
    for (PetscInt field_index = 0; storage->m2 && field_index < storage->field_count; ++field_index) {
        if (storage->m2[field_index]) PetscCall(VecDestroy(&storage->m2[field_index]));
    }
    for (PetscInt pair_index = 0; storage->cm && pair_index < storage->covariance_count; ++pair_index) {
        if (storage->cm[pair_index]) PetscCall(VecDestroy(&storage->cm[pair_index]));
    }
    if (storage->mean) PetscCall(PetscFree(storage->mean));
    if (storage->m2) PetscCall(PetscFree(storage->m2));
    if (storage->cm) PetscCall(PetscFree(storage->cm));
    PetscCall(PetscMemzero(storage, sizeof(*storage)));
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

FindFieldSlot()

static PetscErrorCode FindFieldSlot ( const PicurvWindowDefinition *  definition, PetscInt  field_id, PetscInt *  slot  )

static

Internal helper: locates the storage slot holding one field's running mean.

Local to this translation unit. A covariance member must also appear in the window's field list, because the co-moment update needs that field's running mean; this is where that requirement is enforced.

Definition at line 249 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    for (PetscInt field_index = 0; field_index < definition->field_count; ++field_index) {
        if (definition->fields[field_index].field_id == field_id) {
            *slot = field_index;
            PetscFunctionReturn(0);
        }
    }
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
            "Window '%s' requests a covariance over field '%s', which is not in its field list.",
            definition->name, FieldCanonicalName((FieldId)field_id));
}

Here is the call graph for this function:

Here is the caller graph for this function:

PicurvWindowStoragePayloadCount()

PetscErrorCode PicurvWindowStoragePayloadCount	(	const PicurvWindowStorage *	storage,
		PetscInt *	count
	)

Implementation of PicurvWindowStoragePayloadCount().

Reports how many checkpointable vectors one window's storage holds.

See also

PicurvWindowStoragePayloadCount()

Definition at line 268 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    PetscCheck(storage != NULL && count != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Storage and count output are required.");
    *count = 3 + storage->field_count + storage->covariance_count;
    /* A field without a requested second moment holds no product vector, so the
     * enumeration skips it rather than emitting an empty payload. */
    for (PetscInt field_index = 0; storage->m2 && field_index < storage->field_count; ++field_index) {
        if (storage->m2[field_index]) *count += 1;
    }
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvWindowStoragePayload()

PetscErrorCode PicurvWindowStoragePayload	(	UserCtx *	user,
		const PicurvWindowDefinition *	definition,
		const PicurvWindowStorage *	storage,
		PetscInt	index,
		PicurvStatisticsPayload *	payload
	)

Implementation of PicurvWindowStoragePayload().

Resolves one enumerated payload of a window's storage.

See also

PicurvWindowStoragePayload()

Definition at line 286 of file statistics_accumulator.c.

{
    PetscInt total = 0;
    PetscInt cursor = index;
 
    PetscFunctionBeginUser;
    PetscCheck(user != NULL && definition != NULL && storage != NULL && payload != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context, definition, storage, and payload output are required.");
    PetscCall(PicurvWindowStoragePayloadCount(storage, &total));
    PetscCheck(index >= 0 && index < total, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Payload index %" PetscInt_FMT " is outside [0, %" PetscInt_FMT ").", index, total);
    PetscCall(PetscMemzero(payload, sizeof(*payload)));
 
    /* Occupancy first, then means, then products, then co-moments. The order is the
     * persistence contract: the manifest inventory, the writer, and the reader all
     * walk it identically. */
    if (cursor < 3) {
        static const char *const occupancy_name[3] = {"count", "weight", "weight_sq"};
        Vec occupancy[3];
 
        occupancy[0] = storage->count;
        occupancy[1] = storage->weight;
        occupancy[2] = storage->weight_sq;
        PetscCall(PetscStrncpy(payload->name, occupancy_name[cursor], sizeof(payload->name)));
        payload->vec = occupancy[cursor];
        payload->components = 1;
        payload->role = "occupancy";
        payload->layout = FieldLayoutName(FIELD_LAYOUT_CELL_CENTERED);
        PetscFunctionReturn(0);
    }
    cursor -= 3;
 
    if (cursor < storage->field_count) {
        const FieldDescriptor *descriptor = NULL;
 
        PetscCall(FieldGetDescriptor((FieldId)definition->fields[cursor].field_id, &descriptor));
        PetscCall(PetscSNPrintf(payload->name, sizeof(payload->name), "%s_mean",
                                descriptor->canonical_name));
        payload->vec = storage->mean[cursor];
        payload->components = descriptor->dof;
        payload->role = "mean";
        payload->layout = FieldLayoutName(descriptor->layout);
        PetscFunctionReturn(0);
    }
    cursor -= storage->field_count;
 
    {
        PetscInt product_count = 0;
        PetscInt seen = 0;
 
        for (PetscInt field_index = 0; storage->m2 && field_index < storage->field_count; ++field_index) {
            if (storage->m2[field_index]) ++product_count;
        }
        if (cursor < product_count) {
            for (PetscInt field_index = 0; field_index < storage->field_count; ++field_index) {
                const FieldDescriptor *descriptor = NULL;
 
                if (!storage->m2[field_index]) continue;
                if (seen++ != cursor) continue;
                PetscCall(FieldGetDescriptor((FieldId)definition->fields[field_index].field_id, &descriptor));
                PetscCall(PetscSNPrintf(payload->name, sizeof(payload->name), "%s_m2",
                                        descriptor->canonical_name));
                payload->vec = storage->m2[field_index];
                PetscCall(PicurvProductComponentCount(descriptor->dof, &payload->components));
                payload->role = "second_moment";
                payload->layout = FieldLayoutName(descriptor->layout);
                PetscFunctionReturn(0);
            }
        }
        cursor -= product_count;
    }
 
    {
        const FieldDescriptor *first = NULL;
        const FieldDescriptor *second = NULL;
 
        PetscCheck(cursor >= 0 && cursor < storage->covariance_count, PETSC_COMM_SELF, PETSC_ERR_PLIB,
                   "Payload index %" PetscInt_FMT " fell through the storage enumeration.", index);
        PetscCall(FieldGetDescriptor((FieldId)definition->covariances[cursor].first, &first));
        PetscCall(FieldGetDescriptor((FieldId)definition->covariances[cursor].second, &second));
        PetscCall(PetscSNPrintf(payload->name, sizeof(payload->name), "%s_%s_cm",
                                first->canonical_name, second->canonical_name));
        payload->vec = storage->cm[cursor];
        PetscCall(PicurvCovarianceComponentCount(first->dof, second->dof, &payload->components));
        payload->role = "co_moment";
        payload->layout = FieldLayoutName(first->layout);
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

ParseDerivedKinds()

static PetscErrorCode ParseDerivedKinds ( const char *  outputs, PetscBool  wanted[DERIVED_KIND_COUNT]  )

static

Internal helper: reports which output kinds a recipe requested.

Local to this translation unit.

Definition at line 417 of file statistics_accumulator.c.

{
    char buffer[STATISTICS_DERIVED_OUTPUT_LENGTH];
    char *cursor = buffer;
 
    PetscFunctionBeginUser;
    for (PetscInt k = 0; k < DERIVED_KIND_COUNT; ++k) wanted[k] = PETSC_FALSE;
    PetscCheck(outputs != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Output list is required.");
    PetscCall(PetscStrncpy(buffer, outputs, sizeof(buffer)));
 
    while (cursor && *cursor) {
        char *comma = strchr(cursor, ',');
        PetscBool matched = PETSC_FALSE;
 
        if (comma) *comma = '\0';
        TrimWhitespace(cursor);
        if (cursor[0] != '\0') {
            for (PetscInt k = 0; k < DERIVED_KIND_COUNT; ++k) {
                if (strcmp(cursor, kDerivedKindName[k])) continue;
                wanted[k] = PETSC_TRUE;
                matched = PETSC_TRUE;
                break;
            }
            PetscCheck(matched, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                       "Unknown field-statistics output '%s'. Available outputs are "
                       "mean, reynolds_stress, rms, tke, and flux.", cursor);
        }
        cursor = comma ? comma + 1 : NULL;
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

DerivedKindExtent()

static PetscErrorCode DerivedKindExtent ( const PicurvWindowDefinition *  definition, const PicurvWindowStorage *  storage, DerivedKind  kind, PetscInt *  count  )

static

Internal helper: counts the derived fields each kind contributes.

Local to this translation unit. A kind contributes nothing when the state it needs was never accumulated, so a recipe may name every output without having to know which window carries which moment.

Definition at line 455 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    *count = 0;
    switch (kind) {
    case DERIVED_MEAN:
        *count = storage->field_count;
        break;
    case DERIVED_REYNOLDS_STRESS:
    case DERIVED_RMS:
        for (PetscInt field_index = 0; storage->m2 && field_index < storage->field_count; ++field_index) {
            const FieldDescriptor *descriptor = NULL;
            PetscInt components = 0;
 
            if (!storage->m2[field_index]) continue;
            PetscCall(FieldGetDescriptor((FieldId)definition->fields[field_index].field_id, &descriptor));
            PetscCall(PicurvProductComponentCount(descriptor->dof, &components));
            /* A stress tensor emits every component; an RMS emits only the
             * diagonal, because an off-diagonal has no square root to take. */
            *count += (kind == DERIVED_REYNOLDS_STRESS) ? components : descriptor->dof;
        }
        break;
    case DERIVED_TKE:
        /* Turbulent kinetic energy is the trace of a three-vector's stress tensor,
         * so it exists only for a vector field carrying a second moment. */
        for (PetscInt field_index = 0; storage->m2 && field_index < storage->field_count; ++field_index) {
            const FieldDescriptor *descriptor = NULL;
 
            if (!storage->m2[field_index]) continue;
            PetscCall(FieldGetDescriptor((FieldId)definition->fields[field_index].field_id, &descriptor));
            if (descriptor->dof == 3) *count += 1;
        }
        break;
    default:
        *count = storage->covariance_count;
        break;
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PicurvWindowDerivedCount()

PetscErrorCode PicurvWindowDerivedCount	(	const PicurvWindowDefinition *	definition,
		const PicurvWindowStorage *	storage,
		const char *	outputs,
		PetscInt *	count
	)

Implementation of PicurvWindowDerivedCount().

Reports how many derived fields a requested output set produces.

See also

PicurvWindowDerivedCount()

Definition at line 501 of file statistics_accumulator.c.

{
    PetscBool wanted[DERIVED_KIND_COUNT];
 
    PetscFunctionBeginUser;
    PetscCheck(definition != NULL && storage != NULL && count != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Definition, storage, and count are required.");
    PetscCall(ParseDerivedKinds(outputs, wanted));
    *count = 0;
    for (PetscInt k = 0; k < DERIVED_KIND_COUNT; ++k) {
        PetscInt extent = 0;
 
        if (!wanted[k]) continue;
        PetscCall(DerivedKindExtent(definition, storage, (DerivedKind)k, &extent));
        *count += extent;
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SafeStandardDeviation()

static PetscErrorCode SafeStandardDeviation ( PetscReal  variance, const char *  label, PetscReal *  result  )

static

Internal helper: takes a square root of a variance that may be barely negative.

Local to this translation unit. Centered accumulation can leave a variance a few ulps below zero when the signal is nearly constant. Clamping is confined to this one place, applies only under a root, and never touches stored state; a genuinely negative variance is a defect and is reported.

Definition at line 529 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    if (variance >= 0.0) {
        *result = PetscSqrtReal(variance);
        PetscFunctionReturn(0);
    }
    PetscCheck(variance >= -PICURV_STATISTICS_VARIANCE_FLOOR, PETSC_COMM_SELF, PETSC_ERR_FP,
               "Derived variance for '%s' is %g, which is too negative to be floating-point "
               "cancellation; the accumulated state is inconsistent.", label, (double)variance);
    *result = 0.0;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

ResolveDerivedIndex()

static PetscErrorCode ResolveDerivedIndex ( const PicurvWindowDefinition *  definition, const PicurvWindowStorage *  storage, const PetscBool  wanted[DERIVED_KIND_COUNT], PetscInt  index, DerivedKind *  kind, PetscInt *  offset  )

static

Internal helper: resolves which kind and member one derived index selects.

Local to this translation unit. Walks the same order DerivedKindExtent counts, so the enumeration and the count cannot disagree.

Definition at line 548 of file statistics_accumulator.c.

{
    PetscFunctionBeginUser;
    for (PetscInt k = 0; k < DERIVED_KIND_COUNT; ++k) {
        PetscInt extent = 0;
 
        if (!wanted[k]) continue;
        PetscCall(DerivedKindExtent(definition, storage, (DerivedKind)k, &extent));
        if (index < extent) {
            *kind = (DerivedKind)k;
            *offset = index;
            PetscFunctionReturn(0);
        }
        index -= extent;
    }
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
            "Derived index is past the end of the requested output set.");
}

Here is the call graph for this function:

Here is the caller graph for this function:

PicurvWindowDerive()

PetscErrorCode PicurvWindowDerive	(	UserCtx *	user,
		const PicurvWindowDefinition *	definition,
		const PicurvWindowStorage *	storage,
		const char *	outputs,
		PetscInt	index,
		Vec	scalar_target,
		Vec	vector_target,
		PicurvDerivedField *	field
	)

Implementation of PicurvWindowDerive().

Derives one output field from centered accumulator state.

See also

PicurvWindowDerive()

Definition at line 576 of file statistics_accumulator.c.

{
    PetscBool wanted[DERIVED_KIND_COUNT];
    DerivedKind kind = DERIVED_MEAN;
    SpatialTargetPlan plan;
    PetscReal ***weight_arr = NULL, ***weight_sq_arr = NULL, ***count_arr = NULL;
    PetscScalar ****source = NULL, ****target = NULL;
    const FieldDescriptor *descriptor = NULL;
    DM source_dm = NULL;
    Vec source_vec = NULL;
    PetscInt offset = 0, slot = 0, member = 0;
    /* Resolved where the branch knows which fields it read; 1.0 leaves the result in
     * solver units, which is also what an undeclared scale falls back to. */
    PetscReal dimensional_factor = 1.0;
    /* The source's component count is not the output's: a scalar RMS reads a
     * six-component tensor, so each needs its own DM. */
    PetscInt source_components = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && definition != NULL && storage != NULL && field != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context, definition, storage, and output are required.");
    PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Derived index %" PetscInt_FMT " is negative.", index);
    PetscCall(PetscMemzero(field, sizeof(*field)));
    /* The output list is parsed once and the extents walked once. Resolving the
     * index is what bounds-checks it, so counting first would repeat both. */
    PetscCall(ParseDerivedKinds(outputs, wanted));
    PetscCall(ResolveDerivedIndex(definition, storage, wanted, index, &kind, &offset));
 
    /* Locate the accumulator slot and component this index selects. */
    switch (kind) {
    case DERIVED_MEAN:
        slot = offset;
        PetscCall(FieldGetDescriptor((FieldId)definition->fields[slot].field_id, &descriptor));
        source_vec = storage->mean[slot];
        field->components = descriptor->dof;
        source_components = descriptor->dof;
        PetscCall(PetscSNPrintf(field->name, sizeof(field->name), "%s_%s_mean",
                                definition->name, descriptor->canonical_name));
        break;
    case DERIVED_TKE:
        for (slot = 0; slot < storage->field_count; ++slot) {
            if (!storage->m2[slot]) continue;
            PetscCall(FieldGetDescriptor((FieldId)definition->fields[slot].field_id, &descriptor));
            if (descriptor->dof != 3) continue;
            if (offset-- == 0) break;
        }
        source_vec = storage->m2[slot];
        field->components = 1;
        PetscCall(PicurvProductComponentCount(descriptor->dof, &source_components));
        PetscCall(PetscSNPrintf(field->name, sizeof(field->name), "%s_%s_tke",
                                definition->name, descriptor->canonical_name));
        break;
    case DERIVED_FLUX:
        slot = offset;
        {
            const FieldDescriptor *first = NULL;
            const FieldDescriptor *second = NULL;
 
            PetscCall(FieldGetDescriptor((FieldId)definition->covariances[slot].first, &first));
            PetscCall(FieldGetDescriptor((FieldId)definition->covariances[slot].second, &second));
            PetscCall(PicurvCovarianceComponentCount(first->dof, second->dof, &field->components));
            source_components = field->components;
            source_vec = storage->cm[slot];
            if (user->simCtx->pps && user->simCtx->pps->dimensionalize) {
                PetscReal first_scale = 1.0, second_scale = 1.0;
 
                /* A co-moment carries the product of its two fields' dimensions. */
                PetscCall(PicurvFieldReferenceScale(user->simCtx, first->canonical_name,
                                                    &first_scale, NULL, 0));
                PetscCall(PicurvFieldReferenceScale(user->simCtx, second->canonical_name,
                                                    &second_scale, NULL, 0));
                dimensional_factor = first_scale * second_scale;
            }
            PetscCall(PetscSNPrintf(field->name, sizeof(field->name), "%s_%s_%s_flux",
                                    definition->name, first->canonical_name, second->canonical_name));
        }
        break;
    default:
        /* Stress and RMS both walk the fields carrying a product; stress emits every
         * symmetric component, RMS only the diagonal. */
        for (slot = 0; slot < storage->field_count; ++slot) {
            PetscInt extent = 0;
 
            if (!storage->m2[slot]) continue;
            PetscCall(FieldGetDescriptor((FieldId)definition->fields[slot].field_id, &descriptor));
            PetscCall(PicurvProductComponentCount(descriptor->dof, &extent));
            if (kind == DERIVED_RMS) extent = descriptor->dof;
            if (offset < extent) { member = offset; break; }
            offset -= extent;
        }
        source_vec = storage->m2[slot];
        field->components = 1;
        PetscCall(PicurvProductComponentCount(descriptor->dof, &source_components));
        if (kind == DERIVED_RMS) {
            PetscCall(PetscSNPrintf(field->name, sizeof(field->name), "%s_%s_rms%s",
                                    definition->name, descriptor->canonical_name,
                                    descriptor->dof == 1 ? "" : kAxisName[member]));
        } else if (descriptor->dof == 1) {
            PetscCall(PetscSNPrintf(field->name, sizeof(field->name), "%s_%s_variance",
                                    definition->name, descriptor->canonical_name));
        } else {
            char label[8];
 
            PetscCall(ProductComponentLabel(member, label, sizeof(label)));
            PetscCall(PetscSNPrintf(field->name, sizeof(field->name), "%s_%s_R_%s",
                                    definition->name, descriptor->canonical_name, label));
        }
        break;
    }
 
    PetscCall(SpatialTargetPlanCreate(user, (FieldId)definition->fields[0].field_id,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    PetscCall(PicurvStatisticsComponentDM(user, source_components, &source_dm));
    {
        Vec destination = (field->components == 1) ? scalar_target : vector_target;
        DM destination_dm = NULL;
 
        PetscCheck(destination != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
                   "Derived field '%s' needs a %d-component destination.",
                   field->name, (int)field->components);
        PetscCall(PicurvStatisticsComponentDM(user, field->components, &destination_dm));
        PetscCall(VecZeroEntries(destination));
 
        PetscCall(DMDAVecGetArrayRead(user->da, storage->weight, &weight_arr));
        PetscCall(DMDAVecGetArrayRead(user->da, storage->weight_sq, &weight_sq_arr));
        PetscCall(DMDAVecGetArrayRead(user->da, storage->count, &count_arr));
        PetscCall(DMDAVecGetArrayDOFRead(source_dm, source_vec, &source));
        PetscCall(DMDAVecGetArrayDOF(destination_dm, destination, &target));
 
        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) {
                    PicurvCoMomentState pair;
 
                    /* A point the window never sampled has no average at all, so it
                     * is left at zero rather than divided by a zero weight. */
                    if (weight_arr[k][j][i] <= 0.0) continue;
                    pair.count = count_arr[k][j][i];
                    pair.weight = weight_arr[k][j][i];
                    pair.weight_sq = weight_sq_arr[k][j][i];
                    pair.mean_x = 0.0;
                    pair.mean_y = 0.0;
 
                    if (kind == DERIVED_MEAN) {
                        for (PetscInt c = 0; c < field->components; ++c) {
                            target[k][j][i][c] = source[k][j][i][c];
                        }
                    } else if (kind == DERIVED_TKE) {
                        /* The trace of the stress tensor, halved: components 0, 3
                         * and 5 are xx, yy and zz in the stored symmetric order. */
                        PetscReal trace = 0.0;
 
                        for (PetscInt c = 0; c < 3; ++c) {
                            pair.cm = source[k][j][i][ProductDiagonalIndex(c)];
                            trace += PicurvCoMomentStateCovariance(&pair);
                        }
                        target[k][j][i][0] = 0.5 * trace;
                    } else if (kind == DERIVED_RMS) {
                        PetscReal deviation = 0.0;
 
                        pair.cm = source[k][j][i][(descriptor->dof == 1)
                                                  ? 0 : ProductDiagonalIndex(member)];
                        PetscCall(SafeStandardDeviation(PicurvCoMomentStateCovariance(&pair),
                                                        field->name, &deviation));
                        target[k][j][i][0] = deviation;
                    } else if (kind == DERIVED_FLUX) {
                        for (PetscInt c = 0; c < field->components; ++c) {
                            pair.cm = source[k][j][i][c];
                            target[k][j][i][c] = PicurvCoMomentStateCovariance(&pair);
                        }
                    } else {
                        pair.cm = source[k][j][i][member];
                        target[k][j][i][0] = PicurvCoMomentStateCovariance(&pair);
                    }
                }
            }
        }
 
        PetscCall(DMDAVecRestoreArrayDOF(destination_dm, destination, &target));
        PetscCall(DMDAVecRestoreArrayDOFRead(source_dm, source_vec, &source));
        PetscCall(DMDAVecRestoreArrayRead(user->da, storage->count, &count_arr));
        PetscCall(DMDAVecRestoreArrayRead(user->da, storage->weight_sq, &weight_sq_arr));
        PetscCall(DMDAVecRestoreArrayRead(user->da, storage->weight, &weight_arr));
 
        /* Applied once, here, so the VTK field and the convergence-history CSV - which
         * both read this staging vector - cannot end up in different unit systems. */
        if (user->simCtx->pps && user->simCtx->pps->dimensionalize) {
            if (kind != DERIVED_FLUX && descriptor) {
                PetscReal base = 1.0;
 
                PetscCall(PicurvFieldReferenceScale(user->simCtx, descriptor->canonical_name,
                                                    &base, NULL, 0));
                dimensional_factor = PetscPowRealInt(base, kDerivedKindScaleExponent[kind]);
            }
            if (PetscAbsReal(dimensional_factor - 1.0) > PETSC_MACHINE_EPSILON) {
                PetscCall(VecScale(destination, dimensional_factor));
                LOG_ALLOW(GLOBAL, LOG_DEBUG, "Scaled derived '%s' by %.4e.\n",
                          kDerivedKindName[kind], (double)dimensional_factor);
            }
        }
    }
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PicurvWindowSpatialMean()

PetscErrorCode PicurvWindowSpatialMean	(	UserCtx *	user,
		const PicurvWindowDefinition *	definition,
		const PicurvWindowStorage *	storage,
		Vec	field,
		PetscReal *	mean
	)

Implementation of PicurvWindowSpatialMean().

Reports the spatial mean of a derived field over the points a window sampled.

See also

PicurvWindowSpatialMean()

Definition at line 793 of file statistics_accumulator.c.

{
    SpatialTargetPlan plan;
    /* Every direction collapses: a window mean is one number for the domain. */
    const PetscBool   whole_domain[3] = {PETSC_TRUE, PETSC_TRUE, PETSC_TRUE};
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && definition != NULL && storage != NULL && field != NULL && mean != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context, definition, storage, field, and output are required.");
    *mean = 0.0;
    if (definition->field_count == 0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    PetscCall(SpatialTargetPlanCreate(user, (FieldId)definition->fields[0].field_id,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    /* No denominator means each admitted point counts once, so the result is the mean
     * of the field over those points rather than a volume-weighted average. The
     * accumulated weight is passed as the inclusion mask: a point the window never
     * sampled holds a zero that means "never measured", and averaging it in would
     * scale the answer down by the fraction of the domain the window never covered. */
    PetscCall(PicurvSpatialRatioAverage(user, &plan, field, NULL, storage->weight,
                                        whole_domain, PETSC_COMM_WORLD, NULL, mean));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PicurvWindowValidFractionRange()

PetscErrorCode PicurvWindowValidFractionRange	(	UserCtx *	user,
		const PicurvWindowDefinition *	definition,
		const PicurvWindowStorage *	storage,
		PetscInt	sample_count,
		PetscReal *	minimum,
		PetscReal *	maximum
	)

Implementation of PicurvWindowValidFractionRange().

Reports the range of per-point valid fraction across a window's domain.

See also

PicurvWindowValidFractionRange()

Definition at line 828 of file statistics_accumulator.c.

{
    SpatialTargetPlan plan;
    PetscReal ***nvert = NULL;
    PetscReal ***count_arr = NULL;
    PetscReal local_min = PETSC_MAX_REAL;
    PetscReal local_max = 0.0;
    PetscReal reduced_min = 0.0, reduced_max = 0.0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && definition != NULL && storage != NULL &&
               minimum != NULL && maximum != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Context, definition, storage, and outputs are required.");
    *minimum = 1.0;
    *maximum = 0.0;
    if (definition->field_count == 0 || sample_count <= 0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    PetscCall(SpatialTargetPlanCreate(user, (FieldId)definition->fields[0].field_id,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
    PetscCall(DMDAVecGetArrayRead(user->da, user->Nvert, &nvert));
    PetscCall(DMDAVecGetArrayRead(user->da, storage->count, &count_arr));
    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) {
                /* The mask is evaluated at the current state, but a point excluded
                 * now may have contributed earlier, so its stored count is what
                 * decides its fraction rather than its present eligibility. */
                const PetscReal fraction = count_arr[k][j][i] / (PetscReal)sample_count;
 
                local_min = PetscMin(local_min, fraction);
                local_max = PetscMax(local_max, fraction);
            }
        }
    }
    PetscCall(DMDAVecRestoreArrayRead(user->da, storage->count, &count_arr));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->Nvert, &nvert));
 
    /* A rank owning no targeted point must not drag the minimum down, so its
     * sentinel is neutral under the reduction rather than zero. */
    PetscCallMPI(MPI_Allreduce(&local_min, &reduced_min, 1, MPIU_REAL, MPIU_MIN, PETSC_COMM_WORLD));
    PetscCallMPI(MPI_Allreduce(&local_max, &reduced_max, 1, MPIU_REAL, MPIU_MAX, PETSC_COMM_WORLD));
    *minimum = (reduced_min == PETSC_MAX_REAL) ? 1.0 : reduced_min;
    *maximum = reduced_max;
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PicurvWindowAccumulate()

PetscErrorCode PicurvWindowAccumulate	(	UserCtx *	user,
		const PicurvWindowDefinition *	definition,
		PicurvWindowStorage *	storage,
		PetscReal	weight
	)

Implementation of PicurvWindowAccumulate().

Applies one accepted completed state to a window's accumulators.

See also

PicurvWindowAccumulate()

Definition at line 884 of file statistics_accumulator.c.

{
    SpatialTargetPlan plan;
    PetscReal ***nvert = NULL;
    PetscReal ***count_arr = NULL, ***weight_arr = NULL, ***weight_sq_arr = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(user != NULL && definition != NULL && storage != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Context, definition, and storage are required.");
    PetscCheck(weight > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Accepted states carry a positive weight, got %g.", (double)weight);
    if (definition->field_count == 0) { PROFILE_FUNCTION_END; PetscFunctionReturn(0); }
 
    /* Every requested field shares the pointwise cell-centered domain, so the plan
     * is resolved once rather than per field. */
    PetscCall(SpatialTargetPlanCreate(user, (FieldId)definition->fields[0].field_id,
                                      PICURV_STATISTICS_MASK_FLUID, &plan));
 
    PetscCall(DMDAVecGetArrayRead(user->da, user->Nvert, &nvert));
    PetscCall(DMDAVecGetArray(user->da, storage->count, &count_arr));
    PetscCall(DMDAVecGetArray(user->da, storage->weight, &weight_arr));
    PetscCall(DMDAVecGetArray(user->da, storage->weight_sq, &weight_sq_arr));
 
    /* Pass one advances per-point occupancy, so every field pass afterwards can
     * recover the pre-state weight uniformly as (stored weight - this weight).
     * Occupancy is per point and per accepted state, never per field. */
    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) {
                if (!SpatialTargetPlanMaskAllows(&plan, nvert[k][j][i])) continue;
                count_arr[k][j][i]     += 1.0;
                weight_arr[k][j][i]    += weight;
                weight_sq_arr[k][j][i] += weight * weight;
            }
        }
    }
 
    /* Pass two updates every cross-field co-moment. It runs before any mean is
     * written back, because a co-moment needs the pre-update mean of *both*
     * members, and pass three overwrites them field by field. */
    for (PetscInt pair = 0; pair < definition->covariance_count; ++pair) {
        FieldView view_a, view_b;
        PetscScalar ****src_a = NULL, ****mean_a = NULL;
        PetscScalar ****src_b = NULL, ****mean_b = NULL;
        PetscScalar ****co_moment = NULL;
        DM pair_dm = NULL;
        PetscInt slot_a = 0, slot_b = 0, dof_a = 0, dof_b = 0, components = 0;
 
        PetscCall(FindFieldSlot(definition, definition->covariances[pair].first, &slot_a));
        PetscCall(FindFieldSlot(definition, definition->covariances[pair].second, &slot_b));
        PetscCall(FieldGetView(user, (FieldId)definition->covariances[pair].first, &view_a));
        PetscCall(FieldGetView(user, (FieldId)definition->covariances[pair].second, &view_b));
        dof_a = view_a.descriptor->dof;
        dof_b = view_b.descriptor->dof;
        PetscCall(PicurvCovarianceComponentCount(dof_a, dof_b, &components));
        PetscCall(PicurvStatisticsComponentDM(user, components, &pair_dm));
 
        /* A component-indexed view serves every degree of freedom, so the loop
         * below needs no scalar-versus-vector branching. */
        PetscCall(DMDAVecGetArrayDOFRead(view_a.dm, view_a.global_vec, &src_a));
        PetscCall(DMDAVecGetArrayDOFRead(view_a.dm, storage->mean[slot_a], &mean_a));
        PetscCall(DMDAVecGetArrayDOFRead(view_b.dm, view_b.global_vec, &src_b));
        PetscCall(DMDAVecGetArrayDOFRead(view_b.dm, storage->mean[slot_b], &mean_b));
        PetscCall(DMDAVecGetArrayDOF(pair_dm, storage->cm[pair], &co_moment));
 
        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) {
                    if (!SpatialTargetPlanMaskAllows(&plan, nvert[k][j][i])) continue;
                    for (PetscInt c = 0; c < components; ++c) {
                        /* A scalar member contributes its single value against every
                         * component of a vector member, which is what makes a
                         * vector-scalar covariance a three-component object. */
                        const PetscInt component_a = (dof_a == 3) ? c : 0;
                        const PetscInt component_b = (dof_b == 3) ? c : 0;
                        PicurvCoMomentState state;
 
                        state.count = count_arr[k][j][i] - 1.0;
                        state.weight = weight_arr[k][j][i] - weight;
                        state.weight_sq = weight_sq_arr[k][j][i] - weight * weight;
                        state.mean_x = mean_a[k][j][i][component_a];
                        state.mean_y = mean_b[k][j][i][component_b];
                        state.cm = co_moment[k][j][i][c];
                        PetscCall(PicurvCoMomentStateUpdate(&state,
                                                            src_a[k][j][i][component_a],
                                                            src_b[k][j][i][component_b], weight));
                        co_moment[k][j][i][c] = state.cm;
                    }
                }
            }
        }
 
        PetscCall(DMDAVecRestoreArrayDOF(pair_dm, storage->cm[pair], &co_moment));
        PetscCall(DMDAVecRestoreArrayDOFRead(view_b.dm, storage->mean[slot_b], &mean_b));
        PetscCall(DMDAVecRestoreArrayDOFRead(view_b.dm, view_b.global_vec, &src_b));
        PetscCall(DMDAVecRestoreArrayDOFRead(view_a.dm, storage->mean[slot_a], &mean_a));
        PetscCall(DMDAVecRestoreArrayDOFRead(view_a.dm, view_a.global_vec, &src_a));
    }
 
    /* Pass three updates each requested field's mean and, when asked, its centered
     * self-product. A three-vector's self-product is the six symmetric co-moments
     * between component pairs, not three per-component variances, so the co-moment
     * kernel drives every product component including the diagonal. */
    for (PetscInt field_index = 0; field_index < definition->field_count; ++field_index) {
        FieldView view;
        PetscScalar ****src = NULL, ****mean = NULL, ****product = NULL;
        DM product_dm = NULL;
        PetscInt dof = 0, components = 0;
        const PetscBool second = definition->fields[field_index].want_second;
 
        PetscCall(FieldGetView(user, (FieldId)definition->fields[field_index].field_id, &view));
        dof = view.descriptor->dof;
        PetscCall(DMDAVecGetArrayDOFRead(view.dm, view.global_vec, &src));
        PetscCall(DMDAVecGetArrayDOF(view.dm, storage->mean[field_index], &mean));
        if (second) {
            PetscCall(PicurvProductComponentCount(dof, &components));
            PetscCall(PicurvStatisticsComponentDM(user, components, &product_dm));
            PetscCall(DMDAVecGetArrayDOF(product_dm, storage->m2[field_index], &product));
        }
 
        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) {
                    PetscReal prior_weight = 0.0;
                    PetscReal prior_weight_sq = 0.0;
                    PetscReal prior_count = 0.0;
 
                    if (!SpatialTargetPlanMaskAllows(&plan, nvert[k][j][i])) continue;
 
                    prior_weight    = weight_arr[k][j][i] - weight;
                    prior_weight_sq = weight_sq_arr[k][j][i] - weight * weight;
                    prior_count     = count_arr[k][j][i] - 1.0;
 
                    /* Products are updated before the means are written back,
                     * because the co-moment kernel needs the pre-update means. */
                    for (PetscInt c = 0; c < components; ++c) {
                        const PetscInt a = (dof == 1) ? 0 : kProductFirst[c];
                        const PetscInt b = (dof == 1) ? 0 : kProductSecond[c];
                        PicurvCoMomentState state;
 
                        state.count = prior_count;
                        state.weight = prior_weight;
                        state.weight_sq = prior_weight_sq;
                        state.mean_x = mean[k][j][i][a];
                        state.mean_y = mean[k][j][i][b];
                        state.cm = product[k][j][i][c];
                        PetscCall(PicurvCoMomentStateUpdate(&state, src[k][j][i][a],
                                                            src[k][j][i][b], weight));
                        product[k][j][i][c] = state.cm;
                    }
 
                    for (PetscInt c = 0; c < dof; ++c) {
                        PicurvMomentState moment;
 
                        moment.count = prior_count;
                        moment.weight = prior_weight;
                        moment.weight_sq = prior_weight_sq;
                        moment.mean = mean[k][j][i][c];
                        moment.m2 = 0.0;
                        PetscCall(PicurvMomentStateUpdate(&moment, src[k][j][i][c], weight));
                        mean[k][j][i][c] = moment.mean;
                    }
                }
            }
        }
 
        if (second) PetscCall(DMDAVecRestoreArrayDOF(product_dm, storage->m2[field_index], &product));
        PetscCall(DMDAVecRestoreArrayDOF(view.dm, storage->mean[field_index], &mean));
        PetscCall(DMDAVecRestoreArrayDOFRead(view.dm, view.global_vec, &src));
    }
 
    PetscCall(DMDAVecRestoreArray(user->da, storage->weight_sq, &weight_sq_arr));
    PetscCall(DMDAVecRestoreArray(user->da, storage->weight, &weight_arr));
    PetscCall(DMDAVecRestoreArray(user->da, storage->count, &count_arr));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->Nvert, &nvert));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

kProductFirst

const PetscInt kProductFirst[6] = {0, 0, 0, 1, 1, 2}

static

Upper-triangular row-major component pairs for a three-vector self-product.

This pair table is the single definition of the symmetric component order. The diagonal positions and the two-axis component labels are both derived from it rather than restated, because accumulation reads the table while derivation reads what follows from it: a hand-written copy that drifted would mislabel every Reynolds stress and take the root of the wrong component, with nothing failing.

Definition at line 30 of file statistics_accumulator.c.

{0, 0, 0, 1, 1, 2};

kProductSecond

const PetscInt kProductSecond[6] = {0, 1, 2, 1, 2, 2}

static

Definition at line 31 of file statistics_accumulator.c.

{0, 1, 2, 1, 2, 2};

kAxisName

const char* const kAxisName[3] = {"x", "y", "z"}

static

Axis labels indexing the pair table above.

Definition at line 48 of file statistics_accumulator.c.

{"x", "y", "z"};

kDerivedKindName

const char* const kDerivedKindName[DERIVED_KIND_COUNT]

static

Initial value:

= {
    "mean", "reynolds_stress", "rms", "tke", "flux"
}

Recipe spellings of the output kinds.

Definition at line 390 of file statistics_accumulator.c.

                                                                {
    "mean", "reynolds_stress", "rms", "tke", "flux"
};

kDerivedKindScaleExponent

const PetscInt kDerivedKindScaleExponent[DERIVED_KIND_COUNT]

static

Initial value:

= {
    1,  
    2,  
    1,  
    2,  
    0   
}

Power of the source field's reference scale each derived kind carries.

A first moment and a standard deviation have the field's own units; a covariance, its trace, and a co-moment flux are quadratic in it. This is what the per-field scaling table alone cannot express, and why derived statistics were previously left non-dimensional rather than scaled by a velocity that would have been wrong for three of the five kinds.

Definition at line 403 of file statistics_accumulator.c.

                                                                      {
    1,  /* mean             - first moment, the field's own units */
    2,  /* reynolds_stress  - covariance of one field with itself */
    1,  /* rms              - standard deviation */
    2,  /* tke              - half the trace of that covariance */
    0   /* flux             - a co-moment of two possibly different fields, so its
                                factor is their product and cannot be an exponent of
                                one scale; DERIVED_FLUX resolves it directly. */
};