AnalyticalSolutions.c File Reference

Implements the analytical solution engine for initializing or driving the simulation. More...

#include "AnalyticalSolutions.h"
#include <petscmath.h>

Include dependency graph for AnalyticalSolutions.c:

Go to the source code of this file.

Macros
#define	__FUNCT__   "SetAnalyticalGridInfo"
#define	__FUNCT__   "AnalyticalSolutionEngine"
#define	__FUNCT__   "SetAnalyticalSolution_TGV3D"
#define	__FUNCT__   "SetAnalyticalSolution_ZeroFlow"
#define	__FUNCT__   "SetAnalyticalSolution_UniformFlow"
#define	__FUNCT__   "SetAnalyticalSolutionForParticles_TGV3D"
#define	__FUNCT__   "SetAnalyticalSolutionForParticles_UniformFlow"
#define	__FUNCT__   "SetAnalyticalSolutionForParticles"
#define	__FUNCT__   "EvaluateAnalyticalScalarProfile"
#define	__FUNCT__   "SetAnalyticalScalarFieldOnParticles"
#define	__FUNCT__   "SetAnalyticalScalarFieldAtCellCenters"

Functions
static PetscErrorCode	SetAnalyticalSolution_TGV3D (SimCtx *simCtx)
	Populate the Eulerian field vectors with the Taylor–Green vortex at the current time.
static PetscErrorCode	SetAnalyticalSolution_ZeroFlow (SimCtx *simCtx)
	Set every analytical Eulerian field to the quiescent zero-flow reference state.
static PetscErrorCode	SetAnalyticalSolution_UniformFlow (SimCtx *simCtx)
	Fill the analytical Eulerian fields with the configured spatially uniform flow state.
static PetscErrorCode	SetAnalyticalSolutionForParticles_TGV3D (Vec tempVec, SimCtx *simCtx)
	Evaluate Taylor–Green velocities at swarm particle positions and store them in tempVec.
static PetscErrorCode	SetAnalyticalSolutionForParticles_UniformFlow (Vec tempVec, SimCtx *simCtx)
	Fill the particle velocity vector with the configured uniform analytical flow.
static PetscErrorCode	EvaluateConfiguredScalarProfile (const VerificationScalarConfig *cfg, PetscReal x, PetscReal y, PetscReal z, PetscReal *value)
	Internal helper that evaluates the configured scalar verification profile.
PetscBool	AnalyticalTypeRequiresCustomGeometry (const char *analytical_type)
	Implementation of AnalyticalTypeRequiresCustomGeometry().
PetscBool	AnalyticalTypeSupportsInterpolationError (const char *analytical_type)
	Implementation of AnalyticalTypeSupportsInterpolationError().
PetscErrorCode	SetAnalyticalGridInfo (UserCtx *user)
	Internal helper implementation: SetAnalyticalGridInfo().
PetscErrorCode	AnalyticalSolutionEngine (SimCtx *simCtx)
	Implementation of AnalyticalSolutionEngine().
PetscErrorCode	SetAnalyticalSolutionForParticles (Vec tempVec, SimCtx *simCtx)
	Implementation of SetAnalyticalSolutionForParticles().
PetscErrorCode	EvaluateAnalyticalScalarProfile (const SimCtx *simCtx, PetscReal x, PetscReal y, PetscReal z, PetscReal t, PetscReal *value)
	Implementation of EvaluateAnalyticalScalarProfile().
PetscErrorCode	SetAnalyticalScalarFieldOnParticles (UserCtx *user, ParticleFieldId particle_field_id)
	Implementation of SetAnalyticalScalarFieldOnParticles().
PetscErrorCode	SetAnalyticalScalarFieldAtCellCenters (UserCtx *user, Vec targetVec)
	Implementation of SetAnalyticalScalarFieldAtCellCenters().

Variables
static const FieldId	kAnalyticalCellFields [] = {FIELD_ID_UCAT, FIELD_ID_P}
	Cell-centred fields every analytical setter finalizes, in the order the post-projection finalization uses for the solved flow.

Detailed Description

Implements the analytical solution engine for initializing or driving the simulation.

This file provides a modular and extensible framework for applying analytical solutions to the Eulerian fields. The primary entry point is AnalyticalSolutionEngine, which acts as a dispatcher based on user configuration.

— DESIGN PHILOSOPHY —

1. Non-Dimensional Core: All calculations within this engine are performed in non-dimensional units (e.g., reference velocity U_ref=1.0, reference length L_ref=1.0). This is critical for consistency with the core numerical solver, which also operates on non-dimensional equations. The simCtx->scaling parameters are intentionally NOT used here; they are reserved for dimensionalization during I/O and post-processing only.
2. Separation of Concerns: The role of this engine is to set the physical state of the fluid at a given time t. This involves:
   • For most types (TGV3D, ZERO_FLOW): setting Ucat and P directly.
   • For curvilinear-aware types (UNIFORM_FLOW): setting Ucont via metric dot products and deriving Ucat through Contra2Cart.
   • Declaring the physical values on the boundaries by populating the boundary condition vector (user->Bcs.Ubcs). It does NOT fill dummy cells itself. After setting the physical state, every setter runs the same sequence the solved flow is finalized with (see FinalizePostProjectionCellFields): UpdateDummyCells on non-periodic faces, SynchronizePeriodicCellFields on periodic ones, UpdateCornerNodes for edges and corners, and the periodic synchronization again because corner averaging can overwrite periodic endpoints. The boundary values are this module's own Ubcs, which is why the setters do not call that finalizer: it would refresh Ubcs from the boundary handlers.
3. Extensibility: The dispatcher design makes it straightforward to add new analytical solutions. A developer only needs to add a new else if condition and a corresponding static implementation function, without modifying any other part of the solver.

Definition in file AnalyticalSolutions.c.

Macro Definition Documentation

__FUNCT__ [1/11]

#define __FUNCT__   "SetAnalyticalGridInfo"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [2/11]

#define __FUNCT__   "AnalyticalSolutionEngine"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [3/11]

#define __FUNCT__   "SetAnalyticalSolution_TGV3D"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [4/11]

#define __FUNCT__   "SetAnalyticalSolution_ZeroFlow"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [5/11]

#define __FUNCT__   "SetAnalyticalSolution_UniformFlow"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [6/11]

#define __FUNCT__   "SetAnalyticalSolutionForParticles_TGV3D"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [7/11]

#define __FUNCT__   "SetAnalyticalSolutionForParticles_UniformFlow"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [8/11]

#define __FUNCT__   "SetAnalyticalSolutionForParticles"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [9/11]

#define __FUNCT__   "EvaluateAnalyticalScalarProfile"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [10/11]

#define __FUNCT__   "SetAnalyticalScalarFieldOnParticles"

Definition at line 83 of file AnalyticalSolutions.c.

__FUNCT__ [11/11]

#define __FUNCT__   "SetAnalyticalScalarFieldAtCellCenters"

Definition at line 83 of file AnalyticalSolutions.c.

Function Documentation

SetAnalyticalSolution_TGV3D()

static PetscErrorCode SetAnalyticalSolution_TGV3D ( SimCtx *  simCtx )

static

Populate the Eulerian field vectors with the Taylor–Green vortex at the current time.

