initialcondition.h File Reference

#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscsys.h>
#include <petscdmcomposite.h>
#include <petscsystypes.h>
#include "variables.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "interpolation.h"
#include "AnalyticalSolutions.h"
#include "ParticleMotion.h"
#include "Boundaries.h"
#include "runloop.h"

Include dependency graph for initialcondition.h:

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

Go to the source code of this file.

Functions
PetscErrorCode	SetInitialInteriorField (UserCtx *user, FieldId field_id)
	Sets the initial values for the INTERIOR of a specified Eulerian field.
PetscErrorCode	PopulateInitialUcont (UserCtx *user)
	Populate Ucont for one fresh-start block from the configured IC mode.
PetscErrorCode	InitializeEulerianState (SimCtx *simCtx)
	High-level orchestrator to set the complete initial state of the Eulerian solver.

Function Documentation

SetInitialInteriorField()

PetscErrorCode SetInitialInteriorField	(	UserCtx *	user,
		FieldId	field_id
	)

Sets the initial values for the INTERIOR of a specified Eulerian field.

This function initializes the interior nodes of Ucont based on the mode selected by simCtx->initialConditionMode.

Supported profiles for "Ucont":

• IC_MODE_ZERO: All interior contravariant components are set to zero.
• IC_MODE_CONSTANT_CARTESIAN: UniformCart2Contra dots the Cartesian vector (InitialConstantContra.x/y/z) with the local metric vectors to fill all three contravariant components correctly across the entire interior.
• IC_MODE_CONSTANT_STREAMWISE: Sets only the contravariant component along the streamwise axis (from flowDirection or the identified INLET face) proportional to icVelocityPhysical * |A_n|.
• IC_MODE_POISEUILLE: Separable parabolic profile in the two cross-stream index directions; centerline speed is icVelocityPhysical; streamwise axis from flowDirection or the identified INLET face.

Parameters

user	The main UserCtx struct, containing all simulation data and configuration.
field_id	Typed identity of the field to initialize.

Returns

PetscErrorCode 0 on success.

Sets the initial values for the INTERIOR of a specified Eulerian field.

Local to this translation unit.

Definition at line 14 of file initialcondition.c.

