initialcondition.c

Go to the documentation of this file.

/**
 * @file initialcondition.c  //  Setup the Initial conditions for different cases. 
 * @brief Test program for DMSwarm interpolation using the fdf-curvIB method.
 **/
 
 #include "initialcondition.h"
 
#undef __FUNCT__
#define __FUNCT__ "SetInitialInteriorField"
/**
 * @brief Internal helper implementation: `SetInitialInteriorField()`.
 * @details Local to this translation unit.
 */
PetscErrorCode SetInitialInteriorField(UserCtx *user, const char *fieldName)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    SimCtx *simCtx = user->simCtx;
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Setting initial INTERIOR field for '%s' with profile %d.\n", fieldName, simCtx->FieldInitialization);
 
    // This function currently only implements logic for Ucont.
    if (strcmp(fieldName, "Ucont") != 0) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Skipping SetInitialInteriorField for non-Ucont field '%s'.\n", fieldName);
 
        PROFILE_FUNCTION_END;
 
        PetscFunctionReturn(0);
    }
 
    // --- 1. Get references to required data and PETSc arrays ---
    DM            fieldDM = user->fda;
    Vec           fieldVec = user->Ucont;
    DMDALocalInfo info;
    ierr = DMDAGetLocalInfo(fieldDM, &info); CHKERRQ(ierr);
 
    Vec      localCoor;
    Cmpnts ***coor_arr;
    ierr = DMGetCoordinatesLocal(user->da, &localCoor); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, localCoor, &coor_arr); CHKERRQ(ierr);
    
    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);
   
    const PetscInt im_phys = info.mx - 1;
    const PetscInt jm_phys = info.my - 1;
    const PetscInt km_phys = info.mz - 1;    
 
    // --- 2. Compute Cell-Center Coordinates (only if needed by the selected profile) ---
    Cmpnts ***cent_coor = NULL;
    PetscInt xs_cell=0, xm_cell=0, ys_cell=0, ym_cell=0, zs_cell=0, zm_cell=0;
    
    if (simCtx->FieldInitialization == 2) { // Profile 2 (Poiseuille) requires cell centers.
        
        ierr = GetOwnedCellRange(&info, 0, &xs_cell, &xm_cell); CHKERRQ(ierr);
        ierr = GetOwnedCellRange(&info, 1, &ys_cell, &ym_cell); CHKERRQ(ierr);
        ierr = GetOwnedCellRange(&info, 2, &zs_cell, &zm_cell); CHKERRQ(ierr);
      
        if (xm_cell > 0 && ym_cell > 0 && zm_cell > 0) {
            ierr = Allocate3DArray(&cent_coor, zm_cell, ym_cell, xm_cell); CHKERRQ(ierr);
            ierr = InterpolateFieldFromCornerToCenter(coor_arr, cent_coor, user); CHKERRQ(ierr);
            LOG_ALLOW(LOCAL, LOG_DEBUG, "Computed temporary cell-center coordinates for Poiseuille profile.\n");
        }
    }
    
    // --- 3. Loop Over Owned Nodes and Apply Initial Condition to Interior ---
    Cmpnts ***ucont_arr;
    ierr = DMDAVecGetArray(fieldDM, fieldVec, &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;
    
    const Cmpnts uin = simCtx->InitialConstantContra; // Max/average velocity from user options
    // Flow into a negative face (e.g., -Zeta at k=0) is in the positive physical direction (+z).
    // Flow into a positive face (e.g., +Zeta at k=mz-1) is in the negative physical direction (-z).
 
    LOG_ALLOW(GLOBAL,LOG_DEBUG,"Initial Constant Flux (Contravariant) = (%.3f, %.3f, %.3f)\n",(double)uin.x, (double)uin.y, (double)uin.z);
    
    const PetscReal flow_direction_sign = (user->identifiedInletBCFace % 2 == 0) ? 1.0 : -1.0;
        
    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;
 
                    // Step A: Determine the magnitude of the desired physical normal velocity.
                    switch (simCtx->FieldInitialization) {
                        case 0: // Zero initial velocity
                            normal_velocity_mag = 0.0;
                            break;
                case 1: /* Constant Normal Velocity */
                    if (user->identifiedInletBCFace == BC_FACE_NEG_X ||
                user->identifiedInletBCFace == BC_FACE_POS_X) {
                normal_velocity_mag = uin.x;
                    } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
               normal_velocity_mag = uin.y;
                    } else {
               normal_velocity_mag = uin.z;
                    }
              break; 
                        case 2: // Poiseuille Normal Velocity
                            {
                                                  // This profile assumes flow is aligned with the k-index direction.
                                // It uses grid indices (i,j) to define the cross-section, which works for bent geometries.
                                PetscReal u0 = 0.0;
                                if (user->identifiedInletBCFace <= BC_FACE_POS_X) u0 = uin.x;
                                else if (user->identifiedInletBCFace <= BC_FACE_POS_Y) u0 = uin.y;
                                else u0 = uin.z; // Assumes Z-like inlet direction
 
                                // Define channel geometry in "index space" based on global grid dimensions
                                // We subtract 2.0 because the interior runs from index 1 to mx-2 (or my-2).
                                const PetscReal i_width  = (PetscReal)(im_phys - 2);
                                const PetscReal j_width  = (PetscReal)(jm_phys - 2);
                                const PetscReal i_center = 1.0 + i_width / 2.0;
                                const PetscReal j_center = 1.0 + j_width / 2.0;
 
                                // Create normalized coordinates for the current point (i,j), ranging from -1 to 1
                                const PetscReal i_norm = (i - i_center) / (i_width / 2.0);
                                const PetscReal j_norm = (j - j_center) / (j_width / 2.0);
 
                                // Apply the parabolic profile for a rectangular/square channel
                                // V(i,j) = V_max * (1 - i_norm^2) * (1 - j_norm^2)
                                const PetscReal profile_i = 1.0 - i_norm * i_norm;
                                const PetscReal profile_j = 1.0 - j_norm * j_norm;
                                normal_velocity_mag = u0 * profile_i * profile_j;
 
                                // Clamp to zero for any points outside the channel (or due to minor float errors)
                                if (normal_velocity_mag < 0.0) {
                                    normal_velocity_mag = 0.0;
                                }
                            }
                            break;
                /*
                                PetscReal r_sq = 0.0;
                                const PetscInt i_local = i - xs_cell, j_local = j - ys_cell, k_local = k - zs_cell;
                                if (cent_coor && i_local >= 0 && i_local < xm_cell && j_local >= 0 && j_local < ym_cell && k_local >= 0 && k_local < zm_cell) {
                                    const Cmpnts* center = &cent_coor[k_local][j_local][i_local];
                                    if (user->identifiedInletBCFace <= BC_FACE_POS_X) r_sq = center->y * center->y + center->z * center->z;
                                    else if (user->identifiedInletBCFace <= BC_FACE_POS_Y) r_sq = center->x * center->x + center->z * center->z;
                                    else r_sq = center->x * center->x + center->y * center->y;
                    // pick the correct contravariant component for center‐line speed 
                    PetscReal u0;
                    if (user->identifiedInletBCFace == BC_FACE_NEG_X ||
                    user->identifiedInletBCFace == BC_FACE_POS_X) {
                      u0 = uin.x;
                    } else if (user->identifiedInletBCFace == BC_FACE_NEG_Y ||
                           user->identifiedInletBCFace == BC_FACE_POS_Y) {
                      u0 = uin.y;
                    } else {
                      u0 = uin.z;
                    }
                    // now form the parabolic profile as before 
                                    normal_velocity_mag = 2.0 * u0 * (1.0 - 4.0 * r_sq);
                                }
                            }
                            break;
                */
                        default:
                            LOG_ALLOW(LOCAL, LOG_WARNING, "Unrecognized FieldInitialization profile %d. Defaulting to zero.\n", simCtx->FieldInitialization);
                            normal_velocity_mag = 0.0;
                            break;
                    }
 
                    // Step B: Apply direction sign and set the correct contravariant component.
                    // The contravariant component U^n = v_n * Area_n, where v_n is the physical normal velocity.
                    if (normal_velocity_mag != 0.0) {
               const PetscReal signed_normal_vel = normal_velocity_mag * flow_direction_sign*user->GridOrientation;
                        
                        if (user->identifiedInletBCFace == BC_FACE_NEG_X || user->identifiedInletBCFace == BC_FACE_POS_X) {
                            const PetscReal area_i = 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_normal_vel * area_i;
 
                LOG_LOOP_ALLOW(GLOBAL,LOG_VERBOSE,k,50," ucont_val.x = %.6f (signed_normal_vel=%.3f × area=%.4f)\n",ucont_val.x, signed_normal_vel, area_i);
                        } 
                        else if (user->identifiedInletBCFace == BC_FACE_NEG_Y || user->identifiedInletBCFace == BC_FACE_POS_Y) {
                            const PetscReal area_j = 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_normal_vel * area_j;
 
                LOG_LOOP_ALLOW(GLOBAL,LOG_VERBOSE,k,50," ucont_val.y = %.6f (signed_normal_vel=%.3f × area=%.4f)\n",ucont_val.y, signed_normal_vel, area_j);
                        } 
                        else { // Z-inlet
                            const PetscReal area_k = 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_normal_vel * area_k;
 
                LOG_LOOP_ALLOW(GLOBAL,LOG_VERBOSE,k,50," i,j,k,ucont_val.z = %d, %d, %d, %.6f (signed_normal_vel=%.3f × area=%.4f)\n",i,j,k,ucont_val.z, signed_normal_vel, area_k);
                        }
                    }
                    ucont_arr[k][j][i] = ucont_val;
                } // end if(is_interior)
            }
        }
    }
    ierr = DMDAVecRestoreArray(fieldDM, fieldVec, &ucont_arr); CHKERRQ(ierr);
 
    // --- 5. Cleanup: Restore arrays and free temporary memory ---
    ierr = DMDAVecRestoreArrayRead(user->fda, localCoor, &coor_arr); CHKERRQ(ierr);
    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);
    
    if (cent_coor) {
        ierr = Deallocate3DArray(cent_coor, zm_cell, ym_cell); CHKERRQ(ierr);
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "FinalizeBlockState"
/**
 * @brief Internal helper implementation: `FinalizeBlockState()`.
 * @details Local to this translation unit.
 */
static PetscErrorCode FinalizeBlockState(UserCtx *user)
{
    PetscErrorCode ierr;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    // This sequence ensures a fully consistent state for a single block.
    ierr = ApplyBoundaryConditions(user); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL,LOG_TRACE," Boundary condition applied.\n");
    // 2. Sync contravariant velocity field.
    ierr = UpdateLocalGhosts(user, "Ucont"); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL,LOG_TRACE," Ucont field ghosts updated.\n");
    
    // 3. Convert to Cartesian velocity.
    ierr = Contra2Cart(user); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL,LOG_TRACE," Converted Ucont to Ucat.\n");
 
    // 4. Sync the new Cartesian velocity field.
    ierr = UpdateLocalGhosts(user, "Ucat"); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL,LOG_TRACE," Ucat field ghosts updated.\n"); 
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
 