Definition at line 246 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    UserCtx       *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    
    // --- NON-DIMENSIONAL TGV Parameters ---
    const PetscReal V0  = 1.0; // Non-dimensional reference velocity.
    const PetscReal rho = 1.0; // Non-dimensional reference density.
    const PetscReal p0  = 0.0; // Non-dimensional reference pressure.
    
    // Kinematic viscosity is derived from the non-dimensional Reynolds number.
    const PetscReal nu  = (simCtx->ren > 0) ? (1.0 / simCtx->ren) : 0.0; 
 
    const PetscReal k   = 1.0; // Wavenumber, assumes a non-dimensional [0, 2*pi] domain.
    const PetscReal t   = simCtx->ti;
 
    LOG_ALLOW(GLOBAL,LOG_TRACE,"TGV Setup: t = %.4f, V0* = %.4f, rho* = %.4f, k = %.4f, p0* = %4.f, nu = %.6f.\n",simCtx->ti,V0,rho,k,p0,nu);
 
    const PetscReal vel_decay = exp(-2.0 * nu * k * k * t);
    const PetscReal prs_decay = exp(-4.0 * nu * k * k * t);
 
    PetscFunctionBeginUser;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        UserCtx* user = &user_finest[bi];
        DMDALocalInfo info = user->info;
        PetscInt      xs = info.xs, xe = info.xs + info.xm;
        PetscInt      ys = info.ys, ye = info.ys + info.ym;
        PetscInt      zs = info.zs, ze = info.zs + info.zm;
        PetscInt      mx = info.mx, my = info.my, mz = info.mz;
 
        Cmpnts      ***ucat, ***ubcs;
        const Cmpnts ***cent, ***cent_x, ***cent_y, ***cent_z;
        PetscReal   ***p;
        
        // Define loop bounds for physical interior cells owned by this rank.
        PetscInt lxs = (xs == 0) ? xs + 1 : xs, lxe = (xe == mx) ? xe - 1 : xe;
        PetscInt lys = (ys == 0) ? ys + 1 : ys, lye = (ye == my) ? ye - 1 : ye;
        PetscInt lzs = (zs == 0) ? zs + 1 : zs, lze = (ze == mz) ? ze - 1 : ze;
 
        // --- Get Arrays ---
        ierr = DMDAVecGetArray(user->fda, user->Ucat, &ucat); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->da, user->P, &p); CHKERRQ(ierr);
        ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->Cent, &cent); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lCentx, &cent_x); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lCenty, &cent_y); CHKERRQ(ierr);
        ierr = DMDAVecGetArrayRead(user->fda, user->lCentz, &cent_z); CHKERRQ(ierr);
 
        // --- Set INTERIOR cell-centered velocity (Ucat) ---
        for (PetscInt k_cell = lzs; k_cell < lze; k_cell++) {
            for (PetscInt j_cell = lys; j_cell < lye; j_cell++) {
                for (PetscInt i_cell = lxs; i_cell < lxe; i_cell++) {
                    const PetscReal cx = cent[k_cell][j_cell][i_cell].x, cy = cent[k_cell][j_cell][i_cell].y, cz = cent[k_cell][j_cell][i_cell].z;
                    ucat[k_cell][j_cell][i_cell].x =  V0 * sin(k*cx) * cos(k*cy) * cos(k*cz) * vel_decay;
                    ucat[k_cell][j_cell][i_cell].y = -V0 * cos(k*cx) * sin(k*cy) * cos(k*cz) * vel_decay;
                    ucat[k_cell][j_cell][i_cell].z =  0.0;
                }
            }
        }
 
 
        // --- Set INTERIOR cell-centered pressure (P) ---
        for (PetscInt k_cell = lzs; k_cell < lze; k_cell++) {
            for (PetscInt j_cell = lys; j_cell < lye; j_cell++) {
                for (PetscInt i_cell = lxs; i_cell < lxe; i_cell++) {
                    const PetscReal cx = cent[k_cell][j_cell][i_cell].x, cy = cent[k_cell][j_cell][i_cell].y;
                    p[k_cell][j_cell][i_cell] = p0 + (rho * V0 * V0 / 4.0) * (cos(2*k*cx) + cos(2*k*cy)) * prs_decay;
                }
            }
        }
        
 
        // --- Set BOUNDARY condition vector for velocity (Ubcs) ---
        if (xs == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) {
            const PetscReal fcx=cent_x[k][j][xs].x, fcy=cent_x[k][j][xs].y, fcz=cent_x[k][j][xs].z;
            ubcs[k][j][xs].x = V0*sin(k*fcx)*cos(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][j][xs].y = -V0*cos(k*fcx)*sin(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][j][xs].z = 0.0;
        }
        if (xe == mx) for (PetscInt k=zs; k<ze; k++) for (PetscInt j=ys; j<ye; j++) {
            const PetscReal fcx=cent_x[k][j][xe-1].x, fcy=cent_x[k][j][xe-1].y, fcz=cent_x[k][j][xe-1].z;
            ubcs[k][j][xe-1].x = V0*sin(k*fcx)*cos(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][j][xe-1].y = -V0*cos(k*fcx)*sin(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][j][xe-1].z = 0.0;
        }
        if (ys == 0) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) {
            const PetscReal fcx=cent_y[k][ys][i].x, fcy=cent_y[k][ys][i].y, fcz=cent_y[k][ys][i].z;
            ubcs[k][ys][i].x = V0*sin(k*fcx)*cos(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][ys][i].y = -V0*cos(k*fcx)*sin(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][ys][i].z = 0.0;
        }
        if (ye == my) for (PetscInt k=zs; k<ze; k++) for (PetscInt i=xs; i<xe; i++) {
            const PetscReal fcx=cent_y[k][ye-1][i].x, fcy=cent_y[k][ye-1][i].y, fcz=cent_y[k][ye-1][i].z;
            ubcs[k][ye-1][i].x = V0*sin(k*fcx)*cos(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][ye-1][i].y = -V0*cos(k*fcx)*sin(k*fcy)*cos(k*fcz)*vel_decay; ubcs[k][ye-1][i].z = 0.0;
        }
        if (zs == 0) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) {
            const PetscReal fcx=cent_z[zs][j][i].x, fcy=cent_z[zs][j][i].y, fcz=cent_z[zs][j][i].z;
            ubcs[zs][j][i].x = V0*sin(k*fcx)*cos(k*fcy)*cos(k*fcz)*vel_decay; ubcs[zs][j][i].y = -V0*cos(k*fcx)*sin(k*fcy)*cos(k*fcz)*vel_decay; ubcs[zs][j][i].z = 0.0;
        }
        if (ze == mz) for (PetscInt j=ys; j<ye; j++) for (PetscInt i=xs; i<xe; i++) {
            const PetscReal fcx=cent_z[ze-1][j][i].x, fcy=cent_z[ze-1][j][i].y, fcz=cent_z[ze-1][j][i].z;
            ubcs[ze-1][j][i].x = V0*sin(k*fcx)*cos(k*fcy)*cos(k*fcz)*vel_decay; ubcs[ze-1][j][i].y = -V0*cos(k*fcx)*sin(k*fcy)*cos(k*fcz)*vel_decay; ubcs[ze-1][j][i].z = 0.0;
        }
 
        // --- Restore all arrays ---
        ierr = DMDAVecRestoreArray(user->fda, user->Ucat, &ucat); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->da, user->P, &p); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->Cent, &cent); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentx, &cent_x); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lCenty, &cent_y); CHKERRQ(ierr);
        ierr = DMDAVecRestoreArrayRead(user->fda, user->lCentz, &cent_z); CHKERRQ(ierr);
 
        // Pre-Dummy cell update synchronization.
        ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT);
        ierr = UpdateLocalGhosts(user, FIELD_ID_P);
 
        // --- Finalize all ghost cell values ---
        ierr = UpdateDummyCells(user, FIELD_ID_UCAT); CHKERRQ(ierr);
        ierr = UpdateDummyCells(user, FIELD_ID_P); CHKERRQ(ierr);
        ierr = SynchronizePeriodicCellFields(user, 2, kAnalyticalCellFields); CHKERRQ(ierr);
        ierr = UpdateCornerNodes(user, FIELD_ID_UCAT); CHKERRQ(ierr);
        ierr = UpdateCornerNodes(user, FIELD_ID_P); CHKERRQ(ierr);
        ierr = SynchronizePeriodicCellFields(user, 2, kAnalyticalCellFields); CHKERRQ(ierr);
 
        // Final Synchronization.
        ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT);
        ierr = UpdateLocalGhosts(user, FIELD_ID_P);
 
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetAnalyticalSolution_ZeroFlow()

