solvers.c

Go to the documentation of this file.

#include "solvers.h" // The new header we will create
 
#undef __FUNCT__
#define __FUNCT__ "FlowSolver"
/**
 * @brief Implementation of \ref FlowSolver().
 * @details Full API contract (arguments, ownership, side effects) is documented with
 *          the header declaration in `include/solvers.h`.
 * @see FlowSolver()
 */
PetscErrorCode FlowSolver(SimCtx *simCtx)
{
    PetscErrorCode ierr;
    UserMG         *usermg = NULL;
    PetscInt       level;
    UserCtx        *user = NULL;
 
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    if (simCtx->mom_solver_type != MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK &&
        simCtx->mom_solver_type != MOMENTUM_SOLVER_EXPLICIT_RK &&
        simCtx->mom_solver_type != MOMENTUM_SOLVER_NEWTON_KRYLOV) {
        SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
                "Unknown momentum solver type %d. Supported values are EXPLICIT_RK, DUALTIME_PICARD_JAMESON_RK, and NEWTON_KRYLOV.",
                simCtx->mom_solver_type);
    }
 
    usermg = &simCtx->usermg;
    level  = usermg->mglevels - 1;
    user   = usermg->mgctx[level].user;
    LOG_ALLOW(GLOBAL, LOG_INFO, "[Step %d] Entering orchestrator...\n", simCtx->step);
 
    /*
    // ========================================================================
    //   SECTION: O-Grid Specific Force Calculations (Legacy Feature)
    // ========================================================================
    // This was a specialized calculation for non-immersed O-grid cases.
    if (simCtx->Ogrid && !simCtx->immersed) {
        LOG_ALLOW(GLOBAL, LOG_DEBUG, "Calculating O-grid forces...\n");
        // Calc_forces_Ogrid(&user[0], simCtx->step, 0);
    }
    */
 
 
    
    // ========================================================================
    //   SECTION: Turbulence Models (LES)
    // ========================================================================
    // These models compute the turbulent eddy viscosity (Nu_t) which is then
    // used by the momentum solver in the diffusion term.
 
 
    if (simCtx->les) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Updating LES subgrid-scale model...\n");
        for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
            // Strain rates are formed from the Cartesian velocity, so it must be
            // reconstructed and its ghost and periodic images refreshed first.
            const FieldId staggered_fields[] = {FIELD_ID_UCONT};
            const FieldId cell_fields[]      = {FIELD_ID_UCAT};
 
            ierr = SynchronizePeriodicStaggeredFields(&user[bi], 1, staggered_fields); CHKERRQ(ierr);
            ierr = Contra2Cart(&user[bi]); CHKERRQ(ierr);
            ierr = SynchronizePeriodicCellFields(&user[bi], 1, cell_fields); CHKERRQ(ierr);
            ierr = UpdateLocalGhosts(&user[bi], FIELD_ID_UCAT); CHKERRQ(ierr);
 
            // Contra2Cart rebuilds the Cartesian field from the contravariant one and so
            // discards any near-wall correction the previous momentum solve applied. The
            // wall model is re-applied here, before the strain rates are formed, so the
            // closure and the momentum equation agree about the near-wall velocity. On
            // the coarse grids a wall model exists to serve, the uncorrected strain is
            // under-predicted, and the eddy viscosity with it. Mirrors the order the
            // boundary pass itself uses: reconstruct, exchange, correct, exchange.
            if (simCtx->wallfunction) {
                ierr = ApplyWallFunction(&user[bi]); CHKERRQ(ierr);
                ierr = SynchronizePeriodicCellFields(&user[bi], 1, cell_fields); CHKERRQ(ierr);
                ierr = UpdateLocalGhosts(&user[bi], FIELD_ID_UCAT); CHKERRQ(ierr);
            }
 
            // Only the dynamic model has a coefficient to recompute, and it honours
            // its own cadence. The constant model's coefficient never changes.
            if (simCtx->les == DYNAMIC_SMAGORINSKY &&
                simCtx->step % simCtx->les_config.dynamic_frequency == 0) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "  Computing dynamic coefficient for block %d.\n", bi);
                ierr = ComputeSmagorinskyConstant(&user[bi]); CHKERRQ(ierr);
            }
 
            ierr = ComputeEddyViscosityLES(&user[bi]); CHKERRQ(ierr);
            ierr = LogLESDiagnostics(&user[bi]); CHKERRQ(ierr);
        }
    }
    
 
    // ========================================================================
    //   SECTION: Momentum Equation Solver
    // ========================================================================
    // This is the core of the time step. It computes an intermediate velocity
    // field by solving the momentum equations.
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Beginning momentum step solve (Solver = %s)...\n", MomentumSolverTypeToString(simCtx->mom_solver_type));
    
    // Since IBM is disabled, we pass NULL for ibm and fsi arguments.
    // ierr = ImpRK(user, NULL, NULL); CHKERRQ(ierr);
    // Add new momentum solver types here only after wiring the enum, parser, docs, and tests.
    if(simCtx->mom_solver_type == MOMENTUM_SOLVER_DUALTIME_PICARD_JAMESON_RK) {
        ierr = MomentumSolver_DualTime_Picard_JamesonRK(user,NULL,NULL); CHKERRQ(ierr);
    } else if(simCtx->mom_solver_type == MOMENTUM_SOLVER_EXPLICIT_RK) {   
        // Since IBM is disabled, we pass NULL for ibm and fsi arguments.
        ierr = MomentumSolver_Explicit_RungeKutta4(user, NULL, NULL); CHKERRQ(ierr);
    } else if (simCtx->mom_solver_type == MOMENTUM_SOLVER_NEWTON_KRYLOV) {
        ierr = MomentumSolver_NewtonKrylov(user, NULL, NULL); CHKERRQ(ierr);
    }
// ========================================================================
//   SECTION: Pressure-Poisson Solver
// ========================================================================
// This step enforces the continuity equation (incompressibility) by solving
// for a pressure correction field.
    
 LOG_ALLOW(GLOBAL, LOG_INFO, "Beginning pressure-Poisson solve (Poisson Flag = %d)...\n", simCtx->poisson);
    
 if (simCtx->poisson == 0) {
   ierr = PoissonSolver_Multigrid(usermg); CHKERRQ(ierr);
 } else {
   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG,
           "Unsupported Poisson solver type %d. The current runtime supports only the multigrid path (poisson = 0).",
           simCtx->poisson);
 }
    
    // ========================================================================
    //   SECTION: Velocity Correction (Projection)
    // ========================================================================
    // The pressure correction is used to update the pressure field and project
    // the intermediate velocity onto a divergence-free space.
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "Applying velocity correction/projection step...\n");
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        ierr = UpdatePressure(&user[bi]); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL,LOG_INFO," Pressure Updated for Block %d.\n",bi);
    
        ierr = ProjectVelocity(&user[bi]); CHKERRQ(ierr);
    
        LOG_ALLOW(GLOBAL,LOG_INFO," Velocity corrected for Block %d.\n",bi);
 
    // Ensure local ghost cells for the final pressure field are correct
    ierr = UpdateLocalGhosts(&user[bi], FIELD_ID_P);
    }
    
    // ========================================================================
    
    // ========================================================================
    //   SECTION: Final Diagnostics and I/O
    // ========================================================================
    
    for (PetscInt bi = 0; bi < simCtx->block_number; bi++) {
        LOG_ALLOW(GLOBAL, LOG_INFO, "Finalizing state & Diagnostics for block %d...\n", bi);
        
        // --- Perform Divergence Check ---
        // This is a diagnostic to verify the quality of the velocity correction.
    ierr = ComputeDivergence(&user[bi]); CHKERRQ(ierr);
 
    // -- Log Continuity metrics ----
    ierr = LOG_CONTINUITY_METRICS(&user[bi]);
 
        /* Reported here rather than beside the LES block, because a wall model is
           configured independently of LES and its last pass this step is the boundary
           pass, not the closure prologue. A no-op when no wall model is active. */
        ierr = LogWallModelDiagnostics(&user[bi]); CHKERRQ(ierr);
        /* The force a driven-periodic controller applied this step. A no-op without one. */
        ierr = LogDrivenFlowDiagnostics(&user[bi]); CHKERRQ(ierr);
        /*
        // --- Immersed Boundary Interpolation (Post-Correction) ---
        // This step would update the velocity values AT the IB nodes to match the
        // newly corrected fluid field. Important for the next time step.
        if (simCtx->immersed) {
            for (PetscInt ibi = 0; ibi < simCtx->NumberOfBodies; ibi++) {
                ibm_interpolation_advanced(&user[bi], &simCtx->ibm[ibi], ibi, 1);
            }
        }
        */
 
        // }
    }
    
    LOG_ALLOW(GLOBAL, LOG_INFO, "orchestrator finished for step %d.\n", simCtx->step);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}