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