static PetscErrorCode SetAnalyticalSolution_ZeroFlow ( SimCtx *  simCtx )

static

Set every analytical Eulerian field to the quiescent zero-flow reference state.

Definition at line 380 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    PetscFunctionBeginUser;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        UserCtx *user = &user_finest[bi];
 
        ierr = VecZeroEntries(user->Ucat);     CHKERRQ(ierr);
        ierr = VecZeroEntries(user->P);        CHKERRQ(ierr);
        ierr = VecZeroEntries(user->Bcs.Ubcs); CHKERRQ(ierr);
 
        // Ghost-cell finalization — identical sequence to TGV3D
        ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_P);    CHKERRQ(ierr);
        ierr = UpdateDummyCells(user, FIELD_ID_UCAT);          CHKERRQ(ierr);
        ierr = UpdateDummyCells(user, FIELD_ID_P);          CHKERRQ(ierr);
        ierr = SynchronizePeriodicCellFields(user, 2, kAnalyticalCellFields); CHKERRQ(ierr);
        ierr = UpdateCornerNodes(user, FIELD_ID_UCAT);         CHKERRQ(ierr);
        ierr = UpdateCornerNodes(user, FIELD_ID_P);         CHKERRQ(ierr);
        ierr = SynchronizePeriodicCellFields(user, 2, kAnalyticalCellFields); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_P);    CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetAnalyticalSolution_UniformFlow()

static PetscErrorCode SetAnalyticalSolution_UniformFlow ( SimCtx *  simCtx )

static

Fill the analytical Eulerian fields with the configured spatially uniform flow state.

Definition at line 415 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    const Cmpnts uniform_velocity = simCtx->AnalyticalUniformVelocity;
    const PetscReal u = uniform_velocity.x;
    const PetscReal v = uniform_velocity.y;
    const PetscReal w = uniform_velocity.z;
 
    PetscFunctionBeginUser;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        UserCtx *user = &user_finest[bi];
        DMDALocalInfo info = user->info;
        PetscInt      xs = info.xs, xe = info.xs + info.xm;
        PetscInt      ys = info.ys, ye = info.ys + info.ym;
        PetscInt      zs = info.zs, ze = info.zs + info.zm;
        PetscInt      mx = info.mx, my = info.my, mz = info.mz;
 
        // --- Step 1: Set Ucont (contravariant flux) via Cart2Contra ---
        ierr = UniformCart2Contra(user, u, v, w); CHKERRQ(ierr);
 
        // --- Step 2: Set Ubcs at boundaries (physical Cartesian velocity) ---
        Cmpnts ***ubcs;
        ierr = DMDAVecGetArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
        if (xs == 0) for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) ubcs[k][j][xs] = uniform_velocity;
        if (xe == mx) for (PetscInt k = zs; k < ze; k++) for (PetscInt j = ys; j < ye; j++) ubcs[k][j][xe - 1] = uniform_velocity;
        if (ys == 0) for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) ubcs[k][ys][i] = uniform_velocity;
        if (ye == my) for (PetscInt k = zs; k < ze; k++) for (PetscInt i = xs; i < xe; i++) ubcs[k][ye - 1][i] = uniform_velocity;
        if (zs == 0) for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) ubcs[zs][j][i] = uniform_velocity;
        if (ze == mz) for (PetscInt j = ys; j < ye; j++) for (PetscInt i = xs; i < xe; i++) ubcs[ze - 1][j][i] = uniform_velocity;
        ierr = DMDAVecRestoreArray(user->fda, user->Bcs.Ubcs, &ubcs); CHKERRQ(ierr);
 
        // --- Step 3: Zero pressure ---
        ierr = VecZeroEntries(user->P); CHKERRQ(ierr);
 
        // --- Step 4: Finalize state — derive Ucat from Ucont via metric inversion ---
        const FieldId staggered_fields[] = {FIELD_ID_UCONT};
        ierr = SynchronizePeriodicStaggeredFields(user, 1, staggered_fields); CHKERRQ(ierr);
        ierr = Contra2Cart(user);                CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT);  CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_P);     CHKERRQ(ierr);
        ierr = UpdateDummyCells(user, FIELD_ID_UCAT);           CHKERRQ(ierr);
        ierr = UpdateDummyCells(user, FIELD_ID_P);           CHKERRQ(ierr);
        ierr = SynchronizePeriodicCellFields(user, 2, kAnalyticalCellFields); CHKERRQ(ierr);
        ierr = UpdateCornerNodes(user, FIELD_ID_UCAT);          CHKERRQ(ierr);
        ierr = UpdateCornerNodes(user, FIELD_ID_P);          CHKERRQ(ierr);
        ierr = SynchronizePeriodicCellFields(user, 2, kAnalyticalCellFields); CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_UCAT);  CHKERRQ(ierr);
        ierr = UpdateLocalGhosts(user, FIELD_ID_P);     CHKERRQ(ierr);
    }
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetAnalyticalSolutionForParticles_TGV3D()

static PetscErrorCode SetAnalyticalSolutionForParticles_TGV3D ( Vec  tempVec, SimCtx *  simCtx  )

static

Evaluate Taylor–Green velocities at swarm particle positions and store them in tempVec.

