1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Time-stepping functions for Navier-Stokes example using PETSc 19 20 #include "../navierstokes.h" 21 #include "../qfunctions/mass.h" 22 23 // Compute mass matrix for explicit scheme 24 PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, CeedData ceed_data, 25 Vec M) { 26 Vec M_loc; 27 CeedQFunction qf_mass; 28 CeedOperator op_mass; 29 CeedVector m_ceed, ones_vec; 30 CeedInt num_comp_q, q_data_size; 31 PetscErrorCode ierr; 32 PetscFunctionBeginUser; 33 34 // CEED Restriction 35 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q); 36 CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size); 37 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &m_ceed, NULL); 38 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &ones_vec, NULL); 39 CeedVectorSetValue(ones_vec, 1.0); 40 41 // CEED QFunction 42 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 43 CeedQFunctionAddInput(qf_mass, "q", num_comp_q, CEED_EVAL_INTERP); 44 CeedQFunctionAddInput(qf_mass, "q_data", q_data_size, CEED_EVAL_NONE); 45 CeedQFunctionAddOutput(qf_mass, "v", num_comp_q, CEED_EVAL_INTERP); 46 47 // CEED Operator 48 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 49 CeedOperatorSetField(op_mass, "q", ceed_data->elem_restr_q, ceed_data->basis_q, 50 CEED_VECTOR_ACTIVE); 51 CeedOperatorSetField(op_mass, "q_data", ceed_data->elem_restr_qd_i, 52 CEED_BASIS_COLLOCATED, ceed_data->q_data); 53 CeedOperatorSetField(op_mass, "v", ceed_data->elem_restr_q, ceed_data->basis_q, 54 CEED_VECTOR_ACTIVE); 55 56 // Place PETSc vector in CEED vector 57 CeedScalar *m; 58 PetscMemType m_mem_type; 59 ierr = DMGetLocalVector(dm, &M_loc); CHKERRQ(ierr); 60 ierr = VecGetArrayAndMemType(M_loc, (PetscScalar **)&m, &m_mem_type); 61 CHKERRQ(ierr); 62 CeedVectorSetArray(m_ceed, MemTypeP2C(m_mem_type), CEED_USE_POINTER, m); 63 64 // Apply CEED Operator 65 CeedOperatorApply(op_mass, ones_vec, m_ceed, CEED_REQUEST_IMMEDIATE); 66 67 // Restore vectors 68 CeedVectorTakeArray(m_ceed, MemTypeP2C(m_mem_type), NULL); 69 ierr = VecRestoreArrayReadAndMemType(M_loc, (const PetscScalar **)&m); 70 CHKERRQ(ierr); 71 72 // Local-to-Global 73 ierr = VecZeroEntries(M); CHKERRQ(ierr); 74 ierr = DMLocalToGlobal(dm, M_loc, ADD_VALUES, M); CHKERRQ(ierr); 75 ierr = DMRestoreLocalVector(dm, &M_loc); CHKERRQ(ierr); 76 77 // Invert diagonally lumped mass vector for RHS function 78 ierr = VecReciprocal(M); CHKERRQ(ierr); 79 80 // Cleanup 81 CeedVectorDestroy(&ones_vec); 82 CeedVectorDestroy(&m_ceed); 83 CeedQFunctionDestroy(&qf_mass); 84 CeedOperatorDestroy(&op_mass); 85 86 PetscFunctionReturn(0); 87 } 88 89 // RHS (Explicit time-stepper) function setup 90 // This is the RHS of the ODE, given as u_t = G(t,u) 91 // This function takes in a state vector Q and writes into G 92 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data) { 93 94 User user = *(User *)user_data; 95 PetscScalar *q, *g; 96 Vec Q_loc, G_loc; 97 PetscMemType q_mem_type, g_mem_type; 98 PetscErrorCode ierr; 99 PetscFunctionBeginUser; 100 101 // Update EulerContext 102 if (user->phys->has_curr_time) user->phys->euler_ctx->curr_time = t; 103 104 // Get local vectors 105 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 106 ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 107 108 // Global-to-local 109 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 110 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 111 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, 0.0, 112 NULL, NULL, NULL); CHKERRQ(ierr); 113 ierr = VecZeroEntries(G_loc); CHKERRQ(ierr); 114 115 // Place PETSc vectors in CEED vectors 116 ierr = VecGetArrayReadAndMemType(Q_loc, (const PetscScalar **)&q, &q_mem_type); 117 CHKERRQ(ierr); 118 ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr); 119 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, q); 120 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 121 122 // Apply CEED operator 123 CeedOperatorApply(user->op_rhs, user->q_ceed, user->g_ceed, 124 CEED_REQUEST_IMMEDIATE); 125 126 // Restore vectors 127 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 128 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 129 ierr = VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q); 130 CHKERRQ(ierr); 131 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 132 133 // Local-to-Global 134 ierr = VecZeroEntries(G); CHKERRQ(ierr); 135 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 136 137 // Inverse of the lumped mass matrix (M is Minv) 138 ierr = VecPointwiseMult(G, G, user->M); CHKERRQ(ierr); 139 140 // Restore vectors 141 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 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, Q_dot_loc, G_loc; 154 PetscMemType q_mem_type, q_dot_mem_type, g_mem_type; 155 PetscErrorCode ierr; 156 PetscFunctionBeginUser; 157 158 // Update EulerContext 159 if (user->phys->has_curr_time) user->phys->euler_ctx->curr_time = t; 160 161 // Get local vectors 162 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 163 ierr = DMGetLocalVector(user->dm, &Q_dot_loc); CHKERRQ(ierr); 164 ierr = DMGetLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 165 166 // Global-to-local 167 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 168 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 169 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Q_loc, 0.0, 170 NULL, NULL, NULL); CHKERRQ(ierr); 171 ierr = VecZeroEntries(Q_dot_loc); CHKERRQ(ierr); 172 ierr = DMGlobalToLocal(user->dm, Q_dot, INSERT_VALUES, Q_dot_loc); 173 CHKERRQ(ierr); 174 ierr = VecZeroEntries(G_loc); CHKERRQ(ierr); 175 176 // Place PETSc vectors in CEED vectors 177 ierr = VecGetArrayReadAndMemType(Q_loc, &q, &q_mem_type); CHKERRQ(ierr); 178 ierr = VecGetArrayReadAndMemType(Q_dot_loc, &q_dot, &q_dot_mem_type); 179 CHKERRQ(ierr); 180 ierr = VecGetArrayAndMemType(G_loc, &g, &g_mem_type); CHKERRQ(ierr); 181 CeedVectorSetArray(user->q_ceed, MemTypeP2C(q_mem_type), CEED_USE_POINTER, 182 (PetscScalar *)q); 183 CeedVectorSetArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), 184 CEED_USE_POINTER, 185 (PetscScalar *)q_dot); 186 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 187 188 // Apply CEED operator 189 CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, 190 CEED_REQUEST_IMMEDIATE); 191 192 // Restore vectors 193 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 194 CeedVectorTakeArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), NULL); 195 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 196 ierr = VecRestoreArrayReadAndMemType(Q_loc, &q); CHKERRQ(ierr); 197 ierr = VecRestoreArrayReadAndMemType(Q_dot_loc, &q_dot); CHKERRQ(ierr); 198 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 199 200 // Local-to-Global 201 ierr = VecZeroEntries(G); CHKERRQ(ierr); 202 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 203 204 // Restore vectors 205 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 206 ierr = DMRestoreLocalVector(user->dm, &Q_dot_loc); CHKERRQ(ierr); 207 ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 208 209 PetscFunctionReturn(0); 210 } 211 212 // User provided TS Monitor 213 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, 214 Vec Q, void *ctx) { 215 User user = ctx; 216 Vec Q_loc; 217 char file_path[PETSC_MAX_PATH_LEN]; 218 PetscViewer viewer; 219 PetscErrorCode ierr; 220 PetscFunctionBeginUser; 221 222 // Print every 'output_freq' steps 223 if (step_no % user->app_ctx->output_freq != 0) 224 PetscFunctionReturn(0); 225 226 // Set up output 227 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 228 ierr = PetscObjectSetName((PetscObject)Q_loc, "StateVec"); CHKERRQ(ierr); 229 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 230 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 231 232 // Output 233 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03D.vtu", 234 user->app_ctx->output_dir, step_no + user->app_ctx->cont_steps); 235 CHKERRQ(ierr); 236 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, 237 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 238 ierr = VecView(Q_loc, viewer); CHKERRQ(ierr); 239 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 240 if (user->dm_viz) { 241 Vec Q_refined, Q_refined_loc; 242 char file_path_refined[PETSC_MAX_PATH_LEN]; 243 PetscViewer viewer_refined; 244 245 ierr = DMGetGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 246 ierr = DMGetLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 247 ierr = PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"); 248 CHKERRQ(ierr); 249 ierr = MatInterpolate(user->interp_viz, Q, Q_refined); CHKERRQ(ierr); 250 ierr = VecZeroEntries(Q_refined_loc); CHKERRQ(ierr); 251 ierr = DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc); 252 CHKERRQ(ierr); 253 ierr = PetscSNPrintf(file_path_refined, sizeof file_path_refined, 254 "%s/nsrefined-%03D.vtu", 255 user->app_ctx->output_dir, step_no + user->app_ctx->cont_steps); 256 CHKERRQ(ierr); 257 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), 258 file_path_refined, 259 FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 260 ierr = VecView(Q_refined_loc, viewer_refined); CHKERRQ(ierr); 261 ierr = DMRestoreLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 262 ierr = DMRestoreGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 263 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 264 } 265 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 266 267 // Save data in a binary file for continuation of simulations 268 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", 269 user->app_ctx->output_dir); CHKERRQ(ierr); 270 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 271 CHKERRQ(ierr); 272 ierr = VecView(Q, viewer); CHKERRQ(ierr); 273 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 274 275 // Save time stamp 276 // Dimensionalize time back 277 time /= user->units->second; 278 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 279 user->app_ctx->output_dir); CHKERRQ(ierr); 280 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 281 CHKERRQ(ierr); 282 #if PETSC_VERSION_GE(3,13,0) 283 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 284 #else 285 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 286 #endif 287 CHKERRQ(ierr); 288 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 289 290 PetscFunctionReturn(0); 291 } 292 293 // TS: Create, setup, and solve 294 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, 295 Vec *Q, PetscScalar *f_time, TS *ts) { 296 297 MPI_Comm comm = user->comm; 298 TSAdapt adapt; 299 PetscScalar final_time; 300 PetscErrorCode ierr; 301 PetscFunctionBeginUser; 302 303 ierr = TSCreate(comm, ts); CHKERRQ(ierr); 304 ierr = TSSetDM(*ts, dm); CHKERRQ(ierr); 305 if (phys->implicit) { 306 ierr = TSSetType(*ts, TSBDF); CHKERRQ(ierr); 307 if (user->op_ifunction) { 308 ierr = TSSetIFunction(*ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 309 } else { // Implicit integrators can fall back to using an RHSFunction 310 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 311 } 312 } else { 313 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 314 "Problem does not provide RHSFunction"); 315 ierr = TSSetType(*ts, TSRK); CHKERRQ(ierr); 316 ierr = TSRKSetType(*ts, TSRK5F); CHKERRQ(ierr); 317 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 318 } 319 ierr = TSSetMaxTime(*ts, 500. * user->units->second); CHKERRQ(ierr); 320 ierr = TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 321 ierr = TSSetTimeStep(*ts, 1.e-2 * user->units->second); CHKERRQ(ierr); 322 if (app_ctx->test_mode) {ierr = TSSetMaxSteps(*ts, 10); CHKERRQ(ierr);} 323 ierr = TSGetAdapt(*ts, &adapt); CHKERRQ(ierr); 324 ierr = TSAdaptSetStepLimits(adapt, 325 1.e-12 * user->units->second, 326 1.e2 * user->units->second); CHKERRQ(ierr); 327 ierr = TSSetFromOptions(*ts); CHKERRQ(ierr); 328 if (!app_ctx->cont_steps) { // print initial condition 329 if (!app_ctx->test_mode) { 330 ierr = TSMonitor_NS(*ts, 0, 0., *Q, user); CHKERRQ(ierr); 331 } 332 } else { // continue from time of last output 333 PetscReal time; 334 PetscInt count; 335 PetscViewer viewer; 336 char file_path[PETSC_MAX_PATH_LEN]; 337 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 338 app_ctx->output_dir); CHKERRQ(ierr); 339 ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer); 340 CHKERRQ(ierr); 341 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 342 CHKERRQ(ierr); 343 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 344 ierr = TSSetTime(*ts, time * user->units->second); CHKERRQ(ierr); 345 } 346 if (!app_ctx->test_mode) { 347 ierr = TSMonitorSet(*ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 348 } 349 350 // Solve 351 double start, cpu_time_used; 352 start = MPI_Wtime(); 353 ierr = PetscBarrier((PetscObject) *ts); CHKERRQ(ierr); 354 ierr = TSSolve(*ts, *Q); CHKERRQ(ierr); 355 cpu_time_used = MPI_Wtime() - start; 356 ierr = TSGetSolveTime(*ts, &final_time); CHKERRQ(ierr); 357 *f_time = final_time; 358 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 359 comm); CHKERRQ(ierr); 360 if (!app_ctx->test_mode) { 361 ierr = PetscPrintf(PETSC_COMM_WORLD, 362 "Time taken for solution (sec): %g\n", 363 (double)cpu_time_used); CHKERRQ(ierr); 364 } 365 PetscFunctionReturn(0); 366 } 367