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") 329 || strstr(mat_type, "cusparse")) mem_type = CEED_MEM_DEVICE; 330 CeedOperatorLinearAssemble(user->op_ijacobian, user->coo_values); 331 CeedVectorGetArrayRead(user->coo_values, mem_type, &values); 332 if (coo_vec) { 333 VecPlaceArray(coo_vec, values); 334 VecViewFromOptions(coo_vec, NULL, "-coo_vec_view"); 335 VecDestroy(&coo_vec); 336 } 337 PetscCall(MatSetValuesCOO(J_pre, values, INSERT_VALUES)); 338 CeedVectorRestoreArrayRead(user->coo_values, &values); 339 } 340 PetscFunctionReturn(0); 341 } 342 343 // User provided TS Monitor 344 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, 345 Vec Q, void *ctx) { 346 User user = ctx; 347 Vec Q_loc; 348 char file_path[PETSC_MAX_PATH_LEN]; 349 PetscViewer viewer; 350 PetscErrorCode ierr; 351 PetscFunctionBeginUser; 352 353 // Print every 'output_freq' steps 354 if (step_no % user->app_ctx->output_freq != 0) 355 PetscFunctionReturn(0); 356 357 // Set up output 358 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 359 ierr = PetscObjectSetName((PetscObject)Q_loc, "StateVec"); CHKERRQ(ierr); 360 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 361 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 362 363 // Output 364 ierr = PetscSNPrintf(file_path, sizeof file_path, 365 "%s/ns-%03" PetscInt_FMT ".vtu", 366 user->app_ctx->output_dir, step_no + user->app_ctx->cont_steps); 367 CHKERRQ(ierr); 368 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, 369 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 370 ierr = VecView(Q_loc, viewer); CHKERRQ(ierr); 371 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 372 if (user->dm_viz) { 373 Vec Q_refined, Q_refined_loc; 374 char file_path_refined[PETSC_MAX_PATH_LEN]; 375 PetscViewer viewer_refined; 376 377 ierr = DMGetGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 378 ierr = DMGetLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 379 ierr = PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"); 380 CHKERRQ(ierr); 381 ierr = MatInterpolate(user->interp_viz, Q, Q_refined); CHKERRQ(ierr); 382 ierr = VecZeroEntries(Q_refined_loc); CHKERRQ(ierr); 383 ierr = DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc); 384 CHKERRQ(ierr); 385 ierr = PetscSNPrintf(file_path_refined, sizeof file_path_refined, 386 "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, 387 step_no + user->app_ctx->cont_steps); 388 CHKERRQ(ierr); 389 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), 390 file_path_refined, FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 391 ierr = VecView(Q_refined_loc, viewer_refined); CHKERRQ(ierr); 392 ierr = DMRestoreLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 393 ierr = DMRestoreGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 394 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 395 } 396 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 397 398 // Save data in a binary file for continuation of simulations 399 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", 400 user->app_ctx->output_dir); CHKERRQ(ierr); 401 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 402 CHKERRQ(ierr); 403 ierr = VecView(Q, viewer); CHKERRQ(ierr); 404 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 405 406 // Save time stamp 407 // Dimensionalize time back 408 time /= user->units->second; 409 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 410 user->app_ctx->output_dir); CHKERRQ(ierr); 411 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 412 CHKERRQ(ierr); 413 #if PETSC_VERSION_GE(3,13,0) 414 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 415 #else 416 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 417 #endif 418 CHKERRQ(ierr); 419 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 420 421 PetscFunctionReturn(0); 422 } 423 424 // TS: Create, setup, and solve 425 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, 426 Vec *Q, PetscScalar *f_time, TS *ts) { 427 MPI_Comm comm = user->comm; 428 TSAdapt adapt; 429 PetscScalar final_time; 430 PetscErrorCode ierr; 431 PetscFunctionBeginUser; 432 433 ierr = TSCreate(comm, ts); CHKERRQ(ierr); 434 ierr = TSSetDM(*ts, dm); CHKERRQ(ierr); 435 if (phys->implicit) { 436 ierr = TSSetType(*ts, TSBDF); CHKERRQ(ierr); 437 if (user->op_ifunction) { 438 ierr = TSSetIFunction(*ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 439 } else { // Implicit integrators can fall back to using an RHSFunction 440 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 441 } 442 if (user->op_ijacobian) { 443 ierr = DMTSSetIJacobian(dm, FormIJacobian_NS, &user); CHKERRQ(ierr); 444 } 445 } else { 446 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 447 "Problem does not provide RHSFunction"); 448 ierr = TSSetType(*ts, TSRK); CHKERRQ(ierr); 449 ierr = TSRKSetType(*ts, TSRK5F); CHKERRQ(ierr); 450 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 451 } 452 ierr = TSSetMaxTime(*ts, 500. * user->units->second); CHKERRQ(ierr); 453 ierr = TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 454 ierr = TSSetTimeStep(*ts, 1.e-2 * user->units->second); CHKERRQ(ierr); 455 if (app_ctx->test_mode) {ierr = TSSetMaxSteps(*ts, 10); CHKERRQ(ierr);} 456 ierr = TSGetAdapt(*ts, &adapt); CHKERRQ(ierr); 457 ierr = TSAdaptSetStepLimits(adapt, 458 1.e-12 * user->units->second, 459 1.e2 * user->units->second); CHKERRQ(ierr); 460 ierr = TSSetFromOptions(*ts); CHKERRQ(ierr); 461 user->time = -1.0; // require all BCs and ctx to be updated 462 user->dt = -1.0; 463 if (!app_ctx->cont_steps) { // print initial condition 464 if (!app_ctx->test_mode) { 465 ierr = TSMonitor_NS(*ts, 0, 0., *Q, user); CHKERRQ(ierr); 466 } 467 } else { // continue from time of last output 468 PetscReal time; 469 PetscInt count; 470 PetscViewer viewer; 471 char file_path[PETSC_MAX_PATH_LEN]; 472 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 473 app_ctx->output_dir); CHKERRQ(ierr); 474 ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer); 475 CHKERRQ(ierr); 476 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 477 CHKERRQ(ierr); 478 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 479 ierr = TSSetTime(*ts, time * user->units->second); CHKERRQ(ierr); 480 } 481 if (!app_ctx->test_mode) { 482 ierr = TSMonitorSet(*ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 483 } 484 485 // Solve 486 double start, cpu_time_used; 487 start = MPI_Wtime(); 488 ierr = PetscBarrier((PetscObject) *ts); CHKERRQ(ierr); 489 ierr = TSSolve(*ts, *Q); CHKERRQ(ierr); 490 cpu_time_used = MPI_Wtime() - start; 491 ierr = TSGetSolveTime(*ts, &final_time); CHKERRQ(ierr); 492 *f_time = final_time; 493 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 494 comm); CHKERRQ(ierr); 495 if (!app_ctx->test_mode) { 496 ierr = PetscPrintf(PETSC_COMM_WORLD, 497 "Time taken for solution (sec): %g\n", 498 (double)cpu_time_used); CHKERRQ(ierr); 499 } 500 PetscFunctionReturn(0); 501 } 502