The expression language shared by every configured field, and the keyed hashes its random draws and other reproducible draws are built from. More...

#include <stdint.h>
#include <petscsys.h>
#include "variables.h"

Include dependency graph for expression.h:

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

Go to the source code of this file.

Data Structures
struct	PicurvExpressionDrawKey
	What a random draw is keyed on, besides the particle and its stream. More...

struct	PicurvVelocitySource
	A prescribed velocity, for any target: one compiled expression per value, the frame they are in, and the bulk velocity it is scaled to, if any. More...

struct	PicurvSourceFaceSlots
	What a source reads at each boundary-face slot this rank owns, gathered once. More...

Typedefs
typedef struct PicurvExpression	PicurvExpression
	A compiled expression: postfix code for a small stack evaluator.

Enumerations
enum	PicurvDrawSpace { PICURV_DRAW_SPACE_FIELD = 0 , PICURV_DRAW_SPACE_NAMED_STREAM = 1 , PICURV_DRAW_SPACE_BROWNIAN = 2 , PICURV_DRAW_SPACE_PLACEMENT = 3 , PICURV_DRAW_SPACE_INLET = 4 }
	Disjoint key spaces of the keyed particle draws, so draws made for different purposes never coincide for one particle. More...

enum	PicurvSourceVariable { PICURV_SOURCE_VAR_X = 0 , PICURV_SOURCE_VAR_Y , PICURV_SOURCE_VAR_Z , PICURV_SOURCE_VAR_T , PICURV_SOURCE_VAR_XI , PICURV_SOURCE_VAR_ETA , PICURV_SOURCE_VAR_ZETA , PICURV_SOURCE_VAR_L_XI , PICURV_SOURCE_VAR_L_ETA , PICURV_SOURCE_VAR_L_ZETA , PICURV_SOURCE_VAR_CELL_COUNT , PICURV_SOURCE_VAR_R = PICURV_SOURCE_VAR_CELL_COUNT , PICURV_SOURCE_VAR_COUNT }
	Variables a prescribed velocity source may use, in the order of `values[]`. More...

enum	PicurvSourceFrame { PICURV_SOURCE_FRAME_CARTESIAN = 0 , PICURV_SOURCE_FRAME_CURVILINEAR = 1 , PICURV_SOURCE_FRAME_NORMAL = 2 , PICURV_SOURCE_FRAME_CONTRAVARIANT = 3 }
	How a source's values are read. More...

Functions
uint64_t	PicurvMix64 (uint64_t z)
	splitmix64 finalizer: a bijective mix of 64 bits.

PetscReal	PicurvUnitFromHash (uint64_t hash)
	Uniform on [0, 1) from the top 53 bits of a hash.

uint64_t	PicurvParticleDrawHash (PetscInt64 seed, PicurvDrawSpace space, uint64_t stream, PetscInt64 pid, PetscInt64 salt)
	Hash that keys a particle draw: a pure function of its arguments, so a draw does not depend on the rank a particle lives on, the order particles are visited, or what was drawn before it.

PetscErrorCode	PicurvExpressionCompile (const char text, const char const names, PetscInt name_count, PetscBool allow_random, PicurvExpression *expression)
	Compile one expression against a list of variable names.

PetscErrorCode	PicurvExpressionEvaluate (const PicurvExpression expression, const PetscReal values, const PicurvExpressionDrawKey key, PetscReal result)
	Evaluate a compiled expression.

PetscErrorCode	PicurvExpressionDestroy (PicurvExpression **expression)
	Free a compiled expression.

PetscErrorCode	PicurvExpressionUsesVariable (const PicurvExpression expression, PetscInt index, PetscBool uses)
	Whether a compiled expression reads one of its variables.

PetscErrorCode	PicurvVelocitySourceCreate (const char prefix, PetscBool on_face, PetscBool has_data, PicurvVelocitySource source)
	Read and compile a velocity source from the options database.

PetscErrorCode	PicurvVelocitySourceBulkFactor (const PicurvVelocitySource source, PetscReal physical_time, const PetscReal local[3], const char where, PetscReal *factor)
	Factor on a source's expressions that makes a target's whole mean equal its bulk velocity.

static void	PicurvSourceSetVariables (PetscReal *values, Cmpnts point, Cmpnts fraction, Cmpnts length, PetscReal L_ref)
	Fill a source's position and grid-line variables at one point.

PetscErrorCode	PicurvSourceFaceSlotsCreate (UserCtx user, BCFace face_id, Vec fraction, Vec length, PicurvSourceFaceSlots slots)
	Gather each owned slot of a boundary face: face-centre position, the grid-line fractions and lengths of the interior cell behind it, and its directions.

PetscErrorCode	PicurvSourceFaceSlotsDestroy (PicurvSourceFaceSlots *slots)
	Release what PicurvSourceFaceSlotsCreate() allocated.

PetscErrorCode	PicurvVelocitySourceEvaluate (const PicurvVelocitySource source, const PetscReal values, const Cmpnts edges[3], Cmpnts *value)
	Evaluate a source at one point.

PetscErrorCode	PicurvSourceSetCellPressure (UserCtx user, const PicurvExpression pressure, Vec fraction, Vec length, PetscReal physical_time)
	Set a block's pressure at every interior cell centre from an expression in Pa.

PetscErrorCode	PicurvSourceCellBulkFactor (UserCtx user, const PicurvVelocitySource source, Vec fraction, Vec length, PetscReal physical_time, PetscReal *factor)
	Factor on a source's expressions that makes the volume mean of the whole field, along the source's own direction at each cell, equal its bulk velocity.

PetscErrorCode	PicurvSourceAddFluxes (UserCtx user, const PicurvVelocitySource source, Vec fraction, Vec length, PetscReal physical_time, PetscReal factor, PetscBool boundary_faces)
	Add a source's flux to every owned face of Ucont, evaluating it at each face centre.

PetscErrorCode	PicurvVelocitySourceDestroy (PicurvVelocitySource *source)
	Free a source's compiled expressions and bulk velocity.

Variables
const char *const	kPicurvSourceVariableNames [PICURV_SOURCE_VAR_COUNT]
	Names of PicurvSourceVariable, in order.

Detailed Description

The expression language shared by every configured field, and the keyed hashes its random draws and other reproducible draws are built from.

Callers compile an expression against their own list of variable names and supply the values in that order at evaluation. generators/ic.gen validates the same language; tests/fixtures/expression_conformance.txt holds cases both implementations must agree on.

Definition in file expression.h.

Data Structure Documentation

◆ PicurvExpressionDrawKey

struct PicurvExpressionDrawKey

What a random draw is keyed on, besides the particle and its stream.

Definition at line 67 of file expression.h.

Data Fields
PetscInt64	seed	Run seed (`-particle_random_seed`).
PetscInt64	pid	Particle identity; draws are pure functions of it.
PetscInt64	salt	Event that produced the particle: 0 for the t=0 population.
PetscInt64	field	Field whose expression draws, the default stream's namespace.

◆ PicurvVelocitySource

struct PicurvVelocitySource

