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