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