A prescribed velocity, for any target: one compiled expression per value, the frame they are in, and the bulk velocity it is scaled to, if any.

A target differs only in the points it evaluates the source at and the weight it takes a mean with.

Definition at line 155 of file expression.h.

Collaboration diagram for PicurvVelocitySource:

[legend]

Data Fields
PicurvSourceFrame	frame
PetscInt	count	1 for the normal frame, otherwise 3; 0 for staged data alone.
PicurvExpression *	component[3]
PetscBool	time_dependent	Some component reads `t`.
PicurvExpression *	bulk	Target mean speed [m/s] in `t`, or NULL.
PetscBool	bulk_time_dependent
PetscBool	bulk_is_flow_rate	`bulk` is a volumetric flow rate [m^3/s] through a face instead.

◆ PicurvSourceFaceSlots

struct PicurvSourceFaceSlots

What a source reads at each boundary-face slot this rank owns, gathered once.

Definition at line 221 of file expression.h.

Collaboration diagram for PicurvSourceFaceSlots:

[legend]

Data Fields
PetscInt	normal	Logical axis of the face normal.
PetscInt	a
PetscInt	b	The two tangential axes.
PetscInt	layer	Index along the normal of the face's metric and flux.
PetscInt	lo[2]	First tangential indices of the owned slots.
PetscInt	n[2]	Owned slots along the two tangential axes; 0 off the face.
PetscReal(*	values)[PICURV_SOURCE_VAR_COUNT]	Variables per slot; `t` is the caller's.
Cmpnts(*	edges)[3]	Directions per slot, for the curvilinear frame.

Typedef Documentation

◆ PicurvExpression

typedef struct PicurvExpression PicurvExpression

A compiled expression: postfix code for a small stack evaluator.

Definition at line 64 of file expression.h.

Enumeration Type Documentation

◆ PicurvDrawSpace

enum PicurvDrawSpace

Disjoint key spaces of the keyed particle draws, so draws made for different purposes never coincide for one particle.

Enumerator
PICURV_DRAW_SPACE_FIELD	An expression's default, per-field stream.
PICURV_DRAW_SPACE_NAMED_STREAM	`uniform(k)`/`normal(k)`: a stream shared across fields.
PICURV_DRAW_SPACE_BROWNIAN	Brownian increments; the salt is the step.
PICURV_DRAW_SPACE_PLACEMENT	Seeding points; the salt is the rejection attempt.
PICURV_DRAW_SPACE_INLET	Injection at an inlet: stream 0 counts per face cell, stream 1 distance in.

Definition at line 27 of file expression.h.

             {
    PICURV_DRAW_SPACE_FIELD = 0,        /**< An expression's default, per-field stream. */
    PICURV_DRAW_SPACE_NAMED_STREAM = 1, /**< `uniform(k)`/`normal(k)`: a stream shared across fields. */
    PICURV_DRAW_SPACE_BROWNIAN = 2,     /**< Brownian increments; the salt is the step. */
    PICURV_DRAW_SPACE_PLACEMENT = 3,    /**< Seeding points; the salt is the rejection attempt. */
    PICURV_DRAW_SPACE_INLET = 4         /**< Injection at an inlet: stream 0 counts per face cell, stream 1 distance in. */
} PicurvDrawSpace;

◆ PicurvSourceVariable

enum PicurvSourceVariable

Variables a prescribed velocity source may use, in the order of values[].

Enumerator
PICURV_SOURCE_VAR_X	Physical x of the point [m].
PICURV_SOURCE_VAR_Y	Physical y [m].
PICURV_SOURCE_VAR_Z	Physical z [m].
PICURV_SOURCE_VAR_T	Physical time [s].
PICURV_SOURCE_VAR_XI	Arc-length fraction along the xi grid line, 0 at the - face.
PICURV_SOURCE_VAR_ETA	Arc-length fraction along the eta grid line.
PICURV_SOURCE_VAR_ZETA	Arc-length fraction along the zeta grid line.
PICURV_SOURCE_VAR_L_XI	Physical length of that xi grid line [m].
PICURV_SOURCE_VAR_L_ETA	Physical length of that eta grid line [m].
PICURV_SOURCE_VAR_L_ZETA	Physical length of that zeta grid line [m].
PICURV_SOURCE_VAR_CELL_COUNT	Variables a source reads on every target.
PICURV_SOURCE_VAR_R	On a face, distance from the face's centre [m].
PICURV_SOURCE_VAR_COUNT

Definition at line 126 of file expression.h.

             {
    PICURV_SOURCE_VAR_X = 0,   /**< Physical x of the point [m]. */
    PICURV_SOURCE_VAR_Y,       /**< Physical y [m]. */
    PICURV_SOURCE_VAR_Z,       /**< Physical z [m]. */
    PICURV_SOURCE_VAR_T,       /**< Physical time [s]. */
    PICURV_SOURCE_VAR_XI,      /**< Arc-length fraction along the xi grid line, 0 at the - face. */
    PICURV_SOURCE_VAR_ETA,     /**< Arc-length fraction along the eta grid line. */
    PICURV_SOURCE_VAR_ZETA,    /**< Arc-length fraction along the zeta grid line. */
    PICURV_SOURCE_VAR_L_XI,    /**< Physical length of that xi grid line [m]. */
    PICURV_SOURCE_VAR_L_ETA,   /**< Physical length of that eta grid line [m]. */
    PICURV_SOURCE_VAR_L_ZETA,  /**< Physical length of that zeta grid line [m]. */
    PICURV_SOURCE_VAR_CELL_COUNT, /**< Variables a source reads on every target. */
    PICURV_SOURCE_VAR_R = PICURV_SOURCE_VAR_CELL_COUNT, /**< On a face, distance from the face's centre [m]. */
    PICURV_SOURCE_VAR_COUNT
} PicurvSourceVariable;

◆ PicurvSourceFrame

enum PicurvSourceFrame

How a source's values are read.

Enumerator
PICURV_SOURCE_FRAME_CARTESIAN	Three values: x, y, z velocity components.
PICURV_SOURCE_FRAME_CURVILINEAR	Three values along the unit xi, eta, zeta grid-line directions.
PICURV_SOURCE_FRAME_NORMAL	One value: the speed into the domain through a face.
PICURV_SOURCE_FRAME_CONTRAVARIANT	Three values: the volume flux through the xi, eta, zeta faces; cells only.

Definition at line 143 of file expression.h.

             {
    PICURV_SOURCE_FRAME_CARTESIAN = 0,   /**< Three values: x, y, z velocity components. */
    PICURV_SOURCE_FRAME_CURVILINEAR = 1, /**< Three values along the unit xi, eta, zeta grid-line directions. */
    PICURV_SOURCE_FRAME_NORMAL = 2,      /**< One value: the speed into the domain through a face. */
    PICURV_SOURCE_FRAME_CONTRAVARIANT = 3 /**< Three values: the volume flux through the xi, eta, zeta faces; cells only. */
} PicurvSourceFrame;

Function Documentation

◆ PicurvMix64()

uint64_t PicurvMix64 ( uint64_t z )

splitmix64 finalizer: a bijective mix of 64 bits.

