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