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