AnalyticalSolutions.h File Reference

#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscdmcomposite.h>
#include "variables.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "setup.h"
#include "interpolation.h"

Include dependency graph for AnalyticalSolutions.h:

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

Go to the source code of this file.

Macros
#define	__FUNCT__   "AnalyticalSolutionEngine"

Functions
PetscErrorCode	SetAnalyticalGridInfo (UserCtx *user)
	Sets the grid domain and resolution for analytical solution cases.
PetscBool	AnalyticalTypeRequiresCustomGeometry (const char *analytical_type)
	Reports whether an analytical type requires custom geometry/decomposition logic.
PetscBool	AnalyticalTypeSupportsInterpolationError (const char *analytical_type)
	Reports whether an analytical type has a non-trivial velocity field for which interpolation error measurement is meaningful.
PetscErrorCode	AnalyticalSolutionEngine (SimCtx *simCtx)
	Dispatches to the appropriate analytical solution function based on simulation settings.
PetscErrorCode	SetAnalyticalSolutionForParticles (Vec tempVec, SimCtx *simCtx)
	Applies the analytical solution to particle velocity vector.
PetscErrorCode	EvaluateAnalyticalScalarProfile (const SimCtx *simCtx, PetscReal x, PetscReal y, PetscReal z, PetscReal t, PetscReal *value)
	Evaluates the configured verification scalar profile at one physical point.
PetscErrorCode	SetAnalyticalScalarFieldOnParticles (UserCtx *user, ParticleFieldId particle_field_id)
	Writes the configured verification scalar profile onto a particle swarm scalar field.
PetscErrorCode	SetAnalyticalScalarFieldAtCellCenters (UserCtx *user, Vec targetVec)
	Writes the configured verification scalar profile at physical cell centers into a scalar Vec.

Macro Definition Documentation

__FUNCT__

#define __FUNCT__   "AnalyticalSolutionEngine"

Definition at line 73 of file AnalyticalSolutions.h.

Function Documentation

SetAnalyticalGridInfo()

PetscErrorCode SetAnalyticalGridInfo

(

UserCtx *

user

)

Sets the grid domain and resolution for analytical solution cases.

This function is called when eulerianSource is "analytical". It is responsible for automatically configuring the grid based on the chosen AnalyticalSolutionType.

TGV3D Multi-Block Decomposition

If the analytical solution is "TGV3D", this function automatically decomposes the required [0, 2*PI] physical domain among the available blocks.

• **Single Block (nblk=1):** The single block is assigned the full [0, 2*PI] domain.
• **Multiple Blocks (nblk>1):** It requires that the number of blocks be a perfect square (e.g., 4, 9, 16). It then arranges the blocks in a sqrt(nblk) by sqrt(nblk) grid in the X-Y plane, partitioning the [0, 2*PI] domain in X and Y accordingly. The Z domain for all blocks remains [0, 2*PI]. If nblk is not a perfect square, the simulation is aborted with an error.

Grid resolution (IM/JM/KM) is expected to be pre-populated in user before this function is called.

Parameters

user	Pointer to the UserCtx for a specific block. The function will populate the geometric fields (IM, JM, KM, Min_X, Max_X, etc.) within this struct.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Sets the grid domain and resolution for analytical solution cases.

Local to this translation unit.

Definition at line 88 of file AnalyticalSolutions.c.

{
    SimCtx         *simCtx = user->simCtx;
    PetscInt       nblk = simCtx->block_number;
    PetscInt       block_index = user->_this;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if (!AnalyticalTypeRequiresCustomGeometry(simCtx->AnalyticalSolutionType)) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                 "SetAnalyticalGridInfo called for analytical type '%s' that does not require custom geometry.",
                 simCtx->AnalyticalSolutionType);
    }
    if (user->IM <= 0 || user->JM <= 0 || user->KM <= 0) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE,
                "Analytical grid resolution is not initialized. Ensure IM/JM/KM are preloaded before SetAnalyticalGridInfo.");
    }
 
    if (strcmp(simCtx->AnalyticalSolutionType, "TGV3D") == 0) {
        LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "Rank %d: Configuring grid for TGV3D analytical solution, block %d.\n", simCtx->rank, block_index);
 
        if (nblk == 1) {
            // --- Single Block Case ---
            if (block_index == 0) {
                LOG_ALLOW(GLOBAL, LOG_INFO, "Single block detected. Setting domain to [0, 2*PI].\n");
            }
            user->Min_X = 0.0; user->Max_X = 2.0 * PETSC_PI;
            user->Min_Y = 0.0; user->Max_Y = 2.0 * PETSC_PI;
            user->Min_Z = 0.0; user->Max_Z = 0.2 * PETSC_PI; //2.0 * PETSC_PI;
 
        } else { // --- Multi-Block Case ---
            PetscReal s = sqrt((PetscReal)nblk);
            
            // Validate that nblk is a perfect square.
            if (fabs(s - floor(s)) > 1e-9) {
                SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP,
                         "\n\n*** CONFIGURATION ERROR FOR TGV3D ***\n"
                         "For multi-block TGV3D cases, the number of blocks must be a perfect square (e.g., 4, 9, 16).\n"
                         "You have specified %d blocks. Please adjust `-block_number`.\n", nblk);
            }
            PetscInt blocks_per_dim = (PetscInt)s;
 
            if (block_index == 0) {
                 LOG_ALLOW(GLOBAL, LOG_INFO, "%d blocks detected. Decomposing domain into a %d x %d grid in the X-Y plane.\n", nblk, blocks_per_dim, blocks_per_dim);
            }
 
            // Determine the (row, col) position of this block in the 2D decomposition
            PetscInt row = block_index / blocks_per_dim;
            PetscInt col = block_index % blocks_per_dim;
 
            // Calculate the width/height of each sub-domain
            PetscReal block_width  = (2.0 * PETSC_PI) / (PetscReal)blocks_per_dim;
            PetscReal block_height = (2.0 * PETSC_PI) / (PetscReal)blocks_per_dim;
            
            // Assign this block its specific sub-domain
            user->Min_X = col * block_width;
            user->Max_X = (col + 1) * block_width;
            user->Min_Y = row * block_height;
            user->Max_Y = (row + 1) * block_height;
            user->Min_Z = 0.0;
            user->Max_Z = 2.0 * PETSC_PI; // Z-domain is not decomposed
        }
    }
    /*
     * --- EXTENSIBILITY HOOK ---
     * To add another analytical case with special grid requirements:
     *
     * else if (strcmp(simCtx->AnalyticalSolutionType, "ChannelFlow") == 0) {
     *     // ... implement logic to set domain for ChannelFlow case ...
     * }
     */
    else {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_UNKNOWN_TYPE,
                 "Analytical type '%s' has no custom geometry implementation.",
                 simCtx->AnalyticalSolutionType);
    }
 
    // We can also read stretching ratios, as they are independent of the domain size
    // For simplicity, we assume uniform grid unless specified.
    user->rx = 1.0; user->ry = 1.0; user->rz = 1.0;
    
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Block %d grid resolution set: IM=%d, JM=%d, KM=%d\n",
              simCtx->rank, block_index, user->IM, user->JM, user->KM);
    LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d: Block %d final bounds: X=[%.4f, %.4f], Y=[%.4f, %.4f], Z=[%.4f, %.4f]\n",
              simCtx->rank, block_index, user->Min_X, user->Max_X, user->Min_Y, user->Max_Y, user->Min_Z, user->Max_Z);
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

AnalyticalTypeRequiresCustomGeometry()

PetscBool AnalyticalTypeRequiresCustomGeometry

(

const char *

analytical_type

)

Reports whether an analytical type requires custom geometry/decomposition logic.

Analytical types returning PETSC_TRUE are expected to route through SetAnalyticalGridInfo. Types returning PETSC_FALSE should use the standard programmatic grid parser fallback.

Parameters

analytical_type

Analytical solution type string.

Returns

PETSC_TRUE if custom geometry is required, PETSC_FALSE otherwise.

Reports whether an analytical type requires custom geometry/decomposition logic.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/AnalyticalSolutions.h.