Parameters

[in] z Bits to mix.

Returns: The mixed bits.

splitmix64 finalizer: a bijective mix of 64 bits.

See also: PicurvMix64()

Definition at line 587 of file expression.c.

{
    z += 0x9E3779B97F4A7C15ULL;
    z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
    z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
    return z ^ (z >> 31);
}

Here is the caller graph for this function:

◆ PicurvUnitFromHash()

PetscReal PicurvUnitFromHash ( uint64_t hash )

Uniform on [0, 1) from the top 53 bits of a hash.

Parameters

[in] hash Mixed bits.

Returns: A value in [0, 1).

Uniform on [0, 1) from the top 53 bits of a hash.

See also: PicurvUnitFromHash()

Definition at line 599 of file expression.c.

{
    return (PetscReal)(hash >> 11) * (1.0 / 9007199254740992.0);
}

Here is the caller graph for this function:

◆ PicurvParticleDrawHash()

uint64_t PicurvParticleDrawHash	(	PetscInt64	seed,
		PicurvDrawSpace	space,
		uint64_t	stream,
		PetscInt64	pid,
		PetscInt64	salt
	)

Hash that keys a particle draw: a pure function of its arguments, so a draw does not depend on the rank a particle lives on, the order particles are visited, or what was drawn before it.

Parameters

[in]	seed	Run seed (`-particle_random_seed`).
[in]	space	Purpose of the draw.
[in]	stream	Stream within the space.
[in]	pid	Particle identity.
[in]	salt	Event that distinguishes repeated draws for one particle.

Returns: 64 mixed bits; derive further independent values with PicurvMix64().

Hash that keys a particle draw: a pure function of its arguments, so a draw does not depend on the rank a particle lives on, the order particles are visited, or what was drawn before it.

See also: PicurvParticleDrawHash()

Definition at line 608 of file expression.c.