{
    PetscErrorCode ierr;
    const char *fieldName = FieldCanonicalName(field_id);
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    SimCtx *simCtx = user->simCtx;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Setting initial INTERIOR field for '%s' with mode %d.\n", fieldName, simCtx->initialConditionMode);
 
    // This function currently only implements logic for Ucont.
    if (field_id != FIELD_ID_UCONT) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Skipping SetInitialInteriorField for non-Ucont field '%s'.\n", fieldName);
 
        PROFILE_FUNCTION_END;
 
        PetscFunctionReturn(0);
    }
 
    // --- 1. Get DMDA info and grid dimensions ---
    DMDALocalInfo info;
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    const PetscInt im_phys = info.mx - 1;
    const PetscInt jm_phys = info.my - 1;
    const PetscInt km_phys = info.mz - 1;
 
    const PetscReal u_cart = simCtx->InitialConstantContra.x;
    const PetscReal v_cart = simCtx->InitialConstantContra.y;
    const PetscReal w_cart = simCtx->InitialConstantContra.z;
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "IC cartesian=(%.3f,%.3f,%.3f) ic_velocity_physical=%.3f mode=%d\n",
              (double)u_cart, (double)v_cart, (double)w_cart,
              (double)simCtx->icVelocityPhysical, (int)simCtx->initialConditionMode);
 
    // --- 2. Early dispatch: cartesian Constant delegates to the uniform converter ---
    if (simCtx->initialConditionMode == IC_MODE_CONSTANT_CARTESIAN) {
        ierr = UniformCart2Contra(user, u_cart, v_cart, w_cart); CHKERRQ(ierr);
        PROFILE_FUNCTION_END;
        PetscFunctionReturn(0);
    }
 
    // --- 3. Resolve flow direction for streamwise Constant and Poiseuille ---
    const PetscBool needs_flow_dir = (PetscBool)(
        simCtx->initialConditionMode == IC_MODE_POISEUILLE ||
        simCtx->initialConditionMode == IC_MODE_CONSTANT_STREAMWISE);
    FlowDirection fd          = FLOW_DIR_UNSET;
    PetscInt      flow_axis   = 0;
    PetscReal     flow_dir_sign = 1.0;
 
    if (needs_flow_dir) {
        if (user->inletFaceDefined)
            fd = (FlowDirection)user->identifiedInletBCFace;
        else if (simCtx->flowDirection != FLOW_DIR_UNSET)
            fd = simCtx->flowDirection;
        else
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER,
                    "Streamwise Constant and Poiseuille IC modes require either an INLET face or -flow_direction.");
        flow_axis     = (PetscInt)fd / 2;
        flow_dir_sign = ((PetscInt)fd % 2 == 0) ? 1.0 : -1.0;
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "IC flow_direction=%d (axis=%d sign=%.1f)\n",
                  (int)fd, (int)flow_axis, (double)flow_dir_sign);
    }
 
    // --- 4. Open arrays for non-cartesian modes ---
    Cmpnts ***csi_arr, ***eta_arr, ***zet_arr;
    ierr = DMDAVecGetArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
 
    Cmpnts ***ucont_arr;
    ierr = DMDAVecGetArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
 
    PetscInt i, j, k;
    const PetscInt xs = info.xs, xe = info.xs + info.xm;
    const PetscInt ys = info.ys, ye = info.ys + info.ym;
    const PetscInt zs = info.zs, ze = info.zs + info.zm;
        
    for (k = zs; k < ze; k++) {
        for (j = ys; j < ye; j++) {
            for (i = xs; i < xe; i++) {
 
                // Check to ensure we only set initial conditions for PHYSICAL cells, not ghost cells.
                // Ghost cells (at indices 0 and n) will be set later by ApplyBoundaryConditions.
                //
                // Grid structure: For n physical grid points, DMDA has size n+1
                //   - im_phys = mx - 1 = n (number of coordinate points, also equals number of cells + 1)
                //   - Physical cell indices: [1, im_phys-1] = [1, n-1] (gives n-1 physical cells)
                //   - Ghost cells at boundaries: index 0 and index im_phys (= n)
                //
                // Example: n=25 physical points → im_phys=25
                //   - Physical cells: indices 1..24 (24 cells)
                //   - Ghost cells: indices 0 and 25
                const PetscBool is_interior = (i > 0 && i < im_phys &&
                                               j > 0 && j < jm_phys &&
                                               k > 0 && k < km_phys);
 
                if (is_interior) {
                    Cmpnts ucont_val = {0.0, 0.0, 0.0}; // Default to zero velocity
                    PetscReal normal_velocity_mag = 0.0;
 
                    switch (simCtx->initialConditionMode) {
                        case IC_MODE_ZERO:
                            break;
                        case IC_MODE_CONSTANT_STREAMWISE:
                            normal_velocity_mag = simCtx->icVelocityPhysical;
                            break;
                        case IC_MODE_POISEUILLE:
                            {
                                PetscInt  cs1, cs2, n1, n2;
                                PetscBool per1, per2;
                                if (flow_axis == 0)      { cs1 = j; cs2 = k; n1 = jm_phys; n2 = km_phys;
                                                           per1 = (PetscBool)(simCtx->j_periodic != 0); per2 = (PetscBool)(simCtx->k_periodic != 0); }
                                else if (flow_axis == 1) { cs1 = i; cs2 = k; n1 = im_phys; n2 = km_phys;
                                                           per1 = (PetscBool)(simCtx->i_periodic != 0); per2 = (PetscBool)(simCtx->k_periodic != 0); }
                                else                     { cs1 = i; cs2 = j; n1 = im_phys; n2 = jm_phys;
                                                           per1 = (PetscBool)(simCtx->i_periodic != 0); per2 = (PetscBool)(simCtx->j_periodic != 0); }
                                /* An axis of n nodes has n-1 cells, and cell c is centred at logical
                                   position c - 1/2, so the walls sit at 1/2 and n - 1/2 in that measure.
                                   The parabola vanishes there, on the walls, not at the first cell
                                   centres. A periodic cross-stream axis has no walls to vanish on, so
                                   the profile is uniform along it: a streamwise-periodic, spanwise-
                                   periodic channel then starts from its exact parabola. */
                                const PetscReal half1 = 0.5 * (PetscReal)(n1 - 1);
                                const PetscReal half2 = 0.5 * (PetscReal)(n2 - 1);
                                const PetscReal n1_norm = ((PetscReal)cs1 - 0.5 - half1) / half1;
                                const PetscReal n2_norm = ((PetscReal)cs2 - 0.5 - half2) / half2;
                                const PetscReal f1 = per1 ? 1.0 : (1.0 - n1_norm * n1_norm);
                                const PetscReal f2 = per2 ? 1.0 : (1.0 - n2_norm * n2_norm);
                                normal_velocity_mag = simCtx->icVelocityPhysical * f1 * f2;
                                if (normal_velocity_mag < 0.0) normal_velocity_mag = 0.0;
                            }
                            break;
                        default:
                            LOG_ALLOW(LOCAL, LOG_WARNING, "Unrecognized initial-condition mode %d. Defaulting to zero.\n", simCtx->initialConditionMode);
                            break;
                    }
 
                    // Step B: apply flow direction and set the single contravariant flux component.
                    if (normal_velocity_mag != 0.0) {
                        const PetscReal signed_vel = normal_velocity_mag * flow_dir_sign * user->GridOrientation;
                        if (flow_axis == 0) {
                            const PetscReal area = sqrt(csi_arr[k][j][i].x * csi_arr[k][j][i].x +
                                                        csi_arr[k][j][i].y * csi_arr[k][j][i].y +
                                                        csi_arr[k][j][i].z * csi_arr[k][j][i].z);
                            ucont_val.x = signed_vel * area;
                        } else if (flow_axis == 1) {
                            const PetscReal area = sqrt(eta_arr[k][j][i].x * eta_arr[k][j][i].x +
                                                        eta_arr[k][j][i].y * eta_arr[k][j][i].y +
                                                        eta_arr[k][j][i].z * eta_arr[k][j][i].z);
                            ucont_val.y = signed_vel * area;
                        } else {
                            const PetscReal area = sqrt(zet_arr[k][j][i].x * zet_arr[k][j][i].x +
                                                        zet_arr[k][j][i].y * zet_arr[k][j][i].y +
                                                        zet_arr[k][j][i].z * zet_arr[k][j][i].z);
                            ucont_val.z = signed_vel * area;
                        }
                    }
                    ucont_arr[k][j][i] = ucont_val;
                } // end if(is_interior)
            }
        }
    }
    ierr = DMDAVecRestoreArray(user->fda, user->Ucont, &ucont_arr); CHKERRQ(ierr);
 
    // --- 5. Restore arrays ---
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lCsi, &csi_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lEta, &eta_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayRead(user->fda, user->lZet, &zet_arr); CHKERRQ(ierr);
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PopulateInitialUcont()