Definition at line 475 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    PetscInt nLocal;
    PetscReal *data;
    
    PetscFunctionBeginUser;
 
    // TGV3D parameters (matching your Eulerian implementation)
    const PetscReal V0  = 1.0;
    const PetscReal k   = 1.0;
    const PetscReal nu  = (simCtx->ren > 0) ? (1.0 / simCtx->ren) : 0.0;
    const PetscReal t   = simCtx->ti;
    const PetscReal vel_decay = exp(-2.0 * nu * k * k * t);
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "TGV3D Particles: t=%.4f, V0=%.4f, k=%.4f, nu=%.6f\n", t, V0, k, nu);
    
    ierr = VecGetLocalSize(tempVec, &nLocal); CHKERRQ(ierr);
    ierr = VecGetArray(tempVec, &data); CHKERRQ(ierr);
    
    // Process particles: data is interleaved [x0,y0,z0, x1,y1,z1, ...]
    for (PetscInt i = 0; i < nLocal; i += 3) {
        const PetscReal x = data[i];
        const PetscReal y = data[i+1];
        const PetscReal z = data[i+2];
        
        // TGV3D velocity field
        data[i]   =  V0 * sin(k*x) * cos(k*y) * cos(k*z) * vel_decay;  // u
        data[i+1] = -V0 * cos(k*x) * sin(k*y) * cos(k*z) * vel_decay;  // v
        data[i+2] =  0.0;                                                // w
    }
    
    ierr = VecRestoreArray(tempVec, &data); CHKERRQ(ierr);
    
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

SetAnalyticalSolutionForParticles_UniformFlow()

static PetscErrorCode SetAnalyticalSolutionForParticles_UniformFlow ( Vec  tempVec, SimCtx *  simCtx  )

static

Fill the particle velocity vector with the configured uniform analytical flow.

Definition at line 517 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    PetscInt nLocal;
    PetscReal *data;
    const Cmpnts uniform_velocity = simCtx->AnalyticalUniformVelocity;
 
    PetscFunctionBeginUser;
 
    ierr = VecGetLocalSize(tempVec, &nLocal); CHKERRQ(ierr);
    ierr = VecGetArray(tempVec, &data); CHKERRQ(ierr);
    for (PetscInt i = 0; i < nLocal; i += 3) {
        data[i] = uniform_velocity.x;
        data[i + 1] = uniform_velocity.y;
        data[i + 2] = uniform_velocity.z;
    }
    ierr = VecRestoreArray(tempVec, &data); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

EvaluateConfiguredScalarProfile()

static PetscErrorCode EvaluateConfiguredScalarProfile ( const VerificationScalarConfig *  cfg, PetscReal  x, PetscReal  y, PetscReal  z, PetscReal *  value  )

static

Internal helper that evaluates the configured scalar verification profile.

Local to this translation unit.

Definition at line 574 of file AnalyticalSolutions.c.

{
    PetscFunctionBeginUser;
    if (!cfg) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "VerificationScalarConfig cannot be NULL.");
    if (!value) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Scalar output pointer cannot be NULL.");
 
    if (strcmp(cfg->profile, "CONSTANT") == 0) {
        *value = cfg->value;
    } else if (strcmp(cfg->profile, "LINEAR_X") == 0) {
        *value = cfg->phi0 + cfg->slope_x * x;
    } else if (strcmp(cfg->profile, "SIN_PRODUCT") == 0) {
        *value = cfg->amplitude *
                 PetscSinReal(cfg->kx * x) *
                 PetscSinReal(cfg->ky * y) *
                 PetscSinReal(cfg->kz * z);
    } else {
        SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
                "Unsupported verification scalar profile '%s'.", cfg->profile);
    }
 
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

AnalyticalTypeRequiresCustomGeometry()

PetscBool AnalyticalTypeRequiresCustomGeometry

(

const char *

analytical_type

)

Implementation of AnalyticalTypeRequiresCustomGeometry().

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

)

Implementation of AnalyticalTypeSupportsInterpolationError().

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:

SetAnalyticalGridInfo()

PetscErrorCode SetAnalyticalGridInfo

(

UserCtx *

user

)

Internal helper implementation: SetAnalyticalGridInfo().

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:

AnalyticalSolutionEngine()

PetscErrorCode AnalyticalSolutionEngine

(

SimCtx *

simCtx

)

Implementation of AnalyticalSolutionEngine().

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
	)

Implementation of SetAnalyticalSolutionForParticles().

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
	)

Implementation of EvaluateAnalyticalScalarProfile().

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
	)

Implementation of SetAnalyticalScalarFieldOnParticles().

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
	)

Implementation of SetAnalyticalScalarFieldAtCellCenters().

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:

Variable Documentation

kAnalyticalCellFields

const FieldId kAnalyticalCellFields[] = {FIELD_ID_UCAT, FIELD_ID_P}

static

Cell-centred fields every analytical setter finalizes, in the order the post-projection finalization uses for the solved flow.

Definition at line 43 of file AnalyticalSolutions.c.

{FIELD_ID_UCAT, FIELD_ID_P};