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, "qdata", 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, "qdata", 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, (PetscScalar *)q_dot); 185 CeedVectorSetArray(user->g_ceed, MemTypeP2C(g_mem_type), CEED_USE_POINTER, g); 186 187 // Apply CEED operator 188 CeedOperatorApply(user->op_ifunction, user->q_ceed, user->g_ceed, 189 CEED_REQUEST_IMMEDIATE); 190 191 // Restore vectors 192 CeedVectorTakeArray(user->q_ceed, MemTypeP2C(q_mem_type), NULL); 193 CeedVectorTakeArray(user->q_dot_ceed, MemTypeP2C(q_dot_mem_type), NULL); 194 CeedVectorTakeArray(user->g_ceed, MemTypeP2C(g_mem_type), NULL); 195 ierr = VecRestoreArrayReadAndMemType(Q_loc, &q); CHKERRQ(ierr); 196 ierr = VecRestoreArrayReadAndMemType(Q_dot_loc, &q_dot); CHKERRQ(ierr); 197 ierr = VecRestoreArrayAndMemType(G_loc, &g); CHKERRQ(ierr); 198 199 // Local-to-Global 200 ierr = VecZeroEntries(G); CHKERRQ(ierr); 201 ierr = DMLocalToGlobal(user->dm, G_loc, ADD_VALUES, G); CHKERRQ(ierr); 202 203 // Restore vectors 204 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 205 ierr = DMRestoreLocalVector(user->dm, &Q_dot_loc); CHKERRQ(ierr); 206 ierr = DMRestoreLocalVector(user->dm, &G_loc); CHKERRQ(ierr); 207 208 PetscFunctionReturn(0); 209 } 210 211 // User provided TS Monitor 212 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, 213 Vec Q, void *ctx) { 214 User user = ctx; 215 Vec Q_loc; 216 char file_path[PETSC_MAX_PATH_LEN]; 217 PetscViewer viewer; 218 PetscErrorCode ierr; 219 PetscFunctionBeginUser; 220 221 // Print every 'output_freq' steps 222 if (step_no % user->app_ctx->output_freq != 0) 223 PetscFunctionReturn(0); 224 225 // Set up output 226 ierr = DMGetLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 227 ierr = PetscObjectSetName((PetscObject)Q_loc, "StateVec"); CHKERRQ(ierr); 228 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 229 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 230 231 // Output 232 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-%03D.vtu", 233 user->app_ctx->output_dir, step_no + user->app_ctx->cont_steps); 234 CHKERRQ(ierr); 235 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), file_path, 236 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 237 ierr = VecView(Q_loc, viewer); CHKERRQ(ierr); 238 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 239 if (user->dm_viz) { 240 Vec Q_refined, Q_refined_loc; 241 char file_path_refined[PETSC_MAX_PATH_LEN]; 242 PetscViewer viewer_refined; 243 244 ierr = DMGetGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 245 ierr = DMGetLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 246 ierr = PetscObjectSetName((PetscObject)Q_refined_loc, "Refined"); 247 CHKERRQ(ierr); 248 ierr = MatInterpolate(user->interp_viz, Q, Q_refined); CHKERRQ(ierr); 249 ierr = VecZeroEntries(Q_refined_loc); CHKERRQ(ierr); 250 ierr = DMGlobalToLocal(user->dm_viz, Q_refined, INSERT_VALUES, Q_refined_loc); 251 CHKERRQ(ierr); 252 ierr = PetscSNPrintf(file_path_refined, sizeof file_path_refined, 253 "%s/nsrefined-%03D.vtu", user->app_ctx->output_dir, 254 step_no + user->app_ctx->cont_steps); 255 CHKERRQ(ierr); 256 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q_refined), 257 file_path_refined, FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 258 ierr = VecView(Q_refined_loc, viewer_refined); CHKERRQ(ierr); 259 ierr = DMRestoreLocalVector(user->dm_viz, &Q_refined_loc); CHKERRQ(ierr); 260 ierr = DMRestoreGlobalVector(user->dm_viz, &Q_refined); CHKERRQ(ierr); 261 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 262 } 263 ierr = DMRestoreLocalVector(user->dm, &Q_loc); CHKERRQ(ierr); 264 265 // Save data in a binary file for continuation of simulations 266 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", 267 user->app_ctx->output_dir); CHKERRQ(ierr); 268 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 269 CHKERRQ(ierr); 270 ierr = VecView(Q, viewer); CHKERRQ(ierr); 271 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 272 273 // Save time stamp 274 // Dimensionalize time back 275 time /= user->units->second; 276 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 277 user->app_ctx->output_dir); CHKERRQ(ierr); 278 ierr = PetscViewerBinaryOpen(user->comm, file_path, FILE_MODE_WRITE, &viewer); 279 CHKERRQ(ierr); 280 #if PETSC_VERSION_GE(3,13,0) 281 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 282 #else 283 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 284 #endif 285 CHKERRQ(ierr); 286 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 287 288 PetscFunctionReturn(0); 289 } 290 291 // TS: Create, setup, and solve 292 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, 293 Vec *Q, PetscScalar *f_time, TS *ts) { 294 MPI_Comm comm = user->comm; 295 TSAdapt adapt; 296 PetscScalar final_time; 297 PetscErrorCode ierr; 298 PetscFunctionBeginUser; 299 300 ierr = TSCreate(comm, ts); CHKERRQ(ierr); 301 ierr = TSSetDM(*ts, dm); CHKERRQ(ierr); 302 if (phys->implicit) { 303 ierr = TSSetType(*ts, TSBDF); CHKERRQ(ierr); 304 if (user->op_ifunction) { 305 ierr = TSSetIFunction(*ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 306 } else { // Implicit integrators can fall back to using an RHSFunction 307 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 308 } 309 } else { 310 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 311 "Problem does not provide RHSFunction"); 312 ierr = TSSetType(*ts, TSRK); CHKERRQ(ierr); 313 ierr = TSRKSetType(*ts, TSRK5F); CHKERRQ(ierr); 314 ierr = TSSetRHSFunction(*ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 315 } 316 ierr = TSSetMaxTime(*ts, 500. * user->units->second); CHKERRQ(ierr); 317 ierr = TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 318 ierr = TSSetTimeStep(*ts, 1.e-2 * user->units->second); CHKERRQ(ierr); 319 if (app_ctx->test_mode) {ierr = TSSetMaxSteps(*ts, 10); CHKERRQ(ierr);} 320 ierr = TSGetAdapt(*ts, &adapt); CHKERRQ(ierr); 321 ierr = TSAdaptSetStepLimits(adapt, 322 1.e-12 * user->units->second, 323 1.e2 * user->units->second); CHKERRQ(ierr); 324 ierr = TSSetFromOptions(*ts); CHKERRQ(ierr); 325 if (!app_ctx->cont_steps) { // print initial condition 326 if (!app_ctx->test_mode) { 327 ierr = TSMonitor_NS(*ts, 0, 0., *Q, user); CHKERRQ(ierr); 328 } 329 } else { // continue from time of last output 330 PetscReal time; 331 PetscInt count; 332 PetscViewer viewer; 333 char file_path[PETSC_MAX_PATH_LEN]; 334 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-time.bin", 335 app_ctx->output_dir); CHKERRQ(ierr); 336 ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer); 337 CHKERRQ(ierr); 338 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 339 CHKERRQ(ierr); 340 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 341 ierr = TSSetTime(*ts, time * user->units->second); CHKERRQ(ierr); 342 } 343 if (!app_ctx->test_mode) { 344 ierr = TSMonitorSet(*ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 345 } 346 347 // Solve 348 double start, cpu_time_used; 349 start = MPI_Wtime(); 350 ierr = PetscBarrier((PetscObject) *ts); CHKERRQ(ierr); 351 ierr = TSSolve(*ts, *Q); CHKERRQ(ierr); 352 cpu_time_used = MPI_Wtime() - start; 353 ierr = TSGetSolveTime(*ts, &final_time); CHKERRQ(ierr); 354 *f_time = final_time; 355 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 356 comm); CHKERRQ(ierr); 357 if (!app_ctx->test_mode) { 358 ierr = PetscPrintf(PETSC_COMM_WORLD, 359 "Time taken for solution (sec): %g\n", 360 (double)cpu_time_used); CHKERRQ(ierr); 361 } 362 PetscFunctionReturn(0); 363 } 364