{
    uint64_t hash = PicurvMix64((uint64_t)seed);
    hash = PicurvMix64(hash ^ (uint64_t)space);
    hash = PicurvMix64(hash ^ stream);
    hash = PicurvMix64(hash ^ (uint64_t)pid);
    return PicurvMix64(hash ^ (uint64_t)salt);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvExpressionCompile()

PetscErrorCode PicurvExpressionCompile	(	const char *	text,
		const char const	names,
		PetscInt	name_count,
		PetscBool	allow_random,
		PicurvExpression **	expression
	)

Compile one expression against a list of variable names.

The grammar is Python's expression syntax restricted to numbers, the given names, pi, arithmetic (+ - * / % **), comparisons (chained ones included), and/or/not, and the functions abs atan2 cos cosh erf exp floor log maximum minimum sin sinh sqrt tan tanh where, plus uniform and normal when allow_random is true. Comparisons and logical operators yield 1 or 0; % takes the sign of the divisor; ** binds tighter than a unary minus on its left and groups to the right.

Parameters

[in]	text	Expression text.
[in]	names	Variable names, in the order evaluation supplies their values.
[in]	name_count	Number of names.
[in]	allow_random	Whether `uniform` and `normal` may be called.
[out]	expression	Compiled expression, owned by the caller.

Returns: Zero on success; PETSC_ERR_ARG_WRONG for a syntax error or an unknown name or function, with the offending position in the message.

Compile one expression against a list of variable names.

See also: PicurvExpressionCompile()

Definition at line 549 of file expression.c.

{
    ExpressionParser parser;
    PetscInt         root, depth = 0;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCheck(text && expression, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Expression text and output are required.");
    PetscCall(PetscMemzero(&parser, sizeof(parser)));
    parser.text = text;
    parser.names = names;
    parser.name_count = name_count;
    parser.allow_random = allow_random;
    parser.capacity = 4 * (PetscInt)strlen(text) + 8;
    PetscCall(PetscMalloc1(parser.capacity, &parser.nodes));
 
    root = ParseOr(&parser);
    SkipSpace(&parser);
    if (!parser.error[0] && parser.text[parser.pos]) ParserFail(&parser, "unexpected text");
    if (parser.error[0]) {
        PetscCall(PetscFree(parser.nodes));
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid expression: %s.", parser.error);
    }
 
    PetscCall(PetscNew(expression));
    (*expression)->length = 0;
    PetscCall(PetscMalloc1(CountCode(&parser, root), &(*expression)->code));
    EmitNode(&parser, root, *expression, &depth);
    PetscCall(PetscFree(parser.nodes));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvExpressionEvaluate()

PetscErrorCode PicurvExpressionEvaluate	(	const PicurvExpression *	expression,
		const PetscReal *	values,
		const PicurvExpressionDrawKey *	key,
		PetscReal *	result
	)

Evaluate a compiled expression.

Parameters

[in]	expression	Compiled expression.
[in]	values	One value per compiled name, in compile order.
[in]	key	Draw key, required when the expression draws; may be NULL otherwise.
[out]	result	Value of the expression; may be non-finite, which the caller judges.

Returns: Zero on success.

Evaluate a compiled expression.

See also: PicurvExpressionEvaluate()

Definition at line 695 of file expression.c.

{
    PetscReal stack[expression->max_depth + 1];
    PetscReal index[EXPRESSION_MAX_SERIES_DEPTH];
    PetscInt  top = 0;
 
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    PetscCheck(!expression->draws || key, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "An expression that draws random numbers needs a draw key.");
    for (PetscInt i = 0; i < expression->length; ++i) {
        const ExpressionInstruction *instruction = &expression->code[i];
        PetscReal a, b, c;
 
        switch (instruction->op) {
        case OP_CONST: stack[top++] = instruction->value; break;
        case OP_VAR: stack[top++] = values[instruction->arg]; break;
        case OP_UNIFORM:
        case OP_NORMAL: stack[top++] = Draw(key, instruction); break;
        case OP_NEG: stack[top - 1] = -stack[top - 1]; break;
        case OP_POS: break;
        case OP_NOT: stack[top - 1] = (stack[top - 1] == 0.0) ? 1.0 : 0.0; break;
        case OP_ABS: stack[top - 1] = PetscAbsReal(stack[top - 1]); break;
        case OP_COS: stack[top - 1] = PetscCosReal(stack[top - 1]); break;
        case OP_EXP: stack[top - 1] = PetscExpReal(stack[top - 1]); break;
        case OP_SIN: stack[top - 1] = PetscSinReal(stack[top - 1]); break;
        case OP_SQRT: stack[top - 1] = PetscSqrtReal(stack[top - 1]); break;
        case OP_TAN: stack[top - 1] = PetscTanReal(stack[top - 1]); break;
        case OP_ERF: stack[top - 1] = (PetscReal)erf((double)stack[top - 1]); break;
        case OP_TANH: stack[top - 1] = PetscTanhReal(stack[top - 1]); break;
        case OP_COSH: stack[top - 1] = PetscCoshReal(stack[top - 1]); break;
        case OP_SINH: stack[top - 1] = PetscSinhReal(stack[top - 1]); break;
        case OP_LOG: stack[top - 1] = PetscLogReal(stack[top - 1]); break;
        case OP_FLOOR: stack[top - 1] = PetscFloorReal(stack[top - 1]); break;
        case OP_BER: stack[top - 1] = Kelvin(stack[top - 1], PETSC_FALSE); break;
        case OP_BEI: stack[top - 1] = Kelvin(stack[top - 1], PETSC_TRUE); break;
        case OP_LOCAL: stack[top++] = index[instruction->arg]; break;
        case OP_SUM_START:
            stack[top++] = 0.0;
            index[instruction->arg] = instruction->value;
            /* An empty range skips the body and its ACC. */
            if (instruction->value > expression->code[i + 1].value) i += 2 + expression->code[i + 1].arg;
            else i += 1;
            break;
        case OP_SUM_ACC: {
            const PetscInt  body = (PetscInt)instruction->value;
            const PetscReal upper = expression->code[i - body - 1].value;
            const PetscReal term = stack[--top];
            stack[top - 1] += term;
            index[instruction->arg] += 1.0;
            if (index[instruction->arg] <= upper) i -= body + 1;
        } break;
        case OP_WHERE:
            c = stack[--top]; b = stack[--top]; a = stack[top - 1];
            stack[top - 1] = (a != 0.0) ? b : c;
            break;
        default:
            b = stack[--top]; a = stack[top - 1];
            switch (instruction->op) {
            case OP_ADD: a = a + b; break;
            case OP_SUB: a = a - b; break;
            case OP_MUL: a = a * b; break;
            case OP_DIV: a = a / b; break;
            /* Python's modulus takes the sign of the divisor. */
            case OP_MOD: a = (b == 0.0) ? NAN : a - b * PetscFloorReal(a / b); break;
            case OP_POW: a = PetscPowReal(a, b); break;
            case OP_EQ: a = (a == b) ? 1.0 : 0.0; break;
            case OP_NE: a = (a != b) ? 1.0 : 0.0; break;
            case OP_LT: a = (a < b) ? 1.0 : 0.0; break;
            case OP_LE: a = (a <= b) ? 1.0 : 0.0; break;
            case OP_GT: a = (a > b) ? 1.0 : 0.0; break;
            case OP_GE: a = (a >= b) ? 1.0 : 0.0; break;
            case OP_AND: a = (a != 0.0 && b != 0.0) ? 1.0 : 0.0; break;
            case OP_OR: a = (a != 0.0 || b != 0.0) ? 1.0 : 0.0; break;
            case OP_ATAN2: a = PetscAtan2Real(a, b); break;
            case OP_MAXIMUM: a = Extreme(a, b, PETSC_TRUE); break;
            case OP_MINIMUM: a = Extreme(a, b, PETSC_FALSE); break;
            case OP_BESSELJ: a = Bessel(a, b, PETSC_FALSE); break;
            case OP_BESSELY: a = Bessel(a, b, PETSC_TRUE); break;
            default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unknown expression opcode %d.", (int)instruction->op);
            }
            stack[top - 1] = a;
        }
    }
    *result = stack[0];
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvExpressionDestroy()

PetscErrorCode PicurvExpressionDestroy ( PicurvExpression ** expression )

Free a compiled expression.

Parameters

[in,out] expression Expression to free; set to NULL.

Returns: Zero on success.

Free a compiled expression.

See also: PicurvExpressionDestroy()

Definition at line 789 of file expression.c.

{
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    if (!expression || !*expression) { PROFILE_HELPER_END; PetscFunctionReturn(0); }
    PetscCall(PetscFree((*expression)->code));
    PetscCall(PetscFree(*expression));
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

◆ PicurvExpressionUsesVariable()

PetscErrorCode PicurvExpressionUsesVariable	(	const PicurvExpression *	expression,
		PetscInt	index,
		PetscBool *	uses
	)

Whether a compiled expression reads one of its variables.

Parameters

[in]	expression	Compiled expression.
[in]	index	Position of the variable in the compile-time name list.
[out]	uses	True when some instruction loads that variable.

Returns: Zero on success.

Whether a compiled expression reads one of its variables.

See also: PicurvExpressionUsesVariable()

Definition at line 805 of file expression.c.

{
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    *uses = PETSC_FALSE;
    for (PetscInt i = 0; i < expression->length && !*uses; ++i) {
        if (expression->code[i].op == OP_VAR && expression->code[i].arg == index) *uses = PETSC_TRUE;
    }
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

◆ PicurvVelocitySourceCreate()

PetscErrorCode PicurvVelocitySourceCreate	(	const char *	prefix,
		PetscBool	on_face,
		PetscBool	has_data,
		PicurvVelocitySource *	source
	)

Read and compile a velocity source from the options database.

Reads -<prefix>_frame (cartesian, curvilinear, on a face normal, otherwise contravariant), one expression per value, -<prefix>_expr_<c>, written in physical units over the names of PicurvSourceVariable, and the optional -<prefix>_bulk_velocity, an expression in t. Random draws are not available.

Parameters

[in]	prefix	Option prefix without its dash, e.g. `ic_source` or `face_source_0`.
[in]	on_face	Whether the target is a face, where the normal frame means something.
[in]	has_data	Whether the target also has staged data, which makes the expressions optional; without a frame, data on a face is a normal speed.
[out]	source	Compiled source; release it with PicurvVelocitySourceDestroy().

Returns: Zero on success; PETSC_ERR_USER_INPUT for a missing option or unknown frame, PETSC_ERR_ARG_WRONG for a malformed expression.

Read and compile a velocity source from the options database.

See also: PicurvVelocitySourceCreate()

Definition at line 830 of file expression.c.

{
    static const char *const kTime[] = {"t"};
    char      option[PETSC_MAX_PATH_LEN], frame[32], text[PICURV_SOURCE_EXPRESSION_LENGTH];
    PetscBool found = PETSC_FALSE;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCall(PetscMemzero(source, sizeof(*source)));
    PetscCall(PetscSNPrintf(option, sizeof(option), "-%s_frame", prefix));
    PetscCall(PetscOptionsGetString(NULL, NULL, option, frame, sizeof(frame), &found));
    PetscCheck(found || has_data, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Velocity source option %s is required.", option);
    if (!found) PetscCall(PetscStrncpy(frame, on_face ? "normal" : "cartesian", sizeof(frame)));
    if (!strcmp(frame, "cartesian")) source->frame = PICURV_SOURCE_FRAME_CARTESIAN;
    else if (!strcmp(frame, "curvilinear")) source->frame = PICURV_SOURCE_FRAME_CURVILINEAR;
    else if (!strcmp(frame, "normal") && on_face) source->frame = PICURV_SOURCE_FRAME_NORMAL;
    else if (!strcmp(frame, "contravariant") && !on_face) source->frame = PICURV_SOURCE_FRAME_CONTRAVARIANT;
    else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "%s is '%s'; expected cartesian, curvilinear, or %s.",
                 option, frame, on_face ? "normal" : "contravariant");
    source->count = (source->frame == PICURV_SOURCE_FRAME_NORMAL) ? 1 : 3;
    for (PetscInt c = 0; c < source->count; ++c) {
        PetscBool uses_t = PETSC_FALSE;
        PetscCall(PetscSNPrintf(option, sizeof(option), "-%s_expr_%d", prefix, (int)c));
        PetscCall(PetscOptionsGetString(NULL, NULL, option, text, sizeof(text), &found));
        if (!found && has_data && c == 0) { source->count = 0; break; }
        PetscCheck(found, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Velocity source option %s is required.", option);
        PetscCall(PicurvExpressionCompile(text, kPicurvSourceVariableNames,
                                          on_face ? PICURV_SOURCE_VAR_COUNT : PICURV_SOURCE_VAR_CELL_COUNT, PETSC_FALSE,
                                          &source->component[c]));
        PetscCall(PicurvExpressionUsesVariable(source->component[c], PICURV_SOURCE_VAR_T, &uses_t));
        source->time_dependent = (PetscBool)(source->time_dependent || uses_t);
    }
    PetscCall(PetscSNPrintf(option, sizeof(option), "-%s_bulk_velocity", prefix));
    PetscCall(PetscOptionsGetString(NULL, NULL, option, text, sizeof(text), &found));
    if (!found && on_face) {
        PetscCall(PetscSNPrintf(option, sizeof(option), "-%s_flow_rate", prefix));
        PetscCall(PetscOptionsGetString(NULL, NULL, option, text, sizeof(text), &source->bulk_is_flow_rate));
        found = source->bulk_is_flow_rate;
    }
    if (found) {
        PetscCall(PicurvExpressionCompile(text, kTime, 1, PETSC_FALSE, &source->bulk));
        PetscCall(PicurvExpressionUsesVariable(source->bulk, 0, &source->bulk_time_dependent));
    }
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvVelocitySourceBulkFactor()

PetscErrorCode PicurvVelocitySourceBulkFactor	(	const PicurvVelocitySource *	source,
		PetscReal	physical_time,
		const PetscReal	local[3],
		const char *	where,
		PetscReal *	factor
	)

Factor on a source's expressions that makes a target's whole mean equal its bulk velocity.

Solves factor * expressions + data = bulk, over the target's points and weight. Collective.

Parameters

[in]	source	Source with a bulk velocity.
[in]	physical_time	Time the bulk velocity is read at [s].
[in]	local	This rank's weighted sums, physical: the expressions' speed, the staged data's speed along it, and the weight, an area [m^2] for a flow rate.
[in]	where	The target, for the error messages.
[out]	factor	Factor on the expressions.

Returns: Zero on success; PETSC_ERR_ARG_WRONGSTATE for no weight or a zero mean.

Factor on a source's expressions that makes a target's whole mean equal its bulk velocity.

See also: PicurvVelocitySourceBulkFactor()

Definition at line 882 of file expression.c.

{
    PetscReal global[3], target;
 
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    PetscCallMPI(MPI_Allreduce(local, global, 3, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD));
    PetscCheck(global[2] > 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
               "%s has no open area or volume to set a bulk velocity on.", where);
    PetscCheck(global[0] != 0.0, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
               "%s: the source's mean speed is zero, so it cannot be scaled to a bulk velocity.", where);
    PetscCall(PicurvExpressionEvaluate(source->bulk, &physical_time, NULL, &target));
    if (source->bulk_is_flow_rate) target /= global[2];
    *factor = (target * global[2] - global[1]) / global[0];
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "%s: bulk velocity %g, source mean %g, data mean %g, factor %g.\n", where,
              (double)target, (double)(global[0] / global[2]), (double)(global[1] / global[2]), (double)*factor);
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvSourceSetVariables()

static void PicurvSourceSetVariables	(	PetscReal *	values,
		Cmpnts	point,
		Cmpnts	fraction,
		Cmpnts	length,
		PetscReal	L_ref
	)

inlinestatic

Fill a source's position and grid-line variables at one point.

Parameters

[out]	values	Array of PICURV_SOURCE_VAR_COUNT; `t` is left to the caller.
[in]	point	Solver-unit position.
[in]	fraction	Grid-line fractions there, from ComputeGridLineFractions().
[in]	length	Solver-unit grid-line lengths there.
[in]	L_ref	Length reference [m].

Definition at line 206 of file expression.h.

{
    values[PICURV_SOURCE_VAR_X] = point.x * L_ref;
    values[PICURV_SOURCE_VAR_Y] = point.y * L_ref;
    values[PICURV_SOURCE_VAR_Z] = point.z * L_ref;
    values[PICURV_SOURCE_VAR_XI] = fraction.x;
    values[PICURV_SOURCE_VAR_ETA] = fraction.y;
    values[PICURV_SOURCE_VAR_ZETA] = fraction.z;
    values[PICURV_SOURCE_VAR_L_XI] = length.x * L_ref;
    values[PICURV_SOURCE_VAR_L_ETA] = length.y * L_ref;
    values[PICURV_SOURCE_VAR_L_ZETA] = length.z * L_ref;
}

Here is the caller graph for this function:

◆ PicurvSourceFaceSlotsCreate()

PetscErrorCode PicurvSourceFaceSlotsCreate	(	UserCtx *	user,
		BCFace	face_id,
		Vec	fraction,
		Vec	length,
		PicurvSourceFaceSlots *	slots
	)

Gather each owned slot of a boundary face: face-centre position, the grid-line fractions and lengths of the interior cell behind it, and its directions.

A curvilinear component is the speed through the grid face, as on every target, so along the face's own axis the direction follows the face metric rather than the cell's edge.

Parameters

[in]	user	Block with current metric and face-centre ghosts.
[in]	face_id	Face.
[in]	fraction	Grid-line fractions, from ComputeGridLineFractions().
[in]	length	Grid-line lengths, from ComputeGridLineFractions().
[out]	slots	The slots; release with PicurvSourceFaceSlotsDestroy().

Returns: PETSc error code.

Gather each owned slot of a boundary face: face-centre position, the grid-line fractions and lengths of the interior cell behind it, and its directions.

See also: PicurvSourceFaceSlotsCreate()

Definition at line 1124 of file expression.c.

{
    const DMDALocalInfo *info = &user->info;
    const PetscInt start[3] = {info->xs, info->ys, info->zs};
    const PetscInt end[3] = {info->xs + info->xm, info->ys + info->ym, info->zs + info->zm};
    const PetscInt nodes[3] = {info->mx, info->my, info->mz};
    const PetscReal L_ref = user->simCtx->scaling.L_ref;
    const PetscBool negative = (PetscBool)((PetscInt)face_id % 2 == 0);
    PetscBool can_service = PETSC_FALSE;
    PetscInt  n_slots;
    Cmpnts    face_centre;
    PetscReal face_area;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCall(CalculateFaceCenterAndArea(user, face_id, &face_centre, &face_area));
    slots->normal = (PetscInt)face_id / 2;
    slots->a = (slots->normal == 0) ? 1 : 0;
    slots->b = (slots->normal == 2) ? 1 : 2;
    const PetscInt normal = slots->normal, a = slots->a, b = slots->b;
    PetscCall(CanRankServiceFace(info, user->IM, user->JM, user->KM, face_id, &can_service));
    slots->layer = negative ? start[normal] : nodes[normal] - 2;
    slots->lo[0] = (start[a] == 0) ? 1 : start[a];
    slots->lo[1] = (start[b] == 0) ? 1 : start[b];
    slots->n[0] = can_service ? ((end[a] == nodes[a]) ? nodes[a] - 1 : end[a]) - slots->lo[0] : 0;
    slots->n[1] = can_service ? ((end[b] == nodes[b]) ? nodes[b] - 1 : end[b]) - slots->lo[1] : 0;
    n_slots = slots->n[0] * slots->n[1];
    PetscCall(PetscMalloc2(n_slots, &slots->values, n_slots, &slots->edges));
 
    if (n_slots > 0) {
        /* The interior cell behind the slot layer. */
        const PetscInt cell = negative ? slots->layer + 1 : slots->layer;
        const Vec centres = (normal == 0) ? user->lCentx : (normal == 1) ? user->lCenty : user->lCentz;
        const Cmpnts ***cent, ***frac, ***len, ***csi, ***eta, ***zet;
        const PetscReal ***aj;
        PetscCall(DMDAVecGetArrayRead(user->fda, centres, &cent));
        PetscCall(DMDAVecGetArrayRead(user->fda, fraction, &frac));
        PetscCall(DMDAVecGetArrayRead(user->fda, length, &len));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->lCsi, &csi));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->lEta, &eta));
        PetscCall(DMDAVecGetArrayRead(user->fda, user->lZet, &zet));
        PetscCall(DMDAVecGetArrayRead(user->da, user->lAj, &aj));
        for (PetscInt q = 0; q < slots->n[1]; ++q)
        for (PetscInt p = 0; p < slots->n[0]; ++p) {
            const PetscInt s = p + slots->n[0] * q;
            PetscInt f[3], c[3];
            f[normal] = slots->layer; f[a] = slots->lo[0] + p; f[b] = slots->lo[1] + q;
            c[normal] = cell; c[a] = f[a]; c[b] = f[b];
            PicurvSourceSetVariables(slots->values[s], cent[f[2]][f[1]][f[0]], frac[c[2]][c[1]][c[0]], len[c[2]][c[1]][c[0]], L_ref);
            slots->values[s][PICURV_SOURCE_VAR_T] = 0.0;
            /* On the face itself the fraction along the normal line is exactly 0 or 1. */
            slots->values[s][PICURV_SOURCE_VAR_XI + normal] = negative ? 0.0 : 1.0;
            {
                const Cmpnts X = cent[f[2]][f[1]][f[0]];
                const PetscReal dx = X.x - face_centre.x, dy = X.y - face_centre.y, dz = X.z - face_centre.z;
                slots->values[s][PICURV_SOURCE_VAR_R] = PetscSqrtReal(dx * dx + dy * dy + dz * dz) * L_ref;
            }
            PetscCall(ComputeCellEdgeVectors(aj[c[2]][c[1]][c[0]], csi[c[2]][c[1]][c[0]], eta[c[2]][c[1]][c[0]], zet[c[2]][c[1]][c[0]], slots->edges[s]));
            const Cmpnts m = (normal == 0) ? csi[f[2]][f[1]][f[0]] : (normal == 1) ? eta[f[2]][f[1]][f[0]] : zet[f[2]][f[1]][f[0]];
            const Cmpnts through = {user->GridOrientation * m.x, user->GridOrientation * m.y, user->GridOrientation * m.z};
            slots->edges[s][normal] = through;
        }
        PetscCall(DMDAVecRestoreArrayRead(user->da, user->lAj, &aj));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, length, &len));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, fraction, &frac));
        PetscCall(DMDAVecRestoreArrayRead(user->fda, centres, &cent));
    }
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvSourceFaceSlotsDestroy()

PetscErrorCode PicurvSourceFaceSlotsDestroy ( PicurvSourceFaceSlots * slots )

Release what PicurvSourceFaceSlotsCreate() allocated.

Parameters

[in,out] slots The slots.

Returns: PETSc error code.

Release what PicurvSourceFaceSlotsCreate() allocated.

See also: PicurvSourceFaceSlotsDestroy()

Definition at line 1203 of file expression.c.

{
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    PetscCall(PetscFree2(slots->values, slots->edges));
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

◆ PicurvVelocitySourceEvaluate()

PetscErrorCode PicurvVelocitySourceEvaluate	(	const PicurvVelocitySource *	source,
		const PetscReal *	values,
		const Cmpnts	edges[3],
		Cmpnts *	value
	)

Evaluate a source at one point.

Parameters

[in]	source	Compiled source.
[in]	values	One value per PicurvSourceVariable.
[in]	edges	Cell edge vectors along xi, eta, zeta (any length); read only in the curvilinear frame.
[out]	value	Physical Cartesian velocity, or for the normal frame the speed in `x`.

Returns: Zero on success; PETSC_ERR_FP for a non-finite value.

Evaluate a source at one point.

See also: PicurvVelocitySourceEvaluate()

Definition at line 907 of file expression.c.

{
    PetscReal v[3] = {0.0, 0.0, 0.0};
 
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    for (PetscInt c = 0; c < source->count; ++c) {
        PetscCall(PicurvExpressionEvaluate(source->component[c], values, NULL, &v[c]));
        PetscCheck(!PetscIsInfOrNanReal(v[c]), PETSC_COMM_SELF, PETSC_ERR_FP,
                   "Velocity source value %" PetscInt_FMT " is %g at (x, y, z) = (%g, %g, %g), t = %g.",
                   c, (double)v[c], (double)values[PICURV_SOURCE_VAR_X], (double)values[PICURV_SOURCE_VAR_Y],
                   (double)values[PICURV_SOURCE_VAR_Z], (double)values[PICURV_SOURCE_VAR_T]);
    }
    if (source->frame == PICURV_SOURCE_FRAME_CURVILINEAR) {
        /* Values are components along the unit grid-line directions. */
        value->x = value->y = value->z = 0.0;
        for (PetscInt a = 0; a < 3; ++a) {
            const PetscReal norm = PetscSqrtReal(edges[a].x * edges[a].x + edges[a].y * edges[a].y + edges[a].z * edges[a].z);
            PetscCheck(norm > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "A cell edge has zero length.");
            value->x += v[a] * edges[a].x / norm;
            value->y += v[a] * edges[a].y / norm;
            value->z += v[a] * edges[a].z / norm;
        }
    } else {
        value->x = v[0];
        value->y = v[1];
        value->z = v[2];
    }
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvSourceSetCellPressure()

PetscErrorCode PicurvSourceSetCellPressure	(	UserCtx *	user,
		const PicurvExpression *	pressure,
		Vec	fraction,
		Vec	length,
		PetscReal	physical_time
	)

Set a block's pressure at every interior cell centre from an expression in Pa.

Dummy cells are left to the caller's boundary pass.

Parameters

[in]	user	Block with current cell centres.
[in]	pressure	Expression over the cell variables, physical.
[in]	fraction	Grid-line fractions, from ComputeGridLineFractions().
[in]	length	Grid-line lengths, from ComputeGridLineFractions().
[in]	physical_time	Time [s].

Returns: PETSc error code.

Set a block's pressure at every interior cell centre from an expression in Pa.

See also: PicurvSourceSetCellPressure()

Definition at line 1086 of file expression.c.

{
    const DMDALocalInfo *info = &user->info;
    const ScalingCtx *scaling = &user->simCtx->scaling;
    const PetscInt xs = info->xs, xe = info->xs + info->xm, mx = info->mx;
    const PetscInt ys = info->ys, ye = info->ys + info->ym, my = info->my;
    const PetscInt zs = info->zs, ze = info->zs + info->zm, mz = info->mz;
    const Cmpnts ***cent, ***frac, ***len;
    PetscReal ***p, values[PICURV_SOURCE_VAR_COUNT];
 
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Cent, &cent));
    PetscCall(DMDAVecGetArrayRead(user->fda, fraction, &frac));
    PetscCall(DMDAVecGetArrayRead(user->fda, length, &len));
    PetscCall(DMDAVecGetArray(user->da, user->P, &p));
    values[PICURV_SOURCE_VAR_T] = physical_time;
    for (PetscInt k = (zs == 0) ? 1 : zs; k < ((ze == mz) ? mz - 1 : ze); ++k)
    for (PetscInt j = (ys == 0) ? 1 : ys; j < ((ye == my) ? my - 1 : ye); ++j)
    for (PetscInt i = (xs == 0) ? 1 : xs; i < ((xe == mx) ? mx - 1 : xe); ++i) {
        PetscReal value;
        PicurvSourceSetVariables(values, cent[k][j][i], frac[k][j][i], len[k][j][i], scaling->L_ref);
        PetscCall(PicurvExpressionEvaluate(pressure, values, NULL, &value));
        p[k][j][i] = value / scaling->P_ref;
    }
    PetscCall(DMDAVecRestoreArray(user->da, user->P, &p));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, length, &len));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, fraction, &frac));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Cent, &cent));
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvSourceCellBulkFactor()

PetscErrorCode PicurvSourceCellBulkFactor	(	UserCtx *	user,
		const PicurvVelocitySource *	source,
		Vec	fraction,
		Vec	length,
		PetscReal	physical_time,
		PetscReal *	factor
	)

Factor on a source's expressions that makes the volume mean of the whole field, along the source's own direction at each cell, equal its bulk velocity.

Collective.

Parameters

[in]	user	Block; its Ucat holds whatever is staged.
[in]	source	Source with a bulk velocity.
[in]	fraction	Grid-line fractions, from ComputeGridLineFractions().
[in]	length	Grid-line lengths, from ComputeGridLineFractions().
[in]	physical_time	Time [s].
[out]	factor	Factor on the expressions.

Returns: PETSc error code.

Factor on a source's expressions that makes the volume mean of the whole field, along the source's own direction at each cell, equal its bulk velocity.

See also: PicurvSourceCellBulkFactor()

Definition at line 944 of file expression.c.

{
    const DMDALocalInfo *info = &user->info;
    const PetscReal L_ref = user->simCtx->scaling.L_ref, U_ref = user->simCtx->scaling.U_ref;
    const PetscInt xs = info->xs, xe = info->xs + info->xm, mx = info->mx;
    const PetscInt ys = info->ys, ye = info->ys + info->ym, my = info->my;
    const PetscInt zs = info->zs, ze = info->zs + info->zm, mz = info->mz;
    PetscReal values[PICURV_SOURCE_VAR_COUNT], local[3] = {0.0, 0.0, 0.0};
    const Cmpnts ***cent, ***frac, ***len, ***csi, ***eta, ***zet, ***ucat;
    const PetscReal ***aj;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Cent, &cent));
    PetscCall(DMDAVecGetArrayRead(user->fda, fraction, &frac));
    PetscCall(DMDAVecGetArrayRead(user->fda, length, &len));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lCsi, &csi));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lEta, &eta));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lZet, &zet));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lAj, &aj));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->Ucat, &ucat));
    values[PICURV_SOURCE_VAR_T] = physical_time;
    for (PetscInt k = (zs == 0) ? 1 : zs; k < ((ze == mz) ? mz - 1 : ze); ++k)
    for (PetscInt j = (ys == 0) ? 1 : ys; j < ((ye == my) ? my - 1 : ye); ++j)
    for (PetscInt i = (xs == 0) ? 1 : xs; i < ((xe == mx) ? mx - 1 : xe); ++i) {
        const PetscReal volume = 1.0 / PetscAbsReal(aj[k][j][i]);
        Cmpnts edges[3], v;
        PicurvSourceSetVariables(values, cent[k][j][i], frac[k][j][i], len[k][j][i], L_ref);
        PetscCall(ComputeCellEdgeVectors(aj[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i], edges));
        PetscCall(PicurvVelocitySourceEvaluate(source, values, edges, &v));
        const PetscReal speed = PetscSqrtReal(v.x * v.x + v.y * v.y + v.z * v.z);
        local[0] += speed * volume;
        if (speed > 0.0) local[1] += U_ref * (ucat[k][j][i].x * v.x + ucat[k][j][i].y * v.y +
                                              ucat[k][j][i].z * v.z) / speed * volume;
        local[2] += volume;
    }
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Ucat, &ucat));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lAj, &aj));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, length, &len));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, fraction, &frac));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->Cent, &cent));
    PetscCall(PicurvVelocitySourceBulkFactor(source, physical_time, local, "The field", factor));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvSourceAddFluxes()

PetscErrorCode PicurvSourceAddFluxes	(	UserCtx *	user,
		const PicurvVelocitySource *	source,
		Vec	fraction,
		Vec	length,
		PetscReal	physical_time,
		PetscReal	factor,
		PetscBool	boundary_faces
	)

Add a source's flux to every owned face of Ucont, evaluating it at each face centre.

A Cartesian velocity is dotted with the face metric; a curvilinear component is the speed through the face, times its area; a contravariant component is the flux itself.

Parameters

[in,out]	user	Block whose Ucont receives the fluxes.
[in]	source	Source with expressions.
[in]	fraction	Grid-line fractions, from ComputeGridLineFractions().
[in]	length	Grid-line lengths, from ComputeGridLineFractions().
[in]	physical_time	Time [s].
[in]	factor	Factor on the expressions.
[in]	boundary_faces	Whether the boundary - faces are set too; each cell's + face always is.

Returns: PETSc error code.

Add a source's flux to every owned face of Ucont, evaluating it at each face centre.

See also: PicurvSourceAddFluxes()

Definition at line 998 of file expression.c.

{
    const DMDALocalInfo *info = &user->info;
    const PetscReal L_ref = user->simCtx->scaling.L_ref, U_ref = user->simCtx->scaling.U_ref;
    const PetscInt start[3] = {info->xs, info->ys, info->zs};
    const PetscInt end[3] = {info->xs + info->xm, info->ys + info->ym, info->zs + info->zm};
    const PetscInt nodes[3] = {info->mx, info->my, info->mz};
    const PetscInt lo[3] = {(start[0] == 0) ? 1 : start[0], (start[1] == 0) ? 1 : start[1], (start[2] == 0) ? 1 : start[2]};
    const PetscInt hi[3] = {(end[0] == nodes[0]) ? nodes[0] - 1 : end[0], (end[1] == nodes[1]) ? nodes[1] - 1 : end[1],
                            (end[2] == nodes[2]) ? nodes[2] - 1 : end[2]};
    const Cmpnts ***frac, ***len, ***csi, ***eta, ***zet, ***centx, ***centy, ***centz;
    const PetscReal ***aj;
    Cmpnts ***ucont;
    PetscReal values[PICURV_SOURCE_VAR_COUNT];
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    PetscCall(DMDAVecGetArrayRead(user->fda, fraction, &frac));
    PetscCall(DMDAVecGetArrayRead(user->fda, length, &len));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lCsi, &csi));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lEta, &eta));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lZet, &zet));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lCentx, &centx));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lCenty, &centy));
    PetscCall(DMDAVecGetArrayRead(user->fda, user->lCentz, &centz));
    PetscCall(DMDAVecGetArrayRead(user->da, user->lAj, &aj));
    PetscCall(DMDAVecGetArray(user->fda, user->Ucont, &ucont));
    values[PICURV_SOURCE_VAR_T] = physical_time;
    for (PetscInt k = lo[2]; k < hi[2]; ++k)
    for (PetscInt j = lo[1]; j < hi[1]; ++j)
    for (PetscInt i = lo[0]; i < hi[0]; ++i) {
        const PetscInt cell[3] = {i, j, k};
        const PetscReal lengths[3] = {len[k][j][i].x, len[k][j][i].y, len[k][j][i].z};
        Cmpnts edges[3];
        PetscCall(ComputeCellEdgeVectors(aj[k][j][i], csi[k][j][i], eta[k][j][i], zet[k][j][i], edges));
        for (PetscInt a = 0; a < 3; ++a) {
            const PetscReal half = 0.5 * PetscSqrtReal(edges[a].x * edges[a].x + edges[a].y * edges[a].y +
                                                       edges[a].z * edges[a].z) / lengths[a];
            /* Each cell owns its + face; the first cell of a line also owns the boundary - face. */
            const PetscInt sides = (boundary_faces && start[a] == 0 && cell[a] == 1) ? 2 : 1;
            for (PetscInt side = 0; side < sides; ++side) {
                const PetscInt offset = side ? -1 : 0;
                const PetscInt fi = i + (a == 0 ? offset : 0), fj = j + (a == 1 ? offset : 0), fk = k + (a == 2 ? offset : 0);
                const Cmpnts metric = (a == 0) ? csi[fk][fj][fi] : (a == 1) ? eta[fk][fj][fi] : zet[fk][fj][fi];
                const Cmpnts centre = (a == 0) ? centx[fk][fj][fi] : (a == 1) ? centy[fk][fj][fi] : centz[fk][fj][fi];
                PetscReal flux;
                PicurvSourceSetVariables(values, centre, frac[k][j][i], len[k][j][i], L_ref);
                values[PICURV_SOURCE_VAR_XI + a] += side ? -half : half;
                if (source->frame == PICURV_SOURCE_FRAME_CONTRAVARIANT) {
                    PetscReal value;
                    PetscCall(PicurvExpressionEvaluate(source->component[a], values, NULL, &value));
                    flux = factor * value / (U_ref * L_ref * L_ref);
                } else if (source->frame == PICURV_SOURCE_FRAME_CURVILINEAR) {
                    const PetscReal area = PetscSqrtReal(metric.x * metric.x + metric.y * metric.y + metric.z * metric.z);
                    PetscReal value;
                    PetscCall(PicurvExpressionEvaluate(source->component[a], values, NULL, &value));
                    flux = factor * (value / U_ref) * user->GridOrientation * area;
                } else {
                    Cmpnts v;
                    PetscCall(PicurvVelocitySourceEvaluate(source, values, edges, &v));
                    flux = (factor * v.x / U_ref) * metric.x + (factor * v.y / U_ref) * metric.y +
                           (factor * v.z / U_ref) * metric.z;
                }
                if (a == 0) ucont[fk][fj][fi].x += flux;
                if (a == 1) ucont[fk][fj][fi].y += flux;
                if (a == 2) ucont[fk][fj][fi].z += flux;
            }
        }
    }
    PetscCall(DMDAVecRestoreArray(user->fda, user->Ucont, &ucont));
    PetscCall(DMDAVecRestoreArrayRead(user->da, user->lAj, &aj));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCentz, &centz));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCenty, &centy));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCentx, &centx));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, length, &len));
    PetscCall(DMDAVecRestoreArrayRead(user->fda, fraction, &frac));
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PicurvVelocitySourceDestroy()