#undef __FUNCT__
#define __FUNCT__ "SetInitialFluidState_FreshStart"
/**
 * @brief Internal helper implementation: `SetInitialFluidState_FreshStart()`.
 * @details Local to this translation unit.
 */
static PetscErrorCode SetInitialFluidState_FreshStart(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->usermg.mglevels - 1].user;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d, Block %d: Setting t=0 state.\n", simCtx->rank, bi);
 
        // 1. Set an initial guess for the INTERIOR of the domain.
        //    Replaces the legacy `if(InitialGuessOne)` block.
 
    LOG_ALLOW(GLOBAL,LOG_TRACE," Initializing Interior Ucont field.\n");
        ierr = SetInitialInteriorField(&user_finest[bi], "Ucont"); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL,LOG_TRACE," Interior Ucont field initialized.\n");
 
    // 2. Apply all boundary conditions, convert to Cartesian, and sync ghosts.
    LOG_ALLOW(GLOBAL,LOG_TRACE," Boundary condition application and state finalization initiated.\n");
    ierr = FinalizeBlockState(&user_finest[bi]); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL,LOG_TRACE," Boundary condition application and state finalization complete.\n");
    }
 
    // If using multiple grid blocks, handle the interface conditions between them.
    if (simCtx->block_number > 1) {
      //   LOG_ALLOW(GLOBAL, LOG_INFO, "Updating multi-block interfaces for t=0.\n");
    //  ierr = Block_Interface_U(user_finest); CHKERRQ(ierr);
        // After interface update, ghost regions might be stale. Refresh them.
        for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
             ierr = UpdateLocalGhosts(&user_finest[bi], "Ucont"); CHKERRQ(ierr);
             ierr = UpdateLocalGhosts(&user_finest[bi], "Ucat"); CHKERRQ(ierr);
        }
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "SetInitialFluidState_Load"
/**
 * @brief Internal helper implementation: `SetInitialFluidState_Load()`.
 * @details Local to this translation unit.
 */
static PetscErrorCode SetInitialFluidState_Load(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    UserCtx *user_finest = simCtx->usermg.mgctx[simCtx->mglevels - 1].user;
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        LOG_ALLOW(LOCAL, LOG_INFO, "Rank %d, Block %d: Reading restart files for step %d.\n",
                  simCtx->rank, bi, simCtx->StartStep);
 
        // ReadSimulationFields handles all file I/O for one block.
        ierr = ReadSimulationFields(&user_finest[bi], simCtx->StartStep); CHKERRQ(ierr);
        
        // Apply Boundary Conditions on Read fields
        ierr = ApplyBoundaryConditions(&user_finest[bi]); CHKERRQ(ierr);
        // After reading from a file, the local ghost regions MUST be updated
        // to ensure consistency across process boundaries for the first time step.
        //ierr = UpdateLocalGhosts(&user_finest[bi], "Ucont"); CHKERRQ(ierr);
        //ierr = UpdateLocalGhosts(&user_finest[bi], "Ucat"); CHKERRQ(ierr);
        //ierr = UpdateLocalGhosts(&user_finest[bi], "P"); CHKERRQ(ierr);
        // ... add ghost updates for any other fields read from file ...
    }
    
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "InitializeEulerianState"
/**
 * @brief Internal helper implementation: `InitializeEulerianState()`.
 * @details Local to this translation unit.
 */
PetscErrorCode InitializeEulerianState(SimCtx *simCtx)
{
    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);
        }
    } 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]); CHKERRQ(ierr);
    }
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "--- Eulerian State Initialized and History Vectors Populated ---\n");
 
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}