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 "../navierstokes.h" 12 13 // Compute mass matrix for explicit scheme 14 PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, CeedData ceed_data, Vec M) { 15 Vec M_loc; 16 CeedQFunction qf_mass; 17 CeedOperator op_mass; 18 CeedVector m_ceed, ones_vec; 19 CeedInt num_comp_q, q_data_size; 20 PetscFunctionBeginUser; 21 22 // CEED Restriction 23 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q); 24 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); 25 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &m_ceed, NULL); 26 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &ones_vec, NULL); 27 CeedVectorSetValue(ones_vec, 1.0); 28 29 // CEED QFunction 30 PetscCall(CreateMassQFunction(ceed, num_comp_q, q_data_size, &qf_mass)); 31 32 // CEED Operator 33 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 34 CeedOperatorSetField(op_mass, "u", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 35 CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 36 CeedOperatorSetField(op_mass, "v", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 37 38 // Place PETSc vector in CEED vector 39 CeedScalar *m; 40 PetscMemType m_mem_type; 41 PetscCall(DMGetLocalVector(dm, &M_loc)); 42 PetscCall(VecGetArrayAndMemType(M_loc, (PetscScalar **)&m, &m_mem_type)); 43 CeedVectorSetArray(m_ceed, MemTypeP2C(m_mem_type), CEED_USE_POINTER, m); 44 45 // Apply CEED Operator 46 CeedOperatorApply(op_mass, ones_vec, m_ceed, CEED_REQUEST_IMMEDIATE); 47 48 // Restore vectors 49 CeedVectorTakeArray(m_ceed, MemTypeP2C(m_mem_type), NULL); 50 PetscCall(VecRestoreArrayReadAndMemType(M_loc, (const PetscScalar **)&m)); 51 52 // Local-to-Global 53 PetscCall(VecZeroEntries(M)); 54 PetscCall(DMLocalToGlobal(dm, M_loc, ADD_VALUES, M)); 55 PetscCall(DMRestoreLocalVector(dm, &M_loc)); 56 57 // Invert diagonally lumped mass vector for RHS function 58 PetscCall(VecReciprocal(M)); 59 60 // Cleanup 61 CeedVectorDestroy(&ones_vec); 62 CeedVectorDestroy(&m_ceed); 63 CeedQFunctionDestroy(&qf_mass); 64 CeedOperatorDestroy(&op_mass); 65 66 PetscFunctionReturn(0); 67 } 68 69 // Insert Boundary values if it's a new time 70 PetscErrorCode UpdateBoundaryValues(User user, Vec Q_loc, PetscReal t) { 71 PetscFunctionBeginUser; 72 if (user->time_bc_set != t) { 73 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t, NULL, NULL, NULL)); 74 user->time_bc_set = t; 75 } 76 PetscFunctionReturn(0); 77 } 78 79 // RHS (Explicit time-stepper) function setup 80 // This is the RHS of the ODE, given as u_t = G(t,u) 81 // This function takes in a state vector Q and writes into G 82 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data) { 83 User user = *(User *)user_data; 84 PetscScalar *q, *g; 85 Vec Q_loc = user->Q_loc, G_loc; 86 PetscMemType q_mem_type, g_mem_type; 87 PetscFunctionBeginUser; 88 89 // Get local vector 90 PetscCall(DMGetLocalVector(user->dm, &G_loc)); 91 92 // Update time dependent data 93 PetscCall(UpdateBoundaryValues(user, Q_loc, t)); 94 if (user->time != t) { 95 if (user->phys->solution_time_label) { 96 CeedOperatorContextSetDouble(user->op_rhs, user->phys->solution_time_label, &t); 97 } 98 user->time = t; 99 } 100 if (user->phys->timestep_size_label) { 101 PetscScalar dt; 102 PetscCall(TSGetTimeStep(ts, &dt)); 103 if (user->dt != dt) { 104 CeedOperatorContextSetDouble(user->op_rhs, user->phys->timestep_size_label, &dt); 105 user->dt = dt; 106 } 107 } 108 109 // Global-to-local 110 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 111 112 // Place PETSc vectors in CEED vectors 113 PetscCall(VecGetArrayReadAndMemType(Q_loc, (const PetscScalar **)&q, &q_mem_type)); 114 PetscCall(VecGetArrayAndMemType(G_loc, &g, &g_mem_type)); 115 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, q); 116 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 117 118 // Apply CEED operator 119 CeedOperatorApply(user->op_rhs, user->q_ceed, user->g_ceed, CEED_REQUEST_IMMEDIATE); 120 121 // Restore vectors 122 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 123 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 124 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q)); 125 PetscCall(VecRestoreArrayAndMemType(G_loc, &g)); 126 127 // Local-to-Global 128 PetscCall(VecZeroEntries(G)); 129 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G)); 130 131 // Inverse of the lumped mass matrix (M is Minv) 132 PetscCall(VecPointwiseMult(G, G, user->M)); 133 134 // Restore vectors 135 PetscCall(DMRestoreLocalVector(user->dm, &G_loc)); 136 137 PetscFunctionReturn(0); 138 } 139 140 // Implicit time-stepper function setup 141 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, void *user_data) { 142 User user = *(User *)user_data; 143 const PetscScalar *q, *q_dot; 144 PetscScalar *g; 145 Vec Q_loc = user->Q_loc, Q_dot_loc = user->Q_dot_loc, G_loc; 146 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type; 147 PetscFunctionBeginUser; 148 149 // Get local vectors 150 PetscCall(DMGetLocalVector(user->dm, &G_loc)); 151 152 // Update time dependent data 153 PetscCall(UpdateBoundaryValues(user, Q_loc, t)); 154 if (user->time != t) { 155 if (user->phys->solution_time_label) { 156 CeedOperatorContextSetDouble(user->op_ifunction, user->phys->solution_time_label, &t); 157 } 158 user->time = t; 159 } 160 if (user->phys->timestep_size_label) { 161 PetscScalar dt; 162 PetscCall(TSGetTimeStep(ts, &dt)); 163 if (user->dt != dt) { 164 CeedOperatorContextSetDouble(user->op_ifunction, user->phys->timestep_size_label, &dt); 165 user->dt = dt; 166 } 167 } 168 169 // Global-to-local 170 PetscCall(DMGlobalToLocalBegin(user->dm, Q, INSERT_VALUES, Q_loc)); 171 PetscCall(DMGlobalToLocalBegin(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc)); 172 PetscCall(DMGlobalToLocalEnd(user->dm, Q, INSERT_VALUES, Q_loc)); 173 PetscCall(DMGlobalToLocalEnd(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc)); 174 175 // Place PETSc vectors in CEED vectors 176 PetscCall(VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type)); 177 PetscCall(VecGetArrayReadAndMemType(Q_dot_loc, &q_dot, &q_dot_mem_type)); 178 PetscCall(VecGetArrayAndMemType(G_loc, &g, &g_mem_type)); 179 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, (PetscScalar *)q); 180 CeedVectorSetArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), CEED_USE_POINTER, (PetscScalar *)q_dot); 181 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 182 183 // Apply CEED operator 184 CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, CEED_REQUEST_IMMEDIATE); 185 186 // Restore vectors 187 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 188 CeedVectorTakeArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), NULL); 189 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 190 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, &q)); 191 PetscCall(VecRestoreArrayReadAndMemType(Q_dot_loc, &q_dot)); 192 PetscCall(VecRestoreArrayAndMemType(G_loc, &g)); 193 194 // Local-to-Global 195 PetscCall(VecZeroEntries(G)); 196 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G)); 197 198 // Restore vectors 199 PetscCall(DMRestoreLocalVector(user->dm, &G_loc)); 200 201 PetscFunctionReturn(0); 202 } 203 204 static PetscErrorCode MatMult_NS_IJacobian(Mat J, Vec Q, Vec G) { 205 User user; 206 const PetscScalar *q; 207 PetscScalar *g; 208 PetscMemType q_mem_type, g_mem_type; 209 PetscFunctionBeginUser; 210 PetscCall(MatShellGetContext(J, &user)); 211 Vec Q_loc = user->Q_dot_loc, // Note - Q_dot_loc has zero BCs 212 G_loc; 213 214 // Get local vectors 215 PetscCall(DMGetLocalVector(user->dm, &G_loc)); 216 217 // Global-to-local 218 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 219 220 // Place PETSc vectors in CEED vectors 221 PetscCall(VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type)); 222 PetscCall(VecGetArrayAndMemType(G_loc, &g, &g_mem_type)); 223 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, (PetscScalar *)q); 224 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 225 226 // Apply CEED operator 227 CeedOperatorApply(user->op_ijacobian, user->q_ceed, user->g_ceed, CEED_REQUEST_IMMEDIATE); 228 229 // Restore vectors 230 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 231 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 232 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, &q)); 233 PetscCall(VecRestoreArrayAndMemType(G_loc, &g)); 234 235 // Local-to-Global 236 PetscCall(VecZeroEntries(G)); 237 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G)); 238 239 // Restore vectors 240 PetscCall(DMRestoreLocalVector(user->dm, &G_loc)); 241 PetscFunctionReturn(0); 242 } 243 244 PetscErrorCode MatGetDiagonal_NS_IJacobian(Mat A, Vec D) { 245 User user; 246 Vec D_loc; 247 PetscScalar *d; 248 PetscMemType mem_type; 249 250 PetscFunctionBeginUser; 251 MatShellGetContext(A, &user); 252 PetscCall(DMGetLocalVector(user->dm, &D_loc)); 253 PetscCall(VecGetArrayAndMemType(D_loc, &d, &mem_type)); 254 CeedVectorSetArray(user->g_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, d); 255 CeedOperatorLinearAssembleDiagonal(user->op_ijacobian, user->g_ceed, CEED_REQUEST_IMMEDIATE); 256 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(mem_type), NULL); 257 PetscCall(VecRestoreArrayAndMemType(D_loc, &d)); 258 PetscCall(VecZeroEntries(D)); 259 PetscCall(DMLocalToGlobal(user->dm, D_loc, ADD_VALUES, D)); 260 PetscCall(DMRestoreLocalVector(user->dm, &D_loc)); 261 VecViewFromOptions(D, NULL, "-diag_vec_view"); 262 PetscFunctionReturn(0); 263 } 264 265 static PetscErrorCode FormPreallocation(User user, PetscBool pbdiagonal, Mat J, CeedVector *coo_values) { 266 PetscCount ncoo; 267 PetscInt *rows, *cols; 268 269 PetscFunctionBeginUser; 270 if (pbdiagonal) { 271 CeedSize l_size; 272 CeedOperatorGetActiveVectorLengths(user->op_ijacobian, &l_size, NULL); 273 ncoo = l_size * 5; 274 rows = malloc(ncoo * sizeof(rows[0])); 275 cols = malloc(ncoo * sizeof(cols[0])); 276 for (PetscCount n = 0; n < l_size / 5; n++) { 277 for (PetscInt i = 0; i < 5; i++) { 278 for (PetscInt j = 0; j < 5; j++) { 279 rows[(n * 5 + i) * 5 + j] = n * 5 + i; 280 cols[(n * 5 + i) * 5 + j] = n * 5 + j; 281 } 282 } 283 } 284 } else { 285 PetscCall(CeedOperatorLinearAssembleSymbolic(user->op_ijacobian, &ncoo, &rows, &cols)); 286 } 287 PetscCall(MatSetPreallocationCOOLocal(J, ncoo, rows, cols)); 288 free(rows); 289 free(cols); 290 CeedVectorCreate(user->ceed, ncoo, coo_values); 291 PetscFunctionReturn(0); 292 } 293 294 static PetscErrorCode FormSetValues(User user, PetscBool pbdiagonal, Mat J, CeedVector coo_values) { 295 CeedMemType mem_type = CEED_MEM_HOST; 296 const PetscScalar *values; 297 MatType mat_type; 298 299 PetscFunctionBeginUser; 300 PetscCall(MatGetType(J, &mat_type)); 301 if (strstr(mat_type, "kokkos") || strstr(mat_type, "cusparse")) mem_type = CEED_MEM_DEVICE; 302 if (user->app_ctx->pmat_pbdiagonal) { 303 CeedOperatorLinearAssemblePointBlockDiagonal(user->op_ijacobian, coo_values, CEED_REQUEST_IMMEDIATE); 304 } else { 305 CeedOperatorLinearAssemble(user->op_ijacobian, coo_values); 306 } 307 CeedVectorGetArrayRead(coo_values, mem_type, &values); 308 PetscCall(MatSetValuesCOO(J, values, INSERT_VALUES)); 309 CeedVectorRestoreArrayRead(coo_values, &values); 310 PetscFunctionReturn(0); 311 } 312 313 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) { 314 User user = *(User *)user_data; 315 PetscBool J_is_shell, J_is_mffd, J_pre_is_shell; 316 PetscFunctionBeginUser; 317 if (user->phys->ijacobian_time_shift_label) CeedOperatorContextSetDouble(user->op_ijacobian, user->phys->ijacobian_time_shift_label, &shift); 318 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd)); 319 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATSHELL, &J_is_shell)); 320 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATSHELL, &J_pre_is_shell)); 321 if (!user->matrices_set_up) { 322 if (J_is_shell) { 323 PetscCall(MatShellSetContext(J, user)); 324 PetscCall(MatShellSetOperation(J, MATOP_MULT, (void (*)(void))MatMult_NS_IJacobian)); 325 PetscCall(MatShellSetOperation(J, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_NS_IJacobian)); 326 PetscCall(MatSetUp(J)); 327 } 328 if (!J_pre_is_shell) { 329 PetscCall(FormPreallocation(user, user->app_ctx->pmat_pbdiagonal, J_pre, &user->coo_values_pmat)); 330 } 331 if (J != J_pre && !J_is_shell && !J_is_mffd) { 332 PetscCall(FormPreallocation(user, PETSC_FALSE, J, &user->coo_values_amat)); 333 } 334 user->matrices_set_up = true; 335 } 336 if (!J_pre_is_shell) { 337 PetscCall(FormSetValues(user, user->app_ctx->pmat_pbdiagonal, J_pre, user->coo_values_pmat)); 338 } 339 if (user->coo_values_amat) { 340 PetscCall(FormSetValues(user, PETSC_FALSE, J, user->coo_values_amat)); 341 } else if (J_is_mffd) { 342 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 343 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 344 } 345 PetscFunctionReturn(0); 346 } 347 348 PetscErrorCode WriteOutput(User user, Vec Q, PetscInt step_no, PetscScalar time) { 349 Vec Q_loc; 350 char file_path[PETSC_MAX_PATH_LEN]; 351 PetscViewer viewer; 352 PetscFunctionBeginUser; 353 354 if (user->app_ctx->checkpoint_vtk) { 355 // Set up output 356 PetscCall(DMGetLocalVector(user->dm, &Q_loc)); 357 PetscCall(PetscObjectSetName((PetscObject)Q_loc, "StateVec")); 358 PetscCall(VecZeroEntries(Q_loc)); 359 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 360 361 // Output 362 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no)); 363 364 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, FILE_MODE_WRITE, &viewer)); 365 PetscCall(VecView(Q_loc, viewer)); 366 PetscCall(PetscViewerDestroy(&viewer)); 367 if (user->dm_viz) { 368 Vec Q_refined, Q_refined_loc; 369 char file_path_refined[PETSC_MAX_PATH_LEN]; 370 PetscViewer viewer_refined; 371 372 PetscCall(DMGetGlobalVector(user->dm_viz, &Q_refined)); 373 PetscCall(DMGetLocalVector(user->dm_viz, &Q_refined_loc)); 374 PetscCall(PetscObjectSetName((PetscObject)Q_refined_loc, "Refined")); 375 376 PetscCall(MatInterpolate(user->interp_viz, Q, Q_refined)); 377 PetscCall(VecZeroEntries(Q_refined_loc)); 378 PetscCall(DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc)); 379 380 PetscCall( 381 PetscSNPrintf(file_path_refined, sizeof file_path_refined, "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no)); 382 383 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), file_path_refined, FILE_MODE_WRITE, &viewer_refined)); 384 PetscCall(VecView(Q_refined_loc, viewer_refined)); 385 PetscCall(DMRestoreLocalVector(user->dm_viz, &Q_refined_loc)); 386 PetscCall(DMRestoreGlobalVector(user->dm_viz, &Q_refined)); 387 PetscCall(PetscViewerDestroy(&viewer_refined)); 388 } 389 PetscCall(DMRestoreLocalVector(user->dm, &Q_loc)); 390 } 391 392 // Save data in a binary file for continuation of simulations 393 if (user->app_ctx->add_stepnum2bin) { 394 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution-%" PetscInt_FMT ".bin", user->app_ctx->output_dir, step_no)); 395 } else { 396 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", user->app_ctx->output_dir)); 397 } 398 PetscCall(PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer)); 399 400 PetscInt token = FLUIDS_FILE_TOKEN; 401 PetscCall(PetscViewerBinaryWrite(viewer, &token, 1, PETSC_INT)); 402 PetscCall(PetscViewerBinaryWrite(viewer, &step_no, 1, PETSC_INT)); 403 time /= user->units->second; // Dimensionalize time back 404 PetscCall(PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL)); 405 PetscCall(VecView(Q, viewer)); 406 PetscCall(PetscViewerDestroy(&viewer)); 407 PetscFunctionReturn(0); 408 } 409 410 // User provided TS Monitor 411 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) { 412 User user = ctx; 413 PetscFunctionBeginUser; 414 415 // Print every 'checkpoint_interval' steps 416 if (user->app_ctx->checkpoint_interval <= 0 || step_no % user->app_ctx->checkpoint_interval != 0 || 417 (user->app_ctx->cont_steps == step_no && step_no != 0)) 418 PetscFunctionReturn(0); 419 420 PetscCall(WriteOutput(user, Q, step_no, time)); 421 422 PetscFunctionReturn(0); 423 } 424 425 // TS: Create, setup, and solve 426 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, Vec *Q, PetscScalar *f_time, TS *ts) { 427 MPI_Comm comm = user->comm; 428 TSAdapt adapt; 429 PetscScalar final_time; 430 PetscFunctionBeginUser; 431 432 PetscCall(TSCreate(comm, ts)); 433 PetscCall(TSSetDM(*ts, dm)); 434 if (phys->implicit) { 435 PetscCall(TSSetType(*ts, TSBDF)); 436 if (user->op_ifunction) { 437 PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &user)); 438 } else { // Implicit integrators can fall back to using an RHSFunction 439 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user)); 440 } 441 if (user->op_ijacobian) { 442 PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &user)); 443 if (app_ctx->amat_type) { 444 Mat Pmat, Amat; 445 PetscCall(DMCreateMatrix(dm, &Pmat)); 446 PetscCall(DMSetMatType(dm, app_ctx->amat_type)); 447 PetscCall(DMCreateMatrix(dm, &Amat)); 448 PetscCall(TSSetIJacobian(*ts, Amat, Pmat, NULL, NULL)); 449 PetscCall(MatDestroy(&Amat)); 450 PetscCall(MatDestroy(&Pmat)); 451 } 452 } 453 } else { 454 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction"); 455 PetscCall(TSSetType(*ts, TSRK)); 456 PetscCall(TSRKSetType(*ts, TSRK5F)); 457 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user)); 458 } 459 PetscCall(TSSetMaxTime(*ts, 500. * user->units->second)); 460 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER)); 461 PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE)); 462 PetscCall(TSSetTimeStep(*ts, 1.e-2 * user->units->second)); 463 if (app_ctx->test_mode) { 464 PetscCall(TSSetMaxSteps(*ts, 10)); 465 } 466 PetscCall(TSGetAdapt(*ts, &adapt)); 467 PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * user->units->second, 1.e2 * user->units->second)); 468 PetscCall(TSSetFromOptions(*ts)); 469 user->time = user->time_bc_set = -1.0; // require all BCs and ctx to be updated 470 user->dt = -1.0; 471 if (!app_ctx->cont_steps) { // print initial condition 472 if (!app_ctx->test_mode) { 473 PetscCall(TSMonitor_NS(*ts, 0, 0., *Q, user)); 474 } 475 } else { // continue from time of last output 476 PetscInt count; 477 PetscViewer viewer; 478 479 if (app_ctx->cont_time <= 0) { // Legacy files did not include step number and time 480 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_time_file, FILE_MODE_READ, &viewer)); 481 PetscCall(PetscViewerBinaryRead(viewer, &app_ctx->cont_time, 1, &count, PETSC_REAL)); 482 PetscCall(PetscViewerDestroy(&viewer)); 483 PetscCheck(app_ctx->cont_steps != -1, comm, PETSC_ERR_ARG_INCOMP, 484 "-continue step number not specified, but checkpoint file does not contain a step number (likely written by older code version)"); 485 } 486 PetscCall(TSSetTime(*ts, app_ctx->cont_time * user->units->second)); 487 PetscCall(TSSetStepNumber(*ts, app_ctx->cont_steps)); 488 } 489 if (!app_ctx->test_mode) { 490 PetscCall(TSMonitorSet(*ts, TSMonitor_NS, user, NULL)); 491 if (app_ctx->turb_spanstats_enable) { 492 PetscCall(TSMonitorSet(*ts, TSMonitor_Statistics, user, NULL)); 493 user->spanstats.prev_time = app_ctx->cont_time * user->units->second; 494 } 495 } 496 497 // Solve 498 PetscReal start_time; 499 PetscInt start_step; 500 PetscCall(TSGetTime(*ts, &start_time)); 501 PetscCall(TSGetStepNumber(*ts, &start_step)); 502 503 PetscPreLoadBegin(PETSC_FALSE, "Fluids Solve"); 504 PetscCall(TSSetTime(*ts, start_time)); 505 PetscCall(TSSetStepNumber(*ts, start_step)); 506 if (PetscPreLoadingOn) { 507 // LCOV_EXCL_START 508 SNES snes; 509 Vec Q_preload; 510 PetscReal rtol; 511 PetscCall(VecDuplicate(*Q, &Q_preload)); 512 PetscCall(VecCopy(*Q, Q_preload)); 513 PetscCall(TSGetSNES(*ts, &snes)); 514 PetscCall(SNESGetTolerances(snes, NULL, &rtol, NULL, NULL, NULL)); 515 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, .99, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 516 PetscCall(TSSetSolution(*ts, *Q)); 517 PetscCall(TSStep(*ts)); 518 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 519 PetscCall(VecDestroy(&Q_preload)); 520 // LCOV_EXCL_STOP 521 } else { 522 PetscCall(PetscBarrier((PetscObject)*ts)); 523 PetscCall(TSSolve(*ts, *Q)); 524 } 525 PetscPreLoadEnd(); 526 527 PetscCall(TSGetSolveTime(*ts, &final_time)); 528 *f_time = final_time; 529 530 if (!app_ctx->test_mode) { 531 PetscInt step_no; 532 PetscCall(TSGetStepNumber(*ts, &step_no)); 533 if (user->app_ctx->checkpoint_interval > 0 || user->app_ctx->checkpoint_interval == -1) { 534 PetscCall(WriteOutput(user, *Q, step_no, final_time)); 535 } 536 537 PetscLogEvent stage_id; 538 PetscStageLog stage_log; 539 540 PetscCall(PetscLogStageGetId("Fluids Solve", &stage_id)); 541 PetscCall(PetscLogGetStageLog(&stage_log)); 542 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time taken for solution (sec): %g\n", stage_log->stageInfo[stage_id].perfInfo.time)); 543 } 544 PetscFunctionReturn(0); 545 } 546