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/newtonian_state.h" 13 14 // Compute mass matrix for explicit scheme 15 PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, CeedData ceed_data, Vec M) { 16 Vec M_loc; 17 CeedQFunction qf_mass; 18 CeedOperator op_mass; 19 CeedVector m_ceed, ones_vec; 20 CeedInt num_comp_q, q_data_size; 21 PetscFunctionBeginUser; 22 23 // CEED Restriction 24 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q); 25 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); 26 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &m_ceed, NULL); 27 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &ones_vec, NULL); 28 CeedVectorSetValue(ones_vec, 1.0); 29 30 // CEED QFunction 31 PetscCall(CreateMassQFunction(ceed, num_comp_q, q_data_size, &qf_mass)); 32 33 // CEED Operator 34 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 35 CeedOperatorSetField(op_mass, "u", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 36 CeedOperatorSetField(op_mass, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data); 37 CeedOperatorSetField(op_mass, "v", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE); 38 39 // Place PETSc vector in CEED vector 40 CeedScalar *m; 41 PetscMemType m_mem_type; 42 PetscCall(DMGetLocalVector(dm, &M_loc)); 43 PetscCall(VecGetArrayAndMemType(M_loc, (PetscScalar **)&m, &m_mem_type)); 44 CeedVectorSetArray(m_ceed, MemTypeP2C(m_mem_type), CEED_USE_POINTER, m); 45 46 // Apply CEED Operator 47 CeedOperatorApply(op_mass, ones_vec, m_ceed, CEED_REQUEST_IMMEDIATE); 48 49 // Restore vectors 50 CeedVectorTakeArray(m_ceed, MemTypeP2C(m_mem_type), NULL); 51 PetscCall(VecRestoreArrayReadAndMemType(M_loc, (const PetscScalar **)&m)); 52 53 // Local-to-Global 54 PetscCall(VecZeroEntries(M)); 55 PetscCall(DMLocalToGlobal(dm, M_loc, ADD_VALUES, M)); 56 PetscCall(DMRestoreLocalVector(dm, &M_loc)); 57 58 // Invert diagonally lumped mass vector for RHS function 59 PetscCall(VecReciprocal(M)); 60 61 // Cleanup 62 CeedVectorDestroy(&ones_vec); 63 CeedVectorDestroy(&m_ceed); 64 CeedQFunctionDestroy(&qf_mass); 65 CeedOperatorDestroy(&op_mass); 66 67 PetscFunctionReturn(0); 68 } 69 70 // Insert Boundary values if it's a new time 71 PetscErrorCode UpdateBoundaryValues(User user, Vec Q_loc, PetscReal t) { 72 PetscFunctionBeginUser; 73 if (user->time_bc_set != t) { 74 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, t, NULL, NULL, NULL)); 75 user->time_bc_set = t; 76 } 77 PetscFunctionReturn(0); 78 } 79 80 PetscErrorCode UpdateContextLabel(PetscScalar *previous_value, PetscScalar update_value, CeedOperator op, CeedContextFieldLabel label) { 81 PetscFunctionBeginUser; 82 83 if (*previous_value != update_value) { 84 if (label) { 85 CeedOperatorContextSetDouble(op, label, &update_value); 86 } 87 *previous_value = update_value; 88 } 89 PetscFunctionReturn(0); 90 } 91 92 // RHS (Explicit time-stepper) function setup 93 // This is the RHS of the ODE, given as u_t = G(t,u) 94 // This function takes in a state vector Q and writes into G 95 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data) { 96 User user = *(User *)user_data; 97 PetscScalar *q, *g, dt; 98 Vec Q_loc = user->Q_loc, G_loc; 99 PetscMemType q_mem_type, g_mem_type; 100 PetscFunctionBeginUser; 101 102 // Get local vector 103 PetscCall(DMGetLocalVector(user->dm, &G_loc)); 104 105 // Update time dependent data 106 PetscCall(UpdateBoundaryValues(user, Q_loc, t)); 107 PetscCall(UpdateContextLabel(&user->time, t, user->op_rhs, user->phys->solution_time_label)); 108 PetscCall(TSGetTimeStep(ts, &dt)); 109 PetscCall(UpdateContextLabel(&user->dt, dt, user->op_rhs, user->phys->timestep_size_label)); 110 111 // Global-to-local 112 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 113 114 // Place PETSc vectors in CEED vectors 115 PetscCall(VecGetArrayReadAndMemType(Q_loc, (const PetscScalar **)&q, &q_mem_type)); 116 PetscCall(VecGetArrayAndMemType(G_loc, &g, &g_mem_type)); 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, 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 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q)); 127 PetscCall(VecRestoreArrayAndMemType(G_loc, &g)); 128 129 // Local-to-Global 130 PetscCall(VecZeroEntries(G)); 131 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G)); 132 133 // Inverse of the lumped mass matrix (M is Minv) 134 PetscCall(VecPointwiseMult(G, G, user->M)); 135 136 // Restore vectors 137 PetscCall(DMRestoreLocalVector(user->dm, &G_loc)); 138 139 PetscFunctionReturn(0); 140 } 141 142 // Surface forces function setup 143 static PetscErrorCode Surface_Forces_NS(DM dm, Vec G_loc, PetscInt num_walls, const PetscInt walls[], PetscScalar *reaction_force) { 144 DMLabel face_label; 145 const PetscScalar *g; 146 PetscInt dim = 3; 147 MPI_Comm comm; 148 149 PetscFunctionBeginUser; 150 PetscCall(PetscArrayzero(reaction_force, num_walls * dim)); 151 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 152 PetscCall(DMGetLabel(dm, "Face Sets", &face_label)); 153 PetscCall(VecGetArrayRead(G_loc, &g)); 154 for (PetscInt w = 0; w < num_walls; w++) { 155 const PetscInt wall = walls[w]; 156 IS wall_is; 157 PetscCall(DMLabelGetStratumIS(face_label, wall, &wall_is)); 158 if (wall_is) { // There exist such points on this process 159 PetscInt num_points; 160 const PetscInt *points; 161 PetscCall(ISGetSize(wall_is, &num_points)); 162 PetscCall(ISGetIndices(wall_is, &points)); 163 for (PetscInt i = 0; i < num_points; i++) { 164 const PetscInt p = points[i]; 165 const StateConservative *r; 166 PetscCall(DMPlexPointLocalRead(dm, p, g, &r)); 167 if (!r) continue; 168 for (PetscInt j = 0; j < 3; j++) { 169 reaction_force[w * dim + j] -= r->momentum[j]; 170 } 171 } 172 PetscCall(ISRestoreIndices(wall_is, &points)); 173 } 174 PetscCall(ISDestroy(&wall_is)); 175 } 176 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, reaction_force, dim * num_walls, MPIU_SCALAR, MPI_SUM, comm)); 177 // Restore Vectors 178 PetscCall(VecRestoreArrayRead(G_loc, &g)); 179 180 PetscFunctionReturn(0); 181 } 182 183 // Implicit time-stepper function setup 184 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, void *user_data) { 185 User user = *(User *)user_data; 186 const PetscScalar *q, *q_dot; 187 PetscScalar *g, dt; 188 Vec Q_loc = user->Q_loc, Q_dot_loc = user->Q_dot_loc, G_loc; 189 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type; 190 PetscFunctionBeginUser; 191 192 // Get local vectors 193 PetscCall(DMGetNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 194 195 // Update time dependent data 196 PetscCall(UpdateBoundaryValues(user, Q_loc, t)); 197 PetscCall(UpdateContextLabel(&user->time, t, user->op_ifunction, user->phys->solution_time_label)); 198 PetscCall(TSGetTimeStep(ts, &dt)); 199 PetscCall(UpdateContextLabel(&user->dt, dt, user->op_ifunction, user->phys->timestep_size_label)); 200 201 // Global-to-local 202 PetscCall(DMGlobalToLocalBegin(user->dm, Q, INSERT_VALUES, Q_loc)); 203 PetscCall(DMGlobalToLocalBegin(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc)); 204 PetscCall(DMGlobalToLocalEnd(user->dm, Q, INSERT_VALUES, Q_loc)); 205 PetscCall(DMGlobalToLocalEnd(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc)); 206 207 // Place PETSc vectors in CEED vectors 208 PetscCall(VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type)); 209 PetscCall(VecGetArrayReadAndMemType(Q_dot_loc, &q_dot, &q_dot_mem_type)); 210 PetscCall(VecGetArrayAndMemType(G_loc, &g, &g_mem_type)); 211 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, (PetscScalar *)q); 212 CeedVectorSetArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), CEED_USE_POINTER, (PetscScalar *)q_dot); 213 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 214 215 // Apply CEED operator 216 CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, CEED_REQUEST_IMMEDIATE); 217 218 // Restore vectors 219 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 220 CeedVectorTakeArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), NULL); 221 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 222 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, &q)); 223 PetscCall(VecRestoreArrayReadAndMemType(Q_dot_loc, &q_dot)); 224 PetscCall(VecRestoreArrayAndMemType(G_loc, &g)); 225 226 // Local-to-Global 227 PetscCall(VecZeroEntries(G)); 228 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G)); 229 230 // Restore vectors 231 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 232 233 PetscFunctionReturn(0); 234 } 235 236 static PetscErrorCode MatMult_NS_IJacobian(Mat J, Vec Q, Vec G) { 237 User user; 238 const PetscScalar *q; 239 PetscScalar *g; 240 PetscMemType q_mem_type, g_mem_type; 241 PetscFunctionBeginUser; 242 PetscCall(MatShellGetContext(J, &user)); 243 Vec Q_loc = user->Q_dot_loc, // Note - Q_dot_loc has zero BCs 244 G_loc; 245 246 // Get local vectors 247 PetscCall(DMGetLocalVector(user->dm, &G_loc)); 248 249 // Global-to-local 250 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 251 252 // Place PETSc vectors in CEED vectors 253 PetscCall(VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type)); 254 PetscCall(VecGetArrayAndMemType(G_loc, &g, &g_mem_type)); 255 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, (PetscScalar *)q); 256 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 257 258 // Apply CEED operator 259 CeedOperatorApply(user->op_ijacobian, user->q_ceed, user->g_ceed, CEED_REQUEST_IMMEDIATE); 260 261 // Restore vectors 262 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 263 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 264 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, &q)); 265 PetscCall(VecRestoreArrayAndMemType(G_loc, &g)); 266 267 // Local-to-Global 268 PetscCall(VecZeroEntries(G)); 269 PetscCall(DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G)); 270 271 // Restore vectors 272 PetscCall(DMRestoreLocalVector(user->dm, &G_loc)); 273 PetscFunctionReturn(0); 274 } 275 276 PetscErrorCode MatGetDiagonal_NS_IJacobian(Mat A, Vec D) { 277 User user; 278 Vec D_loc; 279 PetscScalar *d; 280 PetscMemType mem_type; 281 282 PetscFunctionBeginUser; 283 MatShellGetContext(A, &user); 284 PetscCall(DMGetLocalVector(user->dm, &D_loc)); 285 PetscCall(VecGetArrayAndMemType(D_loc, &d, &mem_type)); 286 CeedVectorSetArray(user->g_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, d); 287 CeedOperatorLinearAssembleDiagonal(user->op_ijacobian, user->g_ceed, CEED_REQUEST_IMMEDIATE); 288 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(mem_type), NULL); 289 PetscCall(VecRestoreArrayAndMemType(D_loc, &d)); 290 PetscCall(VecZeroEntries(D)); 291 PetscCall(DMLocalToGlobal(user->dm, D_loc, ADD_VALUES, D)); 292 PetscCall(DMRestoreLocalVector(user->dm, &D_loc)); 293 VecViewFromOptions(D, NULL, "-diag_vec_view"); 294 PetscFunctionReturn(0); 295 } 296 297 static PetscErrorCode FormPreallocation(User user, PetscBool pbdiagonal, Mat J, CeedVector *coo_values) { 298 PetscCount ncoo; 299 PetscInt *rows, *cols; 300 301 PetscFunctionBeginUser; 302 if (pbdiagonal) { 303 CeedSize l_size; 304 CeedOperatorGetActiveVectorLengths(user->op_ijacobian, &l_size, NULL); 305 ncoo = l_size * 5; 306 rows = malloc(ncoo * sizeof(rows[0])); 307 cols = malloc(ncoo * sizeof(cols[0])); 308 for (PetscCount n = 0; n < l_size / 5; n++) { 309 for (PetscInt i = 0; i < 5; i++) { 310 for (PetscInt j = 0; j < 5; j++) { 311 rows[(n * 5 + i) * 5 + j] = n * 5 + i; 312 cols[(n * 5 + i) * 5 + j] = n * 5 + j; 313 } 314 } 315 } 316 } else { 317 PetscCall(CeedOperatorLinearAssembleSymbolic(user->op_ijacobian, &ncoo, &rows, &cols)); 318 } 319 PetscCall(MatSetPreallocationCOOLocal(J, ncoo, rows, cols)); 320 free(rows); 321 free(cols); 322 CeedVectorCreate(user->ceed, ncoo, coo_values); 323 PetscFunctionReturn(0); 324 } 325 326 static PetscErrorCode FormSetValues(User user, PetscBool pbdiagonal, Mat J, CeedVector coo_values) { 327 CeedMemType mem_type = CEED_MEM_HOST; 328 const PetscScalar *values; 329 MatType mat_type; 330 331 PetscFunctionBeginUser; 332 PetscCall(MatGetType(J, &mat_type)); 333 if (strstr(mat_type, "kokkos") || strstr(mat_type, "cusparse")) mem_type = CEED_MEM_DEVICE; 334 if (user->app_ctx->pmat_pbdiagonal) { 335 CeedOperatorLinearAssemblePointBlockDiagonal(user->op_ijacobian, coo_values, CEED_REQUEST_IMMEDIATE); 336 } else { 337 CeedOperatorLinearAssemble(user->op_ijacobian, coo_values); 338 } 339 CeedVectorGetArrayRead(coo_values, mem_type, &values); 340 PetscCall(MatSetValuesCOO(J, values, INSERT_VALUES)); 341 CeedVectorRestoreArrayRead(coo_values, &values); 342 PetscFunctionReturn(0); 343 } 344 345 PetscErrorCode FormIJacobian_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, PetscReal shift, Mat J, Mat J_pre, void *user_data) { 346 User user = *(User *)user_data; 347 PetscBool J_is_shell, J_is_mffd, J_pre_is_shell; 348 PetscFunctionBeginUser; 349 if (user->phys->ijacobian_time_shift_label) CeedOperatorContextSetDouble(user->op_ijacobian, user->phys->ijacobian_time_shift_label, &shift); 350 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATMFFD, &J_is_mffd)); 351 PetscCall(PetscObjectTypeCompare((PetscObject)J, MATSHELL, &J_is_shell)); 352 PetscCall(PetscObjectTypeCompare((PetscObject)J_pre, MATSHELL, &J_pre_is_shell)); 353 if (!user->matrices_set_up) { 354 if (J_is_shell) { 355 PetscCall(MatShellSetContext(J, user)); 356 PetscCall(MatShellSetOperation(J, MATOP_MULT, (void (*)(void))MatMult_NS_IJacobian)); 357 PetscCall(MatShellSetOperation(J, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_NS_IJacobian)); 358 PetscCall(MatSetUp(J)); 359 } 360 if (!J_pre_is_shell) { 361 PetscCall(FormPreallocation(user, user->app_ctx->pmat_pbdiagonal, J_pre, &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, user->coo_values_pmat)); 370 } 371 if (user->coo_values_amat) { 372 PetscCall(FormSetValues(user, PETSC_FALSE, J, user->coo_values_amat)); 373 } else if (J_is_mffd) { 374 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 375 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 376 } 377 PetscFunctionReturn(0); 378 } 379 380 PetscErrorCode WriteOutput(User user, Vec Q, PetscInt step_no, PetscScalar time) { 381 Vec Q_loc; 382 char file_path[PETSC_MAX_PATH_LEN]; 383 PetscViewer viewer; 384 PetscFunctionBeginUser; 385 386 if (user->app_ctx->checkpoint_vtk) { 387 // Set up output 388 PetscCall(DMGetLocalVector(user->dm, &Q_loc)); 389 PetscCall(PetscObjectSetName((PetscObject)Q_loc, "StateVec")); 390 PetscCall(VecZeroEntries(Q_loc)); 391 PetscCall(DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc)); 392 393 // Output 394 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no)); 395 396 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, FILE_MODE_WRITE, &viewer)); 397 PetscCall(VecView(Q_loc, viewer)); 398 PetscCall(PetscViewerDestroy(&viewer)); 399 if (user->dm_viz) { 400 Vec Q_refined, Q_refined_loc; 401 char file_path_refined[PETSC_MAX_PATH_LEN]; 402 PetscViewer viewer_refined; 403 404 PetscCall(DMGetGlobalVector(user->dm_viz, &Q_refined)); 405 PetscCall(DMGetLocalVector(user->dm_viz, &Q_refined_loc)); 406 PetscCall(PetscObjectSetName((PetscObject)Q_refined_loc, "Refined")); 407 408 PetscCall(MatInterpolate(user->interp_viz, Q, Q_refined)); 409 PetscCall(VecZeroEntries(Q_refined_loc)); 410 PetscCall(DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc)); 411 412 PetscCall( 413 PetscSNPrintf(file_path_refined, sizeof file_path_refined, "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, step_no)); 414 415 PetscCall(PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), file_path_refined, FILE_MODE_WRITE, &viewer_refined)); 416 PetscCall(VecView(Q_refined_loc, viewer_refined)); 417 PetscCall(DMRestoreLocalVector(user->dm_viz, &Q_refined_loc)); 418 PetscCall(DMRestoreGlobalVector(user->dm_viz, &Q_refined)); 419 PetscCall(PetscViewerDestroy(&viewer_refined)); 420 } 421 PetscCall(DMRestoreLocalVector(user->dm, &Q_loc)); 422 } 423 424 // Save data in a binary file for continuation of simulations 425 if (user->app_ctx->add_stepnum2bin) { 426 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution-%" PetscInt_FMT ".bin", user->app_ctx->output_dir, step_no)); 427 } else { 428 PetscCall(PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", user->app_ctx->output_dir)); 429 } 430 PetscCall(PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer)); 431 432 PetscInt token = FLUIDS_FILE_TOKEN; 433 PetscCall(PetscViewerBinaryWrite(viewer, &token, 1, PETSC_INT)); 434 PetscCall(PetscViewerBinaryWrite(viewer, &step_no, 1, PETSC_INT)); 435 time /= user->units->second; // Dimensionalize time back 436 PetscCall(PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL)); 437 PetscCall(VecView(Q, viewer)); 438 PetscCall(PetscViewerDestroy(&viewer)); 439 PetscFunctionReturn(0); 440 } 441 442 // CSV Monitor 443 PetscErrorCode TSMonitor_WallForce(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) { 444 User user = ctx; 445 Vec G_loc; 446 PetscInt num_wall = user->app_ctx->wall_forces.num_wall, dim = 3; 447 const PetscInt *walls = user->app_ctx->wall_forces.walls; 448 PetscViewer viewer = user->app_ctx->wall_forces.viewer; 449 PetscViewerFormat format = user->app_ctx->wall_forces.viewer_format; 450 PetscScalar *reaction_force; 451 PetscBool iascii; 452 453 PetscFunctionBeginUser; 454 if (!viewer) PetscFunctionReturn(0); 455 PetscCall(DMGetNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 456 PetscCall(PetscMalloc1(num_wall * dim, &reaction_force)); 457 PetscCall(Surface_Forces_NS(user->dm, G_loc, num_wall, walls, reaction_force)); 458 PetscCall(DMRestoreNamedLocalVector(user->dm, "ResidualLocal", &G_loc)); 459 460 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 461 462 if (iascii) { 463 if (format == PETSC_VIEWER_ASCII_CSV && !user->app_ctx->wall_forces.header_written) { 464 PetscCall(PetscViewerASCIIPrintf(viewer, "Step,Time,Wall,ForceX,ForceY,ForceZ\n")); 465 user->app_ctx->wall_forces.header_written = PETSC_TRUE; 466 } 467 for (PetscInt w = 0; w < num_wall; w++) { 468 PetscInt wall = walls[w]; 469 if (format == PETSC_VIEWER_ASCII_CSV) { 470 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ",%g,%" PetscInt_FMT ",%g,%g,%g\n", step_no, time, wall, 471 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2])); 472 473 } else { 474 PetscCall(PetscViewerASCIIPrintf(viewer, "Wall %" PetscInt_FMT " Forces: Force_x = %12g, Force_y = %12g, Force_z = %12g\n", wall, 475 reaction_force[w * dim + 0], reaction_force[w * dim + 1], reaction_force[w * dim + 2])); 476 } 477 } 478 } 479 PetscCall(PetscFree(reaction_force)); 480 PetscFunctionReturn(0); 481 } 482 483 // User provided TS Monitor 484 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, Vec Q, void *ctx) { 485 User user = ctx; 486 PetscFunctionBeginUser; 487 488 // Print every 'checkpoint_interval' steps 489 if (user->app_ctx->checkpoint_interval <= 0 || step_no % user->app_ctx->checkpoint_interval != 0 || 490 (user->app_ctx->cont_steps == step_no && step_no != 0)) 491 PetscFunctionReturn(0); 492 493 PetscCall(WriteOutput(user, Q, step_no, time)); 494 495 PetscFunctionReturn(0); 496 } 497 498 // TS: Create, setup, and solve 499 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, Vec *Q, PetscScalar *f_time, TS *ts) { 500 MPI_Comm comm = user->comm; 501 TSAdapt adapt; 502 PetscScalar final_time; 503 PetscFunctionBeginUser; 504 505 PetscCall(TSCreate(comm, ts)); 506 PetscCall(TSSetDM(*ts, dm)); 507 if (phys->implicit) { 508 PetscCall(TSSetType(*ts, TSBDF)); 509 if (user->op_ifunction) { 510 PetscCall(TSSetIFunction(*ts, NULL, IFunction_NS, &user)); 511 } else { // Implicit integrators can fall back to using an RHSFunction 512 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user)); 513 } 514 if (user->op_ijacobian) { 515 PetscCall(DMTSSetIJacobian(dm, FormIJacobian_NS, &user)); 516 if (app_ctx->amat_type) { 517 Mat Pmat, Amat; 518 PetscCall(DMCreateMatrix(dm, &Pmat)); 519 PetscCall(DMSetMatType(dm, app_ctx->amat_type)); 520 PetscCall(DMCreateMatrix(dm, &Amat)); 521 PetscCall(TSSetIJacobian(*ts, Amat, Pmat, NULL, NULL)); 522 PetscCall(MatDestroy(&Amat)); 523 PetscCall(MatDestroy(&Pmat)); 524 } 525 } 526 } else { 527 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, "Problem does not provide RHSFunction"); 528 PetscCall(TSSetType(*ts, TSRK)); 529 PetscCall(TSRKSetType(*ts, TSRK5F)); 530 PetscCall(TSSetRHSFunction(*ts, NULL, RHS_NS, &user)); 531 } 532 PetscCall(TSSetMaxTime(*ts, 500. * user->units->second)); 533 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER)); 534 PetscCall(TSSetErrorIfStepFails(*ts, PETSC_FALSE)); 535 PetscCall(TSSetTimeStep(*ts, 1.e-2 * user->units->second)); 536 if (app_ctx->test_type != TESTTYPE_NONE) { 537 PetscCall(TSSetMaxSteps(*ts, 10)); 538 } 539 PetscCall(TSGetAdapt(*ts, &adapt)); 540 PetscCall(TSAdaptSetStepLimits(adapt, 1.e-12 * user->units->second, 1.e2 * user->units->second)); 541 PetscCall(TSSetFromOptions(*ts)); 542 user->time = user->time_bc_set = -1.0; // require all BCs and ctx to be updated 543 user->dt = -1.0; 544 if (!app_ctx->cont_steps) { // print initial condition 545 if (app_ctx->test_type == TESTTYPE_NONE) { 546 PetscCall(TSMonitor_NS(*ts, 0, 0., *Q, user)); 547 } 548 } else { // continue from time of last output 549 PetscInt count; 550 PetscViewer viewer; 551 552 if (app_ctx->cont_time <= 0) { // Legacy files did not include step number and time 553 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_time_file, FILE_MODE_READ, &viewer)); 554 PetscCall(PetscViewerBinaryRead(viewer, &app_ctx->cont_time, 1, &count, PETSC_REAL)); 555 PetscCall(PetscViewerDestroy(&viewer)); 556 PetscCheck(app_ctx->cont_steps != -1, comm, PETSC_ERR_ARG_INCOMP, 557 "-continue step number not specified, but checkpoint file does not contain a step number (likely written by older code version)"); 558 } 559 PetscCall(TSSetTime(*ts, app_ctx->cont_time * user->units->second)); 560 PetscCall(TSSetStepNumber(*ts, app_ctx->cont_steps)); 561 } 562 if (app_ctx->test_type == TESTTYPE_NONE) { 563 PetscCall(TSMonitorSet(*ts, TSMonitor_NS, user, NULL)); 564 } 565 if (app_ctx->wall_forces.viewer) { 566 PetscCall(TSMonitorSet(*ts, TSMonitor_WallForce, user, NULL)); 567 } 568 if (app_ctx->turb_spanstats_enable) { 569 PetscCall(TSMonitorSet(*ts, TSMonitor_Statistics, user, NULL)); 570 CeedScalar previous_time = app_ctx->cont_time * user->units->second; 571 CeedOperatorContextSetDouble(user->spanstats.op_stats_collect, user->spanstats.previous_time_label, &previous_time); 572 } 573 574 // Solve 575 PetscReal start_time; 576 PetscInt start_step; 577 PetscCall(TSGetTime(*ts, &start_time)); 578 PetscCall(TSGetStepNumber(*ts, &start_step)); 579 580 PetscPreLoadBegin(PETSC_FALSE, "Fluids Solve"); 581 PetscCall(TSSetTime(*ts, start_time)); 582 PetscCall(TSSetStepNumber(*ts, start_step)); 583 if (PetscPreLoadingOn) { 584 // LCOV_EXCL_START 585 SNES snes; 586 Vec Q_preload; 587 PetscReal rtol; 588 PetscCall(VecDuplicate(*Q, &Q_preload)); 589 PetscCall(VecCopy(*Q, Q_preload)); 590 PetscCall(TSGetSNES(*ts, &snes)); 591 PetscCall(SNESGetTolerances(snes, NULL, &rtol, NULL, NULL, NULL)); 592 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, .99, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 593 PetscCall(TSSetSolution(*ts, *Q)); 594 PetscCall(TSStep(*ts)); 595 PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 596 PetscCall(VecDestroy(&Q_preload)); 597 // LCOV_EXCL_STOP 598 } else { 599 PetscCall(PetscBarrier((PetscObject)*ts)); 600 PetscCall(TSSolve(*ts, *Q)); 601 } 602 PetscPreLoadEnd(); 603 604 PetscCall(TSGetSolveTime(*ts, &final_time)); 605 *f_time = final_time; 606 607 if (app_ctx->test_type == TESTTYPE_NONE) { 608 PetscInt step_no; 609 PetscCall(TSGetStepNumber(*ts, &step_no)); 610 if (user->app_ctx->checkpoint_interval > 0 || user->app_ctx->checkpoint_interval == -1) { 611 PetscCall(WriteOutput(user, *Q, step_no, final_time)); 612 } 613 614 PetscLogEvent stage_id; 615 PetscStageLog stage_log; 616 617 PetscCall(PetscLogStageGetId("Fluids Solve", &stage_id)); 618 PetscCall(PetscLogGetStageLog(&stage_log)); 619 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time taken for solution (sec): %g\n", stage_log->stageInfo[stage_id].perfInfo.time)); 620 } 621 PetscFunctionReturn(0); 622 } 623