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.
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.

Macro Definition Documentation

__FUNCT__

#define __FUNCT__   "AnalyticalSolutionEngine"

Definition at line 61 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.

Definition at line 70 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)) {
        SETERRQ1(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) {
                SETERRQ1(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 {
        SETERRQ1(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.

Definition at line 38 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:

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.

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

Parameters

[in]

simCtx

The main simulation context, containing all configuration and state.

Returns

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

Definition at line 179 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);
    }
    /*
     * --- 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.
        SETERRQ1(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.

Definition at line 485 of file AnalyticalSolutions.c.

{
    PetscErrorCode ierr;
    PetscInt nLocal;
    PetscReal *vels;
 
    PetscFunctionBeginUser;
    
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Type: %s\n", 
              simCtx->AnalyticalSolutionType ? simCtx->AnalyticalSolutionType : "default");
    
    // Check for specific analytical solution types
    if (simCtx->AnalyticalSolutionType && strcmp(simCtx->AnalyticalSolutionType, "TGV3D") == 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Using TGV3D solution.\n");
        ierr = SetAnalyticalSolutionForParticles_TGV3D(tempVec, simCtx); CHKERRQ(ierr);
        return 0;
    }
        
    ierr = VecRestoreArray(tempVec, &vels); CHKERRQ(ierr);
    
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: