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