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