See also

AnalyticalTypeRequiresCustomGeometry()

Definition at line 63 of file AnalyticalSolutions.c.

{
    if (!analytical_type) return PETSC_FALSE;
    return (strcmp(analytical_type, "TGV3D") == 0) ? PETSC_TRUE : PETSC_FALSE;
}

Here is the caller graph for this function:

AnalyticalTypeSupportsInterpolationError()

PetscBool AnalyticalTypeSupportsInterpolationError

(

const char *

analytical_type

)

Reports whether an analytical type has a non-trivial velocity field for which interpolation error measurement is meaningful.

Types with identically zero velocity (e.g. ZERO_FLOW) return PETSC_FALSE because the interpolation error is trivially zero and uninformative.

Parameters

analytical_type

Analytical solution type string.

Returns

PETSC_TRUE if interpolation error is meaningful, PETSC_FALSE otherwise.

Reports whether an analytical type has a non-trivial velocity field for which interpolation error measurement is meaningful.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/AnalyticalSolutions.h.

See also

AnalyticalTypeSupportsInterpolationError()

Definition at line 75 of file AnalyticalSolutions.c.

{
    if (!analytical_type) return PETSC_FALSE;
    if (strcmp(analytical_type, "ZERO_FLOW") == 0 || strcmp(analytical_type, "UNIFORM_FLOW") == 0) return PETSC_FALSE;
    return PETSC_TRUE;
}

Here is the caller graph for this function:

AnalyticalSolutionEngine()

PetscErrorCode AnalyticalSolutionEngine

(

SimCtx *

simCtx

)

Dispatches to the appropriate analytical solution function based on simulation settings.

This function acts as a router. It reads the AnalyticalSolutionType from the simulation context and calls the corresponding private implementation function (e.g., for Taylor-Green Vortex, lid-driven cavity, etc.). This design keeps the main simulation code clean and makes it easy to add new analytical test cases.

Parameters

simCtx

The main simulation context, containing configuration and state.

Returns

PetscErrorCode 0 on success.

Dispatches to the appropriate analytical solution function based on simulation settings.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/AnalyticalSolutions.h.

See also

AnalyticalSolutionEngine()

Definition at line 192 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    // -- Before any operation, here is defensive test to ensure that the Corner->Center Interpolation method works
    //ierr = TestCornerToCenterInterpolation(&(simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1]->user[0])); 
 
    // --- Dispatch based on the string provided by the user ---
    if (strcmp(simCtx->AnalyticalSolutionType, "TGV3D") == 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Applying Analytical Solution: 3D Taylor-Green Vortex (TGV3D).\n");
        ierr = SetAnalyticalSolution_TGV3D(simCtx); CHKERRQ(ierr);
    }
    else if (strcmp(simCtx->AnalyticalSolutionType, "ZERO_FLOW") == 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Applying Analytical Solution: Zero Background Flow (ZERO_FLOW).\n");
        ierr = SetAnalyticalSolution_ZeroFlow(simCtx); CHKERRQ(ierr);
    }
    else if (strcmp(simCtx->AnalyticalSolutionType, "UNIFORM_FLOW") == 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Applying Analytical Solution: Uniform Background Flow (UNIFORM_FLOW).\n");
        ierr = SetAnalyticalSolution_UniformFlow(simCtx); CHKERRQ(ierr);
    }
    /*
     * --- EXTENSIBILITY HOOK ---
     * To add a new analytical solution (e.g., "ChannelFlow"):
     * 1. Add an `else if` block here:
     *
     *    else if (strcmp(simCtx->AnalyticalSolutionType, "ChannelFlow") == 0) {
     *        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Applying Analytical Solution: Channel Flow.\n");
     *        ierr = SetAnalyticalSolution_ChannelFlow(simCtx); CHKERRQ(ierr);
     *    }
     *
     * 2. Implement the static function `SetAnalyticalSolution_ChannelFlow(SimCtx *simCtx)`
     *    below, following the TGV3D pattern.
     */
    else {
        // If the type is unknown, raise a fatal error to prevent silent failures.
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown AnalyticalSolutionType specified: '%s'", simCtx->AnalyticalSolutionType);
    }
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetAnalyticalSolutionForParticles()

PetscErrorCode SetAnalyticalSolutionForParticles	(	Vec	tempVec,
		SimCtx *	simCtx
	)

Applies the analytical solution to particle velocity vector.

Dispatcher function that calls the appropriate analytical solution based on simCtx->AnalyticalSolutionType. Supports multiple solution types.

Parameters

tempVec	The PETSc Vec containing particle positions which will be used to store velocities.
simCtx	The simulation context.

Returns

PetscErrorCode Returns 0 on success.

Applies the analytical solution to particle velocity vector.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/AnalyticalSolutions.h.

See also

SetAnalyticalSolutionForParticles()

Definition at line 546 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    const char *analytical_type = simCtx->AnalyticalSolutionType[0] ? simCtx->AnalyticalSolutionType : "default";
 
    PetscFunctionBeginUser;
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Type: %s\n", analytical_type);
    
    // Check for specific analytical solution types
    if (strcmp(simCtx->AnalyticalSolutionType, "TGV3D") == 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Using TGV3D solution.\n");
        ierr = SetAnalyticalSolutionForParticles_TGV3D(tempVec, simCtx); CHKERRQ(ierr);
        return 0;
    }
    if (strcmp(simCtx->AnalyticalSolutionType, "UNIFORM_FLOW") == 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Using UNIFORM_FLOW solution.\n");
        ierr = SetAnalyticalSolutionForParticles_UniformFlow(tempVec, simCtx); CHKERRQ(ierr);
        return 0;
    }
    
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

EvaluateAnalyticalScalarProfile()

PetscErrorCode EvaluateAnalyticalScalarProfile	(	const SimCtx *	simCtx,
		PetscReal	x,
		PetscReal	y,
		PetscReal	z,
		PetscReal	t,
		PetscReal *	value
	)

Evaluates the configured verification scalar profile at one physical point.

Parameters

[in]	simCtx	Simulation context providing the scalar verification profile.
[in]	x	Physical x coordinate.
[in]	y	Physical y coordinate.
[in]	z	Physical z coordinate.
[in]	t	Simulation time.
[out]	value	Evaluated scalar value.

Returns

PetscErrorCode 0 on success.

Evaluates the configured verification scalar profile at one physical point.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/AnalyticalSolutions.h.

See also

EvaluateAnalyticalScalarProfile()

Definition at line 609 of file AnalyticalSolutions.c.

{
    PetscFunctionBeginUser;
    (void)t;
    if (!simCtx) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "SimCtx cannot be NULL.");
    if (!simCtx->verificationScalar.enabled) {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE,
                "EvaluateAnalyticalScalarProfile requires verification scalar mode to be enabled.");
    }
    PetscCall(EvaluateConfiguredScalarProfile(&simCtx->verificationScalar, x, y, z, value));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetAnalyticalScalarFieldOnParticles()

PetscErrorCode SetAnalyticalScalarFieldOnParticles	(	UserCtx *	user,
		ParticleFieldId	particle_field_id
	)

Writes the configured verification scalar profile onto a particle swarm scalar field.

Parameters

[in,out]	user	Block-local context providing particle positions.
[in]	particle_field_id	Typed identity of the scalar swarm field to overwrite.

Returns

PetscErrorCode 0 on success.

Writes the configured verification scalar profile onto a particle swarm scalar field.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/AnalyticalSolutions.h.

See also

SetAnalyticalScalarFieldOnParticles()

Definition at line 635 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    PetscInt       nlocal = 0;
    PetscReal     *positions = NULL;
    PetscReal     *scalar_values = NULL;
    const ParticleFieldDescriptor *descriptor = NULL;
    const char     *swarm_field_name = NULL;
 
    PetscFunctionBeginUser;
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx cannot be NULL.");
    if (!user->swarm) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "UserCtx->swarm is NULL.");
 
    ierr = ParticleFieldGetDescriptor(particle_field_id, &descriptor); CHKERRQ(ierr);
    PetscCheck(descriptor->components == 1 && descriptor->data_type == PETSC_REAL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
               "Analytical scalar assignment requires a one-component PETSC_REAL particle field; '%s' has %d components and type %s.",
               descriptor->canonical_name, descriptor->components, PetscDataTypes[descriptor->data_type]);
    swarm_field_name = descriptor->canonical_name;
 
    ierr = DMSwarmGetLocalSize(user->swarm, &nlocal); CHKERRQ(ierr);
    if (nlocal == 0) PetscFunctionReturn(0);
 
    ierr = DMSwarmGetField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void **)&positions); CHKERRQ(ierr);
    ierr = DMSwarmGetField(user->swarm, swarm_field_name, NULL, NULL, (void **)&scalar_values); CHKERRQ(ierr);
 
    for (PetscInt p = 0; p < nlocal; ++p) {
        PetscReal value = 0.0;
        ierr = EvaluateAnalyticalScalarProfile(user->simCtx,
                                               positions[3 * p + 0],
                                               positions[3 * p + 1],
                                               positions[3 * p + 2],
                                               user->simCtx->ti,
                                               &value); CHKERRQ(ierr);
        scalar_values[p] = value;
    }
 
    ierr = DMSwarmRestoreField(user->swarm, swarm_field_name, NULL, NULL, (void **)&scalar_values); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION), NULL, NULL, (void **)&positions); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetAnalyticalScalarFieldAtCellCenters()

PetscErrorCode SetAnalyticalScalarFieldAtCellCenters	(	UserCtx *	user,
		Vec	targetVec
	)

Writes the configured verification scalar profile at physical cell centers into a scalar Vec.

Parameters

[in]	user	Block-local context providing the cell-center coordinates.
[in,out]	targetVec	Scalar Vec on user->da storage to fill with analytical reference values.

Returns

PetscErrorCode 0 on success.

Writes the configured verification scalar profile at physical cell centers into a scalar Vec.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/AnalyticalSolutions.h.

See also

SetAnalyticalScalarFieldAtCellCenters()

Definition at line 685 of file AnalyticalSolutions.c.

{
    PetscErrorCode  ierr;
    PetscReal     ***target = NULL;
    const Cmpnts  ***cent = NULL;
    DMDALocalInfo   info;
    PetscInt        xs, xe, ys, ye, zs, ze, mx, my, mz;
    PetscInt        lxs, lxe, lys, lye, lzs, lze;
 
    PetscFunctionBeginUser;
    if (!user) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "UserCtx cannot be NULL.");
    if (!targetVec) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "targetVec cannot be NULL.");
 
    info = user->info;
    xs = info.xs; xe = info.xs + info.xm;
    ys = info.ys; ye = info.ys + info.ym;
    zs = info.zs; ze = info.zs + info.zm;
    mx = info.mx; my = info.my; mz = info.mz;
    lxs = (xs == 0) ? xs + 1 : xs; lxe = (xe == mx) ? xe - 1 : xe;
    lys = (ys == 0) ? ys + 1 : ys; lye = (ye == my) ? ye - 1 : ye;
    lzs = (zs == 0) ? zs + 1 : zs; lze = (ze == mz) ? ze - 1 : ze;
 
    ierr = VecSet(targetVec, 0.0); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->da, targetVec, &target); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->Cent, &cent); CHKERRQ(ierr);
 
    for (PetscInt k = lzs; k < lze; ++k) {
        for (PetscInt j = lys; j < lye; ++j) {
            for (PetscInt i = lxs; i < lxe; ++i) {
                PetscReal value = 0.0;
                ierr = EvaluateAnalyticalScalarProfile(user->simCtx,
                                                       cent[k][j][i].x,
                                                       cent[k][j][i].y,
                                                       cent[k][j][i].z,
                                                       user->simCtx->ti,
                                                       &value); CHKERRQ(ierr);
                target[k][j][i] = value;
            }
        }
    }
 
    ierr = DMDAVecRestoreArrayRead(user->fda, user->Cent, &cent); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->da, targetVec, &target); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: