1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Time-stepping functions for Navier-Stokes example using PETSc 10 11 #include <ceed.h> 12 #include <petscdmplex.h> 13 #include <petscts.h> 14 15 #include "../navierstokes.h" 16 #include "../qfunctions/newtonian_state.h" 17 18 // Insert Boundary values if it's a new time 19 PetscErrorCode UpdateBoundaryValues(User user, Vec Q_loc, PetscReal t) { 20 PetscFunctionBeginUser; 21 if (user->time_bc_set != t) { 22 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t, NULL, NULL, NULL)); 23 user->time_bc_set = t; 24 } 25 PetscFunctionReturn(PETSC_SUCCESS); 26 } 27 28 // RHS (Explicit time-stepper) function setup 29 // This is the RHS of the ODE, given as u_t = G(t,u) 30 // This function takes in a state vector Q and writes into G 31 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data) { 32 User user = *(User *)user_data; 33 Ceed ceed = user->ceed; 34 PetscScalar dt; 35 Vec Q_loc = user->Q_loc; 36 PetscMemType q_mem_type; 37 38 PetscFunctionBeginUser; 39 // Update time dependent data 40 PetscCall(UpdateBoundaryValues(user, Q_loc, t)); 41 if (user->phys->solution_time_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_rhs_ctx->op, user->phys->solution_time_label, &t)); 42 PetscCall(TSGetTimeStep(ts, &dt)); 43 if (user->phys->timestep_size_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_rhs_ctx->op, user->phys->timestep_size_label, &dt)); 44 45 PetscCall(ApplyCeedOperatorGlobalToGlobal(Q, G, user->op_rhs_ctx)); 46 47 PetscCall(VecReadPetscToCeed(Q_loc, &q_mem_type, user->q_ceed)); 48 49 // Inverse of the mass matrix 50 PetscCall(KSPSolve(user->mass_ksp, G, G)); 51 52 PetscCall(VecReadCeedToPetsc(user->q_ceed, q_mem_type, Q_loc)); 53 PetscFunctionReturn(PETSC_SUCCESS); 54 } 55 56 // Surface forces function setup 57 static PetscErrorCode Surface_Forces_NS(DM dm, Vec G_loc, PetscInt num_walls, const PetscInt walls[], PetscScalar *reaction_force) { 58 DMLabel face_label; 59 const PetscScalar *g; 60 PetscInt dof, dim = 3; 61 MPI_Comm comm; 62 PetscSection s; 63 64 PetscFunctionBeginUser; 65 PetscCall(PetscArrayzero(reaction_force, num_walls * dim)); 66 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 67 PetscCall(DMGetLabel(dm, "Face Sets", &face_label)); 68 PetscCall(VecGetArrayRead(G_loc, &g)); 69 for (PetscInt w = 0; w < num_walls; w++) { 70 const PetscInt wall = walls[w]; 71 IS wall_is; 72 PetscCall(DMGetLocalSection(dm, &s)); 73 PetscCall(DMLabelGetStratumIS(face_label, wall, &wall_is)); 74 if (wall_is) { // There exist such points on this process 75 PetscInt num_points; 76 PetscInt num_comp = 0; 77 const PetscInt *points; 78 PetscCall(PetscSectionGetFieldComponents(s, 0, &num_comp)); 79 PetscCall(ISGetSize(wall_is, &num_points)); 80 PetscCall(ISGetIndices(wall_is, &points)); 81 for (PetscInt i = 0; i < num_points; i++) { 82 const PetscInt p = points[i]; 83 const StateConservative *r; 84 PetscCall(DMPlexPointLocalRead(dm, p, g, &r)); 85 PetscCall(PetscSectionGetDof(s, p, &dof)); 86 for (PetscInt node = 0; node < dof / num_comp; node++) { 87 for (PetscInt j = 0; j < 3; j++) { 88 reaction_force[w * dim + j] -= r[node].momentum[j]; 89 } 90 } 91 } 92 PetscCall(ISRestoreIndices(wall_is, &points)); 93 } 94 PetscCall(ISDestroy(&wall_is)); 95 } 96 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, reaction_force, dim * num_walls, MPIU_SCALAR, MPI_SUM, comm)); 97 // Restore Vectors 98 PetscCall(VecRestoreArrayRead(G_loc, &g)); 99 PetscFunctionReturn(PETSC_SUCCESS); 100 } 101 102 // Implicit time-stepper function setup 103 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, void *user_data) { 104 User user = *(User *)user_data; 105 Ceed ceed = user->ceed; 106 PetscScalar dt; 107 Vec Q_loc = user->Q_loc, Q_dot_loc = user->Q_dot_loc, G_loc; 108 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type; 109 110 PetscFunctionBeginUser; 111 // Get local vectors 112 PetscCall(DMGetNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 113 114 // Update time dependent data 115 PetscCall(UpdateBoundaryValues(user, Q_loc, t)); 116 if (user->phys->solution_time_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_ifunction, user->phys->solution_time_label, &t)); 117 PetscCall(TSGetTimeStep(ts, &dt)); 118 if (user->phys->timestep_size_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(user->op_ifunction, user->phys->timestep_size_label, &dt)); 119 120 // Global-to-local 121 PetscCall(DMGlobalToLocalBegin(user->dm, Q, INSERT_VALUES, Q_loc)); 122 PetscCall(DMGlobalToLocalBegin(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc)); 123 PetscCall(DMGlobalToLocalEnd(user->dm, Q, INSERT_VALUES, Q_loc)); 124 PetscCall(DMGlobalToLocalEnd(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc)); 125 126 // Place PETSc vectors in CEED vectors 127 PetscCall(VecReadPetscToCeed(Q_loc, &q_mem_type, user->q_ceed)); 128 PetscCall(VecReadPetscToCeed(Q_dot_loc, &q_dot_mem_type, user->q_dot_ceed)); 129 PetscCall(VecPetscToCeed(G_loc, &g_mem_type, user->g_ceed)); 130 131 // Apply CEED operator 132 PetscCall(PetscLogEventBegin(FLUIDS_CeedOperatorApply, Q, G, 0, 0)); 133 PetscCall(PetscLogGpuTimeBegin()); 134 PetscCallCeed(user->ceed, CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, CEED_REQUEST_IMMEDIATE)); 135 PetscCall(PetscLogGpuTimeEnd()); 136 PetscCall(PetscLogEventEnd(FLUIDS_CeedOperatorApply, Q, G, 0, 0)); 137 138 // Restore vectors 139 PetscCall(VecReadCeedToPetsc(user->q_ceed, q_mem_type, Q_loc)); 140 PetscCall(VecReadCeedToPetsc(user->q_dot_ceed, q_dot_mem_type, Q_dot_loc)); 141 PetscCall(VecCeedToPetsc(user->g_ceed, g_mem_type, G_loc)); 142 143 if (user->app_ctx->sgs_model_type == SGS_MODEL_DATA_DRIVEN) { 144 PetscCall(SgsDDApplyIFunction(user, Q_loc, G_loc)); 145 } 146 147 // Local-to-Global 148 PetscCall(VecZeroEntries(G)); 149 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G)); 150 151 // Restore vectors 152 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) { 157 User user = *(User *)user_data; 158 double shift_double = shift; 159 PetscBool J_is_matceed, J_is_mffd, J_pre_is_matceed, J_pre_is_mffd; 160 161 PetscFunctionBeginUser; 162 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd)); 163 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATCEED, &J_is_matceed)); 164 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATMFFD, &J_pre_is_mffd)); 165 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATCEED, &J_pre_is_matceed)); 166 167 PetscCall(MatCeedSetContextDouble(user->mat_ijacobian, "ijacobian time shift", shift_double)); 168 169 if (J_is_matceed || J_is_mffd) { 170 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 171 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 172 } else PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J)); 173 174 if (J_pre_is_matceed && J != J_pre) { 175 PetscCall(MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY)); 176 PetscCall(MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY)); 177 } else if (!J_pre_is_matceed && !J_pre_is_mffd && J != J_pre) { 178 PetscCall(MatCeedAssembleCOO(user->mat_ijacobian, J_pre)); 179 } 180 PetscFunctionReturn(PETSC_SUCCESS); 181 } 182 183 PetscErrorCode WriteOutput(User user, Vec Q, PetscInt step_no, PetscScalar time) { 184 Vec Q_loc; 185 char file_path[PETSC_MAX_PATH_LEN]; 186 PetscViewer viewer; 187 188 PetscFunctionBeginUser; 189 if (user->app_ctx->checkpoint_vtk) { 190 // Set up output 191 PetscCall(DMGetLocalVector(user->dm, &Q_loc)); 192 PetscCall(PetscObjectSetName((PetscObject)Q_loc, "StateVec")); 193 PetscCall(VecZeroEntries(Q_loc)); 194 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 195 196 // Output 197 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no)); 198 199 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, FILE_MODE_WRITE, &viewer)); 200 PetscCall(VecView(Q_loc, viewer)); 201 PetscCall(PetscViewerDestroy(&viewer)); 202 if (user->dm_viz) { 203 Vec Q_refined, Q_refined_loc; 204 char file_path_refined[PETSC_MAX_PATH_LEN]; 205 PetscViewer viewer_refined; 206 207 PetscCall(DMGetGlobalVector(user->dm_viz, &Q_refined)); 208 PetscCall(DMGetLocalVector(user->dm_viz, &Q_refined_loc)); 209 PetscCall(PetscObjectSetName((PetscObject)Q_refined_loc, "Refined")); 210 211 PetscCall(MatInterpolate(user->interp_viz, Q, Q_refined)); 212 PetscCall(VecZeroEntries(Q_refined_loc)); 213 PetscCall(DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc)); 214 215 PetscCall( 216 PetscSNPrintf(file_path_refined, sizeof file_path_refined, "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no)); 217 218 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), file_path_refined, FILE_MODE_WRITE, &viewer_refined)); 219 PetscCall(VecView(Q_refined_loc, viewer_refined)); 220 PetscCall(DMRestoreLocalVector(user->dm_viz, &Q_refined_loc)); 221 PetscCall(DMRestoreGlobalVector(user->dm_viz, &Q_refined)); 222 PetscCall(PetscViewerDestroy(&viewer_refined)); 223 } 224 PetscCall(DMRestoreLocalVector(user->dm, &Q_loc)); 225 } 226 227 // Save data in a binary file for continuation of simulations 228 if (user->app_ctx->add_stepnum2bin) { 229 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution-%" PetscInt_FMT ".bin", user->app_ctx->output_dir, step_no)); 230 } else { 231 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", user->app_ctx->output_dir)); 232 } 233 PetscCall(PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer)); 234 235 PetscInt32 token = PetscDefined(USE_64BIT_INDICES) ? FLUIDS_FILE_TOKEN_64 : FLUIDS_FILE_TOKEN_32; 236 PetscCall(PetscViewerBinaryWrite(viewer, &token, 1, PETSC_INT32)); 237 PetscCall(PetscViewerBinaryWrite(viewer, &step_no, 1, PETSC_INT)); 238 time /= user->units->second; // Dimensionalize time back 239 PetscCall(PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL)); 240 PetscCall(VecView(Q, viewer)); 241 PetscCall(PetscViewerDestroy(&viewer)); 242 PetscFunctionReturn(PETSC_SUCCESS); 243 } 244 245 // CSV Monitor 246 PetscErrorCode TSMonitor_WallForce(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) { 247 User user = ctx; 248 Vec G_loc; 249 PetscInt num_wall = user->app_ctx->wall_forces.num_wall, dim = 3; 250 const PetscInt *walls = user->app_ctx->wall_forces.walls; 251 PetscViewer viewer = user->app_ctx->wall_forces.viewer; 252 PetscViewerFormat format = user->app_ctx->wall_forces.viewer_format; 253 PetscScalar *reaction_force; 254 PetscBool iascii; 255 256 PetscFunctionBeginUser; 257 if (!viewer) PetscFunctionReturn(PETSC_SUCCESS); 258 PetscCall(DMGetNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 259 PetscCall(PetscMalloc1(num_wall * dim, &reaction_force)); 260 PetscCall(Surface_Forces_NS(user->dm, G_loc, num_wall, walls, reaction_force)); 261 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 262 263 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 264 265 if (iascii) { 266 if (format == PETSC_VIEWER_ASCII_CSV && !user->app_ctx->wall_forces.header_written) { 267 PetscCall(PetscViewerASCIIPrintf(viewer, "Step,Time,Wall,ForceX,ForceY,ForceZ\n")); 268 user->app_ctx->wall_forces.header_written = PETSC_TRUE; 269 } 270 for (PetscInt w = 0; w < num_wall; w++) { 271 PetscInt wall = walls[w]; 272 if (format == PETSC_VIEWER_ASCII_CSV) { 273 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ",%g,%" PetscInt_FMT ",%g,%g,%g\n", step_no, time, wall, 274 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2])); 275 276 } else { 277 PetscCall(PetscViewerASCIIPrintf(viewer, "Wall %" PetscInt_FMT " Forces: Force_x = %12g, Force_y = %12g, Force_z = %12g\n", wall, 278 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2])); 279 } 280 } 281 } 282 PetscCall(PetscFree(reaction_force)); 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285 286 // User provided TS Monitor 287 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) { 288 User user = ctx; 289 290 PetscFunctionBeginUser; 291 // Print every 'checkpoint_interval' steps 292 if (user->app_ctx->checkpoint_interval <= 0 || step_no % user->app_ctx->checkpoint_interval != 0 || 293 (user->app_ctx->cont_steps == step_no && step_no != 0)) { 294 PetscFunctionReturn(PETSC_SUCCESS); 295 } 296 297 PetscCall(WriteOutput(user, Q, step_no, time)); 298 PetscFunctionReturn(PETSC_SUCCESS); 299 } 300 301 // TS: Create, setup, and solve 302 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, Vec *Q, PetscScalar *f_time, TS *ts) { 303 MPI_Comm comm = user->comm; 304 TSAdapt adapt; 305 PetscScalar final_time; 306 307 PetscFunctionBeginUser; 308 PetscCall(TSCreate(comm, ts)); 309 PetscCall(TSSetDM(*ts, dm)); 310 PetscCall(TSSetApplicationContext(*ts, user)); 311 if (phys->implicit) { 312 PetscCall(TSSetType(*ts, TSBDF)); 313 if (user->op_ifunction) { 314 PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &user)); 315 } else { // Implicit integrators can fall back to using an RHSFunction 316 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user)); 317 } 318 if (user->mat_ijacobian) { 319 PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &user)); 320 } 321 } else { 322 PetscCheck(user->op_rhs_ctx, comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction"); 323 PetscCall(TSSetType(*ts, TSRK)); 324 PetscCall(TSRKSetType(*ts, TSRK5F)); 325 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user)); 326 } 327 PetscCall(TSSetMaxTime(*ts, 500. * user->units->second)); 328 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER)); 329 if (app_ctx->test_type == TESTTYPE_NONE) PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE)); 330 PetscCall(TSSetTimeStep(*ts, 1.e-2 * user->units->second)); 331 PetscCall(TSGetAdapt(*ts, &adapt)); 332 PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * user->units->second, 1.e2 * user->units->second)); 333 PetscCall(TSSetFromOptions(*ts)); 334 if (user->mat_ijacobian) { 335 if (app_ctx->amat_type && !strcmp(app_ctx->amat_type, MATSHELL)) { 336 SNES snes; 337 KSP ksp; 338 Mat Pmat, Amat; 339 340 PetscCall(TSGetSNES(*ts, &snes)); 341 PetscCall(SNESGetKSP(snes, &ksp)); 342 PetscCall(CreateSolveOperatorsFromMatCeed(ksp, user->mat_ijacobian, PETSC_FALSE, &Amat, &Pmat)); 343 PetscCall(TSSetIJacobian(*ts, user->mat_ijacobian, Pmat, NULL, NULL)); 344 PetscCall(MatDestroy(&Amat)); 345 PetscCall(MatDestroy(&Pmat)); 346 } 347 } 348 user->time_bc_set = -1.0; // require all BCs be updated 349 if (app_ctx->cont_steps) { // continue from previous timestep data 350 PetscInt count; 351 PetscViewer viewer; 352 353 if (app_ctx->cont_time <= 0) { // Legacy files did not include step number and time 354 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_time_file, FILE_MODE_READ, &viewer)); 355 PetscCall(PetscViewerBinaryRead(viewer, &app_ctx->cont_time, 1, &count, PETSC_REAL)); 356 PetscCall(PetscViewerDestroy(&viewer)); 357 PetscCheck(app_ctx->cont_steps != -1, comm, PETSC_ERR_ARG_INCOMP, 358 "-continue step number not specified, but checkpoint file does not contain a step number (likely written by older code version)"); 359 } 360 PetscCall(TSSetTime(*ts, app_ctx->cont_time * user->units->second)); 361 PetscCall(TSSetStepNumber(*ts, app_ctx->cont_steps)); 362 } 363 if (app_ctx->test_type == TESTTYPE_NONE) { 364 PetscCall(TSMonitorSet(*ts, TSMonitor_NS, user, NULL)); 365 } 366 if (app_ctx->wall_forces.viewer) { 367 PetscCall(TSMonitorSet(*ts, TSMonitor_WallForce, user, NULL)); 368 } 369 if (app_ctx->turb_spanstats_enable) { 370 PetscCall(TSMonitorSet(*ts, TSMonitor_TurbulenceStatistics, user, NULL)); 371 CeedScalar previous_time = app_ctx->cont_time * user->units->second; 372 PetscCallCeed(user->ceed, 373 CeedOperatorSetContextDouble(user->spanstats.op_stats_collect_ctx->op, user->spanstats.previous_time_label, &previous_time)); 374 } 375 if (app_ctx->diff_filter_monitor) PetscCall(TSMonitorSet(*ts, TSMonitor_DifferentialFilter, user, NULL)); 376 377 if (app_ctx->sgs_train_enable) { 378 PetscCall(TSMonitorSet(*ts, TSMonitor_SGS_DD_Training, user, NULL)); 379 PetscCall(TSSetPostStep(*ts, TSPostStep_SGS_DD_Training)); 380 } 381 // Solve 382 PetscReal start_time; 383 PetscInt start_step; 384 PetscCall(TSGetTime(*ts, &start_time)); 385 PetscCall(TSGetStepNumber(*ts, &start_step)); 386 387 PetscCall(PetscLogDefaultBegin()); // So we can use PetscLogStageGetPerfInfo without -log_view 388 PetscPreLoadBegin(PETSC_FALSE, "Fluids Solve"); 389 PetscCall(TSSetTime(*ts, start_time)); 390 PetscCall(TSSetStepNumber(*ts, start_step)); 391 if (PetscPreLoadingOn) { 392 // LCOV_EXCL_START 393 SNES snes; 394 Vec Q_preload; 395 PetscReal rtol; 396 PetscCall(VecDuplicate(*Q, &Q_preload)); 397 PetscCall(VecCopy(*Q, Q_preload)); 398 PetscCall(TSGetSNES(*ts, &snes)); 399 PetscCall(SNESGetTolerances(snes, NULL, &rtol, NULL, NULL, NULL)); 400 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, .99, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 401 PetscCall(TSSetSolution(*ts, Q_preload)); 402 PetscCall(TSStep(*ts)); 403 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 404 PetscCall(VecDestroy(&Q_preload)); 405 // LCOV_EXCL_STOP 406 } else { 407 PetscCall(PetscBarrier((PetscObject)*ts)); 408 PetscCall(TSSolve(*ts, *Q)); 409 } 410 PetscPreLoadEnd(); 411 412 PetscCall(TSGetSolveTime(*ts, &final_time)); 413 *f_time = final_time; 414 415 if (app_ctx->test_type == TESTTYPE_NONE) { 416 PetscInt step_no; 417 PetscCall(TSGetStepNumber(*ts, &step_no)); 418 if (user->app_ctx->checkpoint_interval > 0 || user->app_ctx->checkpoint_interval == -1) { 419 PetscCall(WriteOutput(user, *Q, step_no, final_time)); 420 } 421 422 PetscLogStage stage_id; 423 PetscEventPerfInfo stage_perf; 424 425 PetscCall(PetscLogStageGetId("Fluids Solve", &stage_id)); 426 PetscCall(PetscLogStageGetPerfInfo(stage_id, &stage_perf)); 427 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time taken for solution (sec): %g\n", stage_perf.time)); 428 } 429 PetscFunctionReturn(PETSC_SUCCESS); 430 } 431