PetscErrorCode PicurvVelocitySourceDestroy ( PicurvVelocitySource * source )

Free a source's compiled expressions and bulk velocity.

Parameters

[in,out] source Source to release; left empty.

Returns: Zero on success.

Free a source's compiled expressions and bulk velocity.

See also: PicurvVelocitySourceDestroy()

Definition at line 1216 of file expression.c.

{
    PetscFunctionBeginUser;
    PROFILE_HELPER_BEGIN;
    for (PetscInt c = 0; c < 3; ++c) PetscCall(PicurvExpressionDestroy(&source->component[c]));
    PetscCall(PicurvExpressionDestroy(&source->bulk));
    source->count = 0;
    PROFILE_HELPER_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

◆ kPicurvSourceVariableNames

const char* const kPicurvSourceVariableNames[PICURV_SOURCE_VAR_COUNT]

extern

Names of PicurvSourceVariable, in order.

Definition at line 821 of file expression.c.

821 {

822 "x", "y", "z", "t", "xi", "eta", "zeta", "L_xi", "L_eta", "L_zeta", "r"};

Data Structures

Typedefs

Enumerations

Functions

Variables

Detailed Description

Data Structure Documentation

◆ PicurvExpressionDrawKey

◆ PicurvVelocitySource

◆ PicurvSourceFaceSlots

Typedef Documentation

◆ PicurvExpression

Enumeration Type Documentation

◆ PicurvDrawSpace

◆ PicurvSourceVariable

◆ PicurvSourceFrame

Function Documentation

◆ PicurvMix64()

◆ PicurvUnitFromHash()

◆ PicurvParticleDrawHash()

◆ PicurvExpressionCompile()

◆ PicurvExpressionEvaluate()

◆ PicurvExpressionDestroy()

◆ PicurvExpressionUsesVariable()

◆ PicurvVelocitySourceCreate()

◆ PicurvVelocitySourceBulkFactor()

◆ PicurvSourceSetVariables()

◆ PicurvSourceFaceSlotsCreate()

◆ PicurvSourceFaceSlotsDestroy()

◆ PicurvVelocitySourceEvaluate()

◆ PicurvSourceSetCellPressure()

◆ PicurvSourceCellBulkFactor()

◆ PicurvSourceAddFluxes()

◆ PicurvVelocitySourceDestroy()

Variable Documentation

◆ kPicurvSourceVariableNames