PetscErrorCode PopulateInitialUcont

(

UserCtx *

user

)

Populate Ucont for one fresh-start block from the configured IC mode.

Built-in modes generate Ucont directly. File mode loads either Ucat or Ucont from the staged IC directory and converts Ucat when necessary.

Parameters

[in,out]

user

Block context whose velocity field is populated.

Returns

PETSc error code.

Populate Ucont for one fresh-start block from the configured IC mode.

Definition at line 231 of file initialcondition.c.

{
    PetscErrorCode ierr;
    SimCtx *simCtx = user->simCtx;
 
    PetscFunctionBeginUser;
    if (simCtx->initialConditionMode == IC_MODE_FILE) {
        ierr = LoadInitialUcont(user); CHKERRQ(ierr);
    } else {
        ierr = SetInitialInteriorField(user, FIELD_ID_UCONT); CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InitializeEulerianState()

PetscErrorCode InitializeEulerianState

(

SimCtx *

simCtx

)

High-level orchestrator to set the complete initial state of the Eulerian solver.

This function is called once from main() before the time loop begins. It inspects the simulation context to determine whether to perform a fresh start (t=0) or restart from saved files. It then delegates to the appropriate helper function. Finally, it initializes the solver's history vectors (Ucont_o, P_o, etc.) to ensure the first time step has the necessary data.

Parameters

simCtx

Simulation context controlling the operation.

Returns

PetscErrorCode 0 on success.

High-level orchestrator to set the complete initial state of the Eulerian solver.

Local to this translation unit.

Definition at line 369 of file initialcondition.c.

{
    PetscErrorCode ierr;
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Initializing Eulerian State ---\n");
 
    if (simCtx->StartStep > 0) {
        if(strcmp(simCtx->eulerianSource,"analytical")==0){
            LOG_ALLOW(GLOBAL,LOG_INFO,"Initializing Analytical Solution type: %s (t=%.4f, step=%d).\n",simCtx->AnalyticalSolutionType,simCtx->StartTime,simCtx->StartStep);
            ierr = AnalyticalSolutionEngine(simCtx);
        }
        else{
            LOG_ALLOW(GLOBAL, LOG_INFO, "Starting from RESTART files (t=%.4f, step=%d).\n",
                    simCtx->StartTime, simCtx->StartStep);
            ierr = SetInitialFluidState_Load(simCtx); CHKERRQ(ierr);
        }
        /* Statistics resume from the same bundle the flow state came from, and are
         * restored regardless of the Eulerian source: an analytical restart still
         * continues a window that was accumulating before it. */
        ierr = RestoreFieldStatisticsState(simCtx, simCtx->StartStep); CHKERRQ(ierr);
    } else { // StartStep = 0
        LOG_ALLOW(GLOBAL, LOG_INFO, "Performing a FRESH START (t=0, step=0).\n");
        if(strcmp(simCtx->eulerianSource,"solve")==0){
            ierr = SetInitialFluidState_FreshStart(simCtx); CHKERRQ(ierr);
        }else if(strcmp(simCtx->eulerianSource,"load")==0){
            LOG_ALLOW(GLOBAL,LOG_INFO,"FRESH START in LOAD mode. Reading files (t=%.4f,step=%d).\n",
                      simCtx->StartTime,simCtx->StartStep);
            ierr=SetInitialFluidState_Load(simCtx);CHKERRQ(ierr);
        }else if(strcmp(simCtx->eulerianSource,"analytical")==0){
            LOG_ALLOW(GLOBAL,LOG_INFO,"FRESH START in ANALYTICAL mode. Initializing Analytical Solution type: %s (t=%.4f,step=%d).\n",
                      simCtx->AnalyticalSolutionType,simCtx->StartTime,simCtx->StartStep);
            ierr=AnalyticalSolutionEngine(simCtx);CHKERRQ(ierr);
        }
    }
 
    // This crucial step, taken from the end of the legacy setup, ensures
    // that the history vectors (Ucont_o, Ucont_rm1, etc.) are correctly
    // populated before the first call to the time-stepping loop.
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        ierr = UpdateSolverHistoryVectors(&user_finest[bi],
                                           (PetscBool)(simCtx->StartStep > 0 && simCtx->restartHistoryAvailable)); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Eulerian State Initialized and History Vectors Populated ---\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: