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 solution time 93 if (user->phys->solution_time_label) 94 CeedOperatorContextSetDouble(user->op_rhs, user->phys->solution_time_label, &t); 95 96 // Get local vectors 97 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 98 ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 99 100 // Global-to-local 101 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 102 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 103 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, 0.0, 104 NULL, NULL, NULL); CHKERRQ(ierr); 105 ierr = VecZeroEntries(G_loc); CHKERRQ(ierr); 106 107 // Place PETSc vectors in CEED vectors 108 ierr = VecGetArrayReadAndMemType(Q_loc, (const PetscScalar **)&q, &q_mem_type); 109 CHKERRQ(ierr); 110 ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr); 111 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, q); 112 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 113 114 // Apply CEED operator 115 CeedOperatorApply(user->op_rhs, user->q_ceed, user->g_ceed, 116 CEED_REQUEST_IMMEDIATE); 117 118 // Restore vectors 119 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 120 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 121 ierr = VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q); 122 CHKERRQ(ierr); 123 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 124 125 // Local-to-Global 126 ierr = VecZeroEntries(G); CHKERRQ(ierr); 127 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 128 129 // Inverse of the lumped mass matrix (M is Minv) 130 ierr = VecPointwiseMult(G, G, user->M); CHKERRQ(ierr); 131 132 // Restore vectors 133 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 134 ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 135 136 PetscFunctionReturn(0); 137 } 138 139 // Implicit time-stepper function setup 140 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, 141 void *user_data) { 142 User user = *(User *)user_data; 143 const PetscScalar *q, *q_dot; 144 PetscScalar *g; 145 Vec Q_loc, Q_dot_loc, G_loc; 146 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type; 147 PetscErrorCode ierr; 148 PetscFunctionBeginUser; 149 150 // Update solution time 151 if (user->phys->solution_time_label) 152 CeedOperatorContextSetDouble(user->op_ifunction, 153 user->phys->solution_time_label, &t); 154 155 // Get local vectors 156 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 157 ierr = DMGetLocalVector(user->dm, &Q_dot_loc); CHKERRQ(ierr); 158 ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 159 160 // Global-to-local 161 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 162 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 163 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, 0.0, 164 NULL, NULL, NULL); CHKERRQ(ierr); 165 ierr = VecZeroEntries(Q_dot_loc); CHKERRQ(ierr); 166 ierr = DMGlobalToLocal(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc); 167 CHKERRQ(ierr); 168 ierr = VecZeroEntries(G_loc); CHKERRQ(ierr); 169 170 // Place PETSc vectors in CEED vectors 171 ierr = VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type); CHKERRQ(ierr); 172 ierr = VecGetArrayReadAndMemType(Q_dot_loc, &q_dot, &q_dot_mem_type); 173 CHKERRQ(ierr); 174 ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr); 175 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, 176 (PetscScalar *)q); 177 CeedVectorSetArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), 178 CEED_USE_POINTER, (PetscScalar *)q_dot); 179 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 180 181 // Apply CEED operator 182 CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, 183 CEED_REQUEST_IMMEDIATE); 184 185 // Restore vectors 186 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 187 CeedVectorTakeArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), NULL); 188 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 189 ierr = VecRestoreArrayReadAndMemType(Q_loc, &q); CHKERRQ(ierr); 190 ierr = VecRestoreArrayReadAndMemType(Q_dot_loc, &q_dot); CHKERRQ(ierr); 191 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 192 193 // Local-to-Global 194 ierr = VecZeroEntries(G); CHKERRQ(ierr); 195 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 196 197 // Restore vectors 198 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 199 ierr = DMRestoreLocalVector(user->dm, &Q_dot_loc); CHKERRQ(ierr); 200 ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 201 202 PetscFunctionReturn(0); 203 } 204 205 // User provided TS Monitor 206 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, 207 Vec Q, void *ctx) { 208 User user = ctx; 209 Vec Q_loc; 210 char file_path[PETSC_MAX_PATH_LEN]; 211 PetscViewer viewer; 212 PetscErrorCode ierr; 213 PetscFunctionBeginUser; 214 215 // Print every 'output_freq' steps 216 if (step_no % user->app_ctx->output_freq != 0) 217 PetscFunctionReturn(0); 218 219 // Set up output 220 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 221 ierr = PetscObjectSetName((PetscObject)Q_loc, "StateVec"); CHKERRQ(ierr); 222 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 223 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 224 225 // Output 226 ierr = PetscSNPrintf(file_path, sizeof file_path, 227 "%s/ns-%03" PetscInt_FMT ".vtu", 228 user->app_ctx->output_dir, step_no + user->app_ctx->cont_steps); 229 CHKERRQ(ierr); 230 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, 231 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 232 ierr = VecView(Q_loc, viewer); CHKERRQ(ierr); 233 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 234 if (user->dm_viz) { 235 Vec Q_refined, Q_refined_loc; 236 char file_path_refined[PETSC_MAX_PATH_LEN]; 237 PetscViewer viewer_refined; 238 239 ierr = DMGetGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 240 ierr = DMGetLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 241 ierr = PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"); 242 CHKERRQ(ierr); 243 ierr = MatInterpolate(user->interp_viz, Q, Q_refined); CHKERRQ(ierr); 244 ierr = VecZeroEntries(Q_refined_loc); CHKERRQ(ierr); 245 ierr = DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc); 246 CHKERRQ(ierr); 247 ierr = PetscSNPrintf(file_path_refined, sizeof file_path_refined, 248 "%s/nsrefined-%03" PetscInt_FMT ".vtu", user->app_ctx->output_dir, 249 step_no + user->app_ctx->cont_steps); 250 CHKERRQ(ierr); 251 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), 252 file_path_refined, FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 253 ierr = VecView(Q_refined_loc, viewer_refined); CHKERRQ(ierr); 254 ierr = DMRestoreLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 255 ierr = DMRestoreGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 256 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 257 } 258 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 259 260 // Save data in a binary file for continuation of simulations 261 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", 262 user->app_ctx->output_dir); CHKERRQ(ierr); 263 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 264 CHKERRQ(ierr); 265 ierr = VecView(Q, viewer); CHKERRQ(ierr); 266 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 267 268 // Save time stamp 269 // Dimensionalize time back 270 time /= user->units->second; 271 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 272 user->app_ctx->output_dir); CHKERRQ(ierr); 273 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 274 CHKERRQ(ierr); 275 #if PETSC_VERSION_GE(3,13,0) 276 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 277 #else 278 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 279 #endif 280 CHKERRQ(ierr); 281 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 282 283 PetscFunctionReturn(0); 284 } 285 286 // TS: Create, setup, and solve 287 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, 288 Vec *Q, PetscScalar *f_time, TS *ts) { 289 MPI_Comm comm = user->comm; 290 TSAdapt adapt; 291 PetscScalar final_time; 292 PetscErrorCode ierr; 293 PetscFunctionBeginUser; 294 295 ierr = TSCreate(comm, ts); CHKERRQ(ierr); 296 ierr = TSSetDM(*ts, dm); CHKERRQ(ierr); 297 if (phys->implicit) { 298 ierr = TSSetType(*ts, TSBDF); CHKERRQ(ierr); 299 if (user->op_ifunction) { 300 ierr = TSSetIFunction(*ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 301 } else { // Implicit integrators can fall back to using an RHSFunction 302 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 303 } 304 } else { 305 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 306 "Problem does not provide RHSFunction"); 307 ierr = TSSetType(*ts, TSRK); CHKERRQ(ierr); 308 ierr = TSRKSetType(*ts, TSRK5F); CHKERRQ(ierr); 309 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 310 } 311 ierr = TSSetMaxTime(*ts, 500. * user->units->second); CHKERRQ(ierr); 312 ierr = TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 313 ierr = TSSetTimeStep(*ts, 1.e-2 * user->units->second); CHKERRQ(ierr); 314 if (app_ctx->test_mode) {ierr = TSSetMaxSteps(*ts, 10); CHKERRQ(ierr);} 315 ierr = TSGetAdapt(*ts, &adapt); CHKERRQ(ierr); 316 ierr = TSAdaptSetStepLimits(adapt, 317 1.e-12 * user->units->second, 318 1.e2 * user->units->second); CHKERRQ(ierr); 319 ierr = TSSetFromOptions(*ts); CHKERRQ(ierr); 320 if (!app_ctx->cont_steps) { // print initial condition 321 if (!app_ctx->test_mode) { 322 ierr = TSMonitor_NS(*ts, 0, 0., *Q, user); CHKERRQ(ierr); 323 } 324 } else { // continue from time of last output 325 PetscReal time; 326 PetscInt count; 327 PetscViewer viewer; 328 char file_path[PETSC_MAX_PATH_LEN]; 329 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 330 app_ctx->output_dir); CHKERRQ(ierr); 331 ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer); 332 CHKERRQ(ierr); 333 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 334 CHKERRQ(ierr); 335 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 336 ierr = TSSetTime(*ts, time * user->units->second); CHKERRQ(ierr); 337 } 338 if (!app_ctx->test_mode) { 339 ierr = TSMonitorSet(*ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 340 } 341 342 // Solve 343 double start, cpu_time_used; 344 start = MPI_Wtime(); 345 ierr = PetscBarrier((PetscObject) *ts); CHKERRQ(ierr); 346 ierr = TSSolve(*ts, *Q); CHKERRQ(ierr); 347 cpu_time_used = MPI_Wtime() - start; 348 ierr = TSGetSolveTime(*ts, &final_time); CHKERRQ(ierr); 349 *f_time = final_time; 350 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 351 comm); CHKERRQ(ierr); 352 if (!app_ctx->test_mode) { 353 ierr = PetscPrintf(PETSC_COMM_WORLD, 354 "Time taken for solution (sec): %g\n", 355 (double)cpu_time_used); CHKERRQ(ierr); 356 } 357 PetscFunctionReturn(0); 358 } 359