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 = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t, 97 NULL, NULL, NULL); CHKERRQ(ierr); 98 if (user->phys->solution_time_label) { 99 CeedOperatorContextSetDouble(user->op_rhs, user->phys->solution_time_label, &t); 100 } 101 user->time = t; 102 } 103 if (user->phys->timestep_size_label) { 104 PetscScalar dt; 105 ierr = TSGetTimeStep(ts, &dt); CHKERRQ(ierr); 106 if (user->dt != dt) { 107 CeedOperatorContextSetDouble(user->op_rhs, user->phys->timestep_size_label, 108 &dt); 109 user->dt = dt; 110 } 111 } 112 113 // Global-to-local 114 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 115 116 // Place PETSc vectors in CEED vectors 117 ierr = VecGetArrayReadAndMemType(Q_loc, (const PetscScalar **)&q, &q_mem_type); 118 CHKERRQ(ierr); 119 ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr); 120 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, q); 121 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 122 123 // Apply CEED operator 124 CeedOperatorApply(user->op_rhs, user->q_ceed, user->g_ceed, 125 CEED_REQUEST_IMMEDIATE); 126 127 // Restore vectors 128 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 129 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 130 ierr = VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q); 131 CHKERRQ(ierr); 132 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 133 134 // Local-to-Global 135 ierr = VecZeroEntries(G); CHKERRQ(ierr); 136 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 137 138 // Inverse of the lumped mass matrix (M is Minv) 139 ierr = VecPointwiseMult(G, G, user->M); CHKERRQ(ierr); 140 141 // Restore vectors 142 ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 143 144 PetscFunctionReturn(0); 145 } 146 147 // Implicit time-stepper function setup 148 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, 149 void *user_data) { 150 User user = *(User *)user_data; 151 const PetscScalar *q, *q_dot; 152 PetscScalar *g; 153 Vec Q_loc = user->Q_loc, Q_dot_loc = user->Q_dot_loc, G_loc; 154 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type; 155 PetscErrorCode ierr; 156 PetscFunctionBeginUser; 157 158 // Get local vectors 159 ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 160 161 // Update time dependent data 162 if (user->time != t) { 163 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t, 164 NULL, NULL, NULL); CHKERRQ(ierr); 165 if (user->phys->solution_time_label) { 166 CeedOperatorContextSetDouble(user->op_ifunction, 167 user->phys->solution_time_label, &t); 168 } 169 user->time = t; 170 } 171 if (user->phys->timestep_size_label) { 172 PetscScalar dt; 173 ierr = TSGetTimeStep(ts, &dt); CHKERRQ(ierr); 174 if (user->dt != dt) { 175 CeedOperatorContextSetDouble(user->op_ifunction, 176 user->phys->timestep_size_label, &dt); 177 user->dt = dt; 178 } 179 } 180 181 // Global-to-local 182 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 183 ierr = DMGlobalToLocal(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc); 184 CHKERRQ(ierr); 185 186 // Place PETSc vectors in CEED vectors 187 ierr = VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type); CHKERRQ(ierr); 188 ierr = VecGetArrayReadAndMemType(Q_dot_loc, &q_dot, &q_dot_mem_type); 189 CHKERRQ(ierr); 190 ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr); 191 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, 192 (PetscScalar *)q); 193 CeedVectorSetArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), 194 CEED_USE_POINTER, (PetscScalar *)q_dot); 195 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 196 197 // Apply CEED operator 198 CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, 199 CEED_REQUEST_IMMEDIATE); 200 201 // Restore vectors 202 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 203 CeedVectorTakeArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), NULL); 204 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 205 ierr = VecRestoreArrayReadAndMemType(Q_loc, &q); CHKERRQ(ierr); 206 ierr = VecRestoreArrayReadAndMemType(Q_dot_loc, &q_dot); CHKERRQ(ierr); 207 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 208 209 // Local-to-Global 210 ierr = VecZeroEntries(G); CHKERRQ(ierr); 211 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 212 213 // Restore vectors 214 ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 215 216 PetscFunctionReturn(0); 217 } 218 219 static PetscErrorCode MatMult_NS_IJacobian(Mat J, Vec Q, Vec G) { 220 User user; 221 const PetscScalar *q; 222 PetscScalar *g; 223 PetscMemType q_mem_type, g_mem_type; 224 PetscErrorCode ierr; 225 PetscFunctionBeginUser; 226 ierr = MatShellGetContext(J, &user); CHKERRQ(ierr); 227 Vec Q_loc = user->Q_dot_loc, // Note - Q_dot_loc has zero BCs 228 G_loc; 229 230 // Get local vectors 231 ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 232 233 // Global-to-local 234 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 235 236 // Place PETSc vectors in CEED vectors 237 ierr = VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type); CHKERRQ(ierr); 238 ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr); 239 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, 240 (PetscScalar *)q); 241 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 242 243 // Apply CEED operator 244 CeedOperatorApply(user->op_ijacobian, user->q_ceed, user->g_ceed, 245 CEED_REQUEST_IMMEDIATE); 246 247 // Restore vectors 248 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 249 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 250 ierr = VecRestoreArrayReadAndMemType(Q_loc, &q); CHKERRQ(ierr); 251 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 252 253 // Local-to-Global 254 ierr = VecZeroEntries(G); CHKERRQ(ierr); 255 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 256 257 // Restore vectors 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 static PetscErrorCode FormPreallocation(User user, PetscBool pbdiagonal, Mat J, 285 CeedVector *coo_values) { 286 PetscCount ncoo; 287 PetscInt *rows, *cols; 288 289 PetscFunctionBeginUser; 290 if (pbdiagonal) { 291 CeedSize l_size; 292 CeedOperatorGetActiveVectorLengths(user->op_ijacobian, &l_size, NULL); 293 ncoo = l_size * 5; 294 rows = malloc(ncoo*sizeof(rows[0])); 295 cols = malloc(ncoo*sizeof(cols[0])); 296 for (PetscCount n=0; n<l_size/5; n++) { 297 for (PetscInt i=0; i<5; i++) { 298 for (PetscInt j=0; j<5; j++) { 299 rows[(n*5+i)*5+j] = n * 5 + i; 300 cols[(n*5+i)*5+j] = n * 5 + j; 301 } 302 } 303 } 304 } else { 305 PetscCall(CeedOperatorLinearAssembleSymbolic(user->op_ijacobian, &ncoo, &rows, 306 &cols)); 307 } 308 PetscCall(MatSetPreallocationCOOLocal(J, ncoo, rows, cols)); 309 free(rows); 310 free(cols); 311 CeedVectorCreate(user->ceed, ncoo, coo_values); 312 PetscFunctionReturn(0); 313 } 314 315 static PetscErrorCode FormSetValues(User user, PetscBool pbdiagonal, Mat J, 316 CeedVector coo_values) { 317 CeedMemType mem_type = CEED_MEM_HOST; 318 const PetscScalar *values; 319 MatType mat_type; 320 321 PetscFunctionBeginUser; 322 PetscCall(MatGetType(J, &mat_type)); 323 if (strstr(mat_type, "kokkos") || strstr(mat_type, "cusparse")) 324 mem_type = CEED_MEM_DEVICE; 325 if (user->app_ctx->pmat_pbdiagonal) { 326 CeedOperatorLinearAssemblePointBlockDiagonal(user->op_ijacobian, 327 coo_values, CEED_REQUEST_IMMEDIATE); 328 } else { 329 CeedOperatorLinearAssemble(user->op_ijacobian, coo_values); 330 } 331 CeedVectorGetArrayRead(coo_values, mem_type, &values); 332 PetscCall(MatSetValuesCOO(J, values, INSERT_VALUES)); 333 CeedVectorRestoreArrayRead(coo_values, &values); 334 PetscFunctionReturn(0); 335 } 336 337 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, 338 PetscReal shift, Mat J, Mat J_pre, 339 void *user_data) { 340 User user = *(User *)user_data; 341 PetscBool J_is_shell, J_is_mffd, J_pre_is_shell; 342 PetscFunctionBeginUser; 343 if (user->phys->ijacobian_time_shift_label) 344 CeedOperatorContextSetDouble(user->op_ijacobian, 345 user->phys->ijacobian_time_shift_label, &shift); 346 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd)); 347 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATSHELL, &J_is_shell)); 348 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATSHELL, 349 &J_pre_is_shell)); 350 if (!user->matrices_set_up) { 351 if (J_is_shell) { 352 PetscCall(MatShellSetContext(J, user)); 353 PetscCall(MatShellSetOperation(J, MATOP_MULT, 354 (void (*)(void))MatMult_NS_IJacobian)); 355 PetscCall(MatShellSetOperation(J, MATOP_GET_DIAGONAL, 356 (void (*)(void))MatGetDiagonal_NS_IJacobian)); 357 PetscCall(MatSetUp(J)); 358 } 359 if (!J_pre_is_shell) { 360 PetscCall(FormPreallocation(user,user->app_ctx->pmat_pbdiagonal,J_pre, 361 &user->coo_values_pmat)); 362 } 363 if (J != J_pre && !J_is_shell && !J_is_mffd) { 364 PetscCall(FormPreallocation(user,PETSC_FALSE,J, &user->coo_values_amat)); 365 } 366 user->matrices_set_up = true; 367 } 368 if (!J_pre_is_shell) { 369 PetscCall(FormSetValues(user, user->app_ctx->pmat_pbdiagonal, J_pre, 370 user->coo_values_pmat)); 371 } 372 if (user->coo_values_amat) { 373 PetscCall(FormSetValues(user, PETSC_FALSE, J, user->coo_values_amat)); 374 } else if (J_is_mffd) { 375 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 376 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 377 } 378 PetscFunctionReturn(0); 379 } 380 381 // User provided TS Monitor 382 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, 383 Vec Q, void *ctx) { 384 User user = ctx; 385 Vec Q_loc; 386 char file_path[PETSC_MAX_PATH_LEN]; 387 PetscViewer viewer; 388 PetscErrorCode ierr; 389 PetscFunctionBeginUser; 390 391 // Print every 'output_freq' steps 392 if (step_no % user->app_ctx->output_freq != 0) 393 PetscFunctionReturn(0); 394 395 // Set up output 396 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 397 ierr = PetscObjectSetName((PetscObject)Q_loc, "StateVec"); CHKERRQ(ierr); 398 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 399 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 400 401 // Output 402 ierr = PetscSNPrintf(file_path, sizeof file_path, 403 "%s/ns-%03" PetscInt_FMT ".vtu", 404 user->app_ctx->output_dir, step_no + user->app_ctx->cont_steps); 405 CHKERRQ(ierr); 406 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, 407 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 408 ierr = VecView(Q_loc, viewer); CHKERRQ(ierr); 409 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 410 if (user->dm_viz) { 411 Vec Q_refined, Q_refined_loc; 412 char file_path_refined[PETSC_MAX_PATH_LEN]; 413 PetscViewer viewer_refined; 414 415 ierr = DMGetGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 416 ierr = DMGetLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 417 ierr = PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"); 418 CHKERRQ(ierr); 419 ierr = MatInterpolate(user->interp_viz, Q, Q_refined); CHKERRQ(ierr); 420 ierr = VecZeroEntries(Q_refined_loc); CHKERRQ(ierr); 421 ierr = DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc); 422 CHKERRQ(ierr); 423 ierr = PetscSNPrintf(file_path_refined, sizeof file_path_refined, 424 "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, 425 step_no + user->app_ctx->cont_steps); 426 CHKERRQ(ierr); 427 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), 428 file_path_refined, FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 429 ierr = VecView(Q_refined_loc, viewer_refined); CHKERRQ(ierr); 430 ierr = DMRestoreLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 431 ierr = DMRestoreGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 432 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 433 } 434 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 435 436 // Save data in a binary file for continuation of simulations 437 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", 438 user->app_ctx->output_dir); CHKERRQ(ierr); 439 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 440 CHKERRQ(ierr); 441 ierr = VecView(Q, viewer); CHKERRQ(ierr); 442 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 443 444 // Save time stamp 445 // Dimensionalize time back 446 time /= user->units->second; 447 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 448 user->app_ctx->output_dir); CHKERRQ(ierr); 449 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 450 CHKERRQ(ierr); 451 #if PETSC_VERSION_GE(3,13,0) 452 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 453 #else 454 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 455 #endif 456 CHKERRQ(ierr); 457 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 458 459 PetscFunctionReturn(0); 460 } 461 462 // TS: Create, setup, and solve 463 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, 464 Vec *Q, PetscScalar *f_time, TS *ts) { 465 MPI_Comm comm = user->comm; 466 TSAdapt adapt; 467 PetscScalar final_time; 468 PetscErrorCode ierr; 469 PetscFunctionBeginUser; 470 471 ierr = TSCreate(comm, ts); CHKERRQ(ierr); 472 ierr = TSSetDM(*ts, dm); CHKERRQ(ierr); 473 if (phys->implicit) { 474 ierr = TSSetType(*ts, TSBDF); CHKERRQ(ierr); 475 if (user->op_ifunction) { 476 ierr = TSSetIFunction(*ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 477 } else { // Implicit integrators can fall back to using an RHSFunction 478 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 479 } 480 if (user->op_ijacobian) { 481 ierr = DMTSSetIJacobian(dm, FormIJacobian_NS, &user); CHKERRQ(ierr); 482 if (app_ctx->amat_type) { 483 Mat Pmat,Amat; 484 ierr = DMCreateMatrix(dm, &Pmat); CHKERRQ(ierr); 485 ierr = DMSetMatType(dm, app_ctx->amat_type); CHKERRQ(ierr); 486 ierr = DMCreateMatrix(dm, &Amat); CHKERRQ(ierr); 487 ierr = TSSetIJacobian(*ts, Amat, Pmat, NULL, NULL); CHKERRQ(ierr); 488 ierr = MatDestroy(&Amat); CHKERRQ(ierr); 489 ierr = MatDestroy(&Pmat); CHKERRQ(ierr); 490 } 491 } 492 } else { 493 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 494 "Problem does not provide RHSFunction"); 495 ierr = TSSetType(*ts, TSRK); CHKERRQ(ierr); 496 ierr = TSRKSetType(*ts, TSRK5F); CHKERRQ(ierr); 497 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 498 } 499 ierr = TSSetMaxTime(*ts, 500. * user->units->second); CHKERRQ(ierr); 500 ierr = TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 501 ierr = TSSetTimeStep(*ts, 1.e-2 * user->units->second); CHKERRQ(ierr); 502 if (app_ctx->test_mode) {ierr = TSSetMaxSteps(*ts, 10); CHKERRQ(ierr);} 503 ierr = TSGetAdapt(*ts, &adapt); CHKERRQ(ierr); 504 ierr = TSAdaptSetStepLimits(adapt, 505 1.e-12 * user->units->second, 506 1.e2 * user->units->second); CHKERRQ(ierr); 507 ierr = TSSetFromOptions(*ts); CHKERRQ(ierr); 508 user->time = -1.0; // require all BCs and ctx to be updated 509 user->dt = -1.0; 510 if (!app_ctx->cont_steps) { // print initial condition 511 if (!app_ctx->test_mode) { 512 ierr = TSMonitor_NS(*ts, 0, 0., *Q, user); CHKERRQ(ierr); 513 } 514 } else { // continue from time of last output 515 PetscReal time; 516 PetscInt count; 517 PetscViewer viewer; 518 char file_path[PETSC_MAX_PATH_LEN]; 519 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 520 app_ctx->output_dir); CHKERRQ(ierr); 521 ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer); 522 CHKERRQ(ierr); 523 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 524 CHKERRQ(ierr); 525 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 526 ierr = TSSetTime(*ts, time * user->units->second); CHKERRQ(ierr); 527 } 528 if (!app_ctx->test_mode) { 529 ierr = TSMonitorSet(*ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 530 } 531 532 // Solve 533 PetscScalar start_time; 534 ierr = TSGetTime(*ts, &start_time); CHKERRQ(ierr); 535 536 PetscPreLoadBegin(PETSC_FALSE, "Fluids Solve"); 537 PetscCall(TSSetTime(*ts, start_time)); 538 PetscCall(TSSetStepNumber(*ts, 0)); 539 if (PetscPreLoadingOn) { 540 // LCOV_EXCL_START 541 SNES snes; 542 Vec Q_preload; 543 PetscReal rtol; 544 PetscCall(VecDuplicate(*Q, &Q_preload)); 545 PetscCall(VecCopy(*Q, Q_preload)); 546 PetscCall(TSGetSNES(*ts, &snes)); 547 PetscCall(SNESGetTolerances(snes, NULL, &rtol, NULL, NULL, NULL)); 548 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, .99, PETSC_DEFAULT, 549 PETSC_DEFAULT, PETSC_DEFAULT)); 550 PetscCall(TSSetSolution(*ts, *Q)); 551 PetscCall(TSStep(*ts)); 552 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, rtol, PETSC_DEFAULT, 553 PETSC_DEFAULT, PETSC_DEFAULT)); 554 PetscCall(VecDestroy(&Q_preload)); 555 // LCOV_EXCL_STOP 556 } else { 557 ierr = PetscBarrier((PetscObject) *ts); CHKERRQ(ierr); 558 ierr = TSSolve(*ts, *Q); CHKERRQ(ierr); 559 } 560 PetscPreLoadEnd(); 561 562 PetscCall(TSGetSolveTime(*ts, &final_time)); 563 *f_time = final_time; 564 565 if (!app_ctx->test_mode) { 566 PetscLogEvent stage_id; 567 PetscStageLog stage_log; 568 569 PetscCall(PetscLogStageGetId("Fluids Solve", &stage_id)); 570 PetscCall(PetscLogGetStageLog(&stage_log)); 571 PetscCall(PetscPrintf(PETSC_COMM_WORLD, 572 "Time taken for solution (sec): %g\n", 573 stage_log->stageInfo[stage_id].perfInfo.time)); 574 } 575 PetscFunctionReturn(0); 576 } 577