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 // libCEED + PETSc Example: Navier-Stokes 18 // 19 // This example demonstrates a simple usage of libCEED with PETSc to solve a 20 // Navier-Stokes problem. 21 // 22 // The code is intentionally "raw", using only low-level communication 23 // primitives. 24 // 25 // Build with: 26 // 27 // make [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] navierstokes 28 // 29 // Sample runs: 30 // 31 // ./navierstokes -ceed /cpu/self -problem density_current -degree_volume 1 32 // ./navierstokes -ceed /gpu/occa -problem advection -degree_volume 1 33 // 34 //TESTARGS -ceed {ceed_resource} -test -degree_volume 1 35 36 /// @file 37 /// Navier-Stokes example using PETSc 38 39 const char help[] = "Solve Navier-Stokes using PETSc and libCEED\n"; 40 41 #include <petscts.h> 42 #include <petscdmplex.h> 43 #include <ceed.h> 44 #include <stdbool.h> 45 #include <petscsys.h> 46 #include "common.h" 47 #include "advection.h" 48 #include "advection2d.h" 49 #include "densitycurrent.h" 50 51 // Problem Options 52 typedef enum { 53 NS_DENSITY_CURRENT = 0, 54 NS_ADVECTION = 1, 55 NS_ADVECTION2D = 2, 56 } problemType; 57 static const char *const problemTypes[] = { 58 "density_current", 59 "advection", 60 "advection2d", 61 "problemType","NS_",0 62 }; 63 64 typedef enum { 65 STAB_NONE = 0, 66 STAB_SU = 1, // Streamline Upwind 67 STAB_SUPG = 2, // Streamline Upwind Petrov-Galerkin 68 } StabilizationType; 69 static const char *const StabilizationTypes[] = { 70 "NONE", 71 "SU", 72 "SUPG", 73 "StabilizationType", "STAB_", NULL 74 }; 75 76 // Problem specific data 77 typedef struct { 78 CeedInt dim, qdatasizeVol, qdatasizeSur; 79 CeedQFunctionUser setupVol, setupSur, ics, applyVol_rhs, applySur_rhs, 80 applyVol_ifunction, applySur_ifunction; 81 PetscErrorCode (*bc)(PetscInt, PetscReal, const PetscReal[], PetscInt, 82 PetscScalar[], void *); 83 const char *setupVol_loc, *setupSur_loc, *ics_loc, *applyVol_rhs_loc, 84 *applySur_rhs_loc, *applyVol_ifunction_loc, *applySur_ifunction_loc; 85 const bool non_zero_time; 86 } problemData; 87 88 problemData problemOptions[] = { 89 [NS_DENSITY_CURRENT] = { 90 .dim = 3, 91 .qdatasizeVol = 10, 92 //.qdatasizeSur = 5, 93 .setupVol = Setup_Vol, 94 .setupVol_loc = Setup_Vol_loc, 95 //.setupSur = Setup_Sur, 96 //.setupSur_loc = Setup_Sur_loc, 97 .ics = ICsDC, 98 .ics_loc = ICsDC_loc, 99 .applyVol_rhs = DC_Vol, 100 .applyVol_rhs_loc = DC_Vol_loc, 101 //.applySur_rhs = DC_Sur, 102 //.applySur_rhs_loc = DC_Sur_loc, 103 .applyVol_ifunction = IFunction_DC_Vol, 104 .applyVol_ifunction_loc = IFunction_DC_Vol_loc, 105 //.applySur_ifunction = IFunction_DC_Sur, 106 //.applySur_ifunction_loc = IFunction_DC_Sur_loc, 107 .bc = Exact_DC, 108 .non_zero_time = false, 109 }, 110 [NS_ADVECTION] = { 111 .dim = 3, 112 .qdatasizeVol = 10, 113 //.qdatasizeSur = 5, 114 .setupVol = Setup_Vol, 115 .setupVol_loc = Setup_Vol_loc, 116 //.setupSur = Setup_Sur, 117 //.setupSur_loc = Setup_Sur_loc, 118 .ics = ICsAdvection, 119 .ics_loc = ICsAdvection_loc, 120 .applyVol_rhs = Advection_Vol, 121 .applyVol_rhs_loc = Advection_Vol_loc, 122 //.applySur_rhs = Advection_Sur, 123 //.applySur_rhs_loc = Advection_Sur_loc, 124 .applyVol_ifunction = IFunction_Advection_Vol, 125 .applyVol_ifunction_loc = IFunction_Advection_Vol_loc, 126 //.applySur_ifunction = IFunction_Advection_Sur, 127 //.applySur_ifunction_loc = IFunction_Advection_Sur_loc, 128 .bc = Exact_Advection, 129 .non_zero_time = true, 130 }, 131 [NS_ADVECTION2D] = { 132 .dim = 2, 133 .qdatasizeVol = 5, 134 //.qdatasizeSur = 2, 135 .setupVol = Setup2d_Vol, 136 .setupVol_loc = Setup2d_Vol_loc, 137 //.setupSur = Setup2d_Sur, 138 //.setupSur_loc = Setup2d_Sur_loc, 139 .ics = ICsAdvection2d, 140 .ics_loc = ICsAdvection2d_loc, 141 .applyVol_rhs = Advection2d_Vol, 142 .applyVol_rhs_loc = Advection2d_Vol_loc, 143 //.applySur_rhs = Advection2d_Sur, 144 //.applySur_rhs_loc = Advection2d_Sur_loc, 145 .applyVol_ifunction = IFunction_Advection2d_Vol, 146 .applyVol_ifunction_loc = IFunction_Advection2d_Vol_loc, 147 //.applySur_ifunction = IFunction_Advection2d_Sur, 148 //.applySur_ifunction_loc = IFunction_Advection2d_Sur_loc, 149 .bc = Exact_Advection2d, 150 .non_zero_time = true, 151 }, 152 }; 153 154 // PETSc user data 155 typedef struct User_ *User; 156 typedef struct Units_ *Units; 157 158 struct User_ { 159 MPI_Comm comm; 160 PetscInt outputfreq; 161 DM dm; 162 DM dmviz; 163 Mat interpviz; 164 Ceed ceed; 165 Units units; 166 CeedVector qceed, qdotceed, gceed; 167 CeedOperator op_rhs_vol, op_rhs_sur, op_rhs, 168 op_ifunction_vol, op_ifunction_sur, op_ifunction; 169 Vec M; 170 char outputfolder[PETSC_MAX_PATH_LEN]; 171 PetscInt contsteps; 172 }; 173 174 struct Units_ { 175 // fundamental units 176 PetscScalar meter; 177 PetscScalar kilogram; 178 PetscScalar second; 179 PetscScalar Kelvin; 180 // derived units 181 PetscScalar Pascal; 182 PetscScalar JperkgK; 183 PetscScalar mpersquareds; 184 PetscScalar WpermK; 185 PetscScalar kgpercubicm; 186 PetscScalar kgpersquaredms; 187 PetscScalar Joulepercubicm; 188 }; 189 190 typedef struct SimpleBC_ *SimpleBC; 191 struct SimpleBC_ { 192 PetscInt nwall, nslip[3]; 193 PetscInt walls[10], slips[3][10]; 194 }; 195 196 // Essential BC dofs are encoded in closure indices as -(i+1). 197 static PetscInt Involute(PetscInt i) { 198 return i >= 0 ? i : -(i+1); 199 } 200 201 // Utility function to create local CEED restriction 202 static PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 203 CeedInt height, DMLabel domainLabel, CeedInt value, 204 CeedElemRestriction *Erestrict) { 205 206 PetscSection section; 207 PetscInt p, Nelem, Ndof, *erestrict, eoffset, nfields, dim, 208 depth; 209 DMLabel depthLabel; 210 IS depthIS, iterIS; 211 const PetscInt *iterIndices; 212 PetscErrorCode ierr; 213 Vec Uloc; 214 215 PetscFunctionBeginUser; 216 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 217 ierr = DMGetLocalSection(dm,§ion); CHKERRQ(ierr); 218 ierr = PetscSectionGetNumFields(section, &nfields); CHKERRQ(ierr); 219 PetscInt ncomp[nfields], fieldoff[nfields+1]; 220 fieldoff[0] = 0; 221 for (PetscInt f=0; f<nfields; f++) { 222 ierr = PetscSectionGetFieldComponents(section, f, &ncomp[f]); CHKERRQ(ierr); 223 fieldoff[f+1] = fieldoff[f] + ncomp[f]; 224 } 225 226 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 227 ierr = DMPlexGetDepthLabel(dm, &depthLabel); CHKERRQ(ierr); 228 ierr = DMLabelGetStratumIS(depthLabel, depth - height, &depthIS); CHKERRQ(ierr); 229 if (domainLabel) { 230 IS domainIS; 231 ierr = DMLabelGetStratumIS(domainLabel, value, &domainIS); CHKERRQ(ierr); 232 ierr = ISIntersect(depthIS, domainIS, &iterIS); CHKERRQ(ierr); 233 ierr = ISDestroy(&domainIS); CHKERRQ(ierr); 234 ierr = ISDestroy(&depthIS); CHKERRQ(ierr); 235 } else { 236 iterIS = depthIS; 237 } 238 ierr = ISGetLocalSize(iterIS, &Nelem); CHKERRQ(ierr); 239 ierr = ISGetIndices(iterIS, &iterIndices); CHKERRQ(ierr); 240 ierr = PetscMalloc1(Nelem*PetscPowInt(P, dim), &erestrict); CHKERRQ(ierr); 241 for (p=0,eoffset=0; p<Nelem; p++) { 242 PetscInt c = iterIndices[p]; 243 PetscInt numindices, *indices, nnodes; 244 ierr = DMPlexGetClosureIndices(dm, section, section, c, &numindices, 245 &indices, NULL); CHKERRQ(ierr); 246 if (numindices % fieldoff[nfields]) 247 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 248 "Number of closure indices not compatible with Cell %D", c); 249 nnodes = numindices / fieldoff[nfields]; 250 for (PetscInt i=0; i<nnodes; i++) { 251 // Check that indices are blocked by node and thus can be coalesced as a single field with 252 // fieldoff[nfields] = sum(ncomp) components. 253 for (PetscInt f=0; f<nfields; f++) { 254 for (PetscInt j=0; j<ncomp[f]; j++) { 255 if (Involute(indices[fieldoff[f]*nnodes + i*ncomp[f] + j]) 256 != Involute(indices[i*ncomp[0]]) + fieldoff[f] + j) 257 SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, 258 "Cell %D closure indices not interlaced for node %D field %D component %D", 259 c, i, f, j); 260 } 261 } 262 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 263 PetscInt loc = Involute(indices[i*ncomp[0]]); 264 erestrict[eoffset++] = loc; 265 } 266 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, &numindices, 267 &indices, NULL); CHKERRQ(ierr); 268 } 269 if (eoffset != Nelem*PetscPowInt(P, dim)) 270 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 271 "ElemRestriction of size (%D,%D) initialized %D nodes", Nelem, 272 PetscPowInt(P, dim),eoffset); 273 ierr = ISRestoreIndices(iterIS, &iterIndices); CHKERRQ(ierr); 274 ierr = ISDestroy(&iterIS); CHKERRQ(ierr); 275 276 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 277 ierr = VecGetLocalSize(Uloc, &Ndof); CHKERRQ(ierr); 278 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 279 CeedElemRestrictionCreate(ceed, Nelem, PetscPowInt(P, dim), fieldoff[nfields], 280 1, Ndof, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 281 Erestrict); 282 ierr = PetscFree(erestrict); CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286 static int CreateVectorFromPetscVec(Ceed ceed, Vec p, CeedVector *v) { 287 PetscErrorCode ierr; 288 PetscInt m; 289 290 PetscFunctionBeginUser; 291 ierr = VecGetLocalSize(p, &m); CHKERRQ(ierr); 292 ierr = CeedVectorCreate(ceed, m, v); CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 static int VectorPlacePetscVec(CeedVector c, Vec p) { 297 PetscErrorCode ierr; 298 PetscInt mceed,mpetsc; 299 PetscScalar *a; 300 301 PetscFunctionBeginUser; 302 ierr = CeedVectorGetLength(c, &mceed); CHKERRQ(ierr); 303 ierr = VecGetLocalSize(p, &mpetsc); CHKERRQ(ierr); 304 if (mceed != mpetsc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, 305 "Cannot place PETSc Vec of length %D in CeedVector of length %D",mpetsc,mceed); 306 ierr = VecGetArray(p, &a); CHKERRQ(ierr); 307 CeedVectorSetArray(c, CEED_MEM_HOST, CEED_USE_POINTER, a); 308 PetscFunctionReturn(0); 309 } 310 311 static PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 312 PetscBool insertEssential, Vec Qloc, PetscReal time, Vec faceGeomFVM, 313 Vec cellGeomFVM, Vec gradFVM) { 314 PetscErrorCode ierr; 315 Vec Qbc; 316 317 PetscFunctionBegin; 318 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 319 ierr = VecAXPY(Qloc, 1., Qbc); CHKERRQ(ierr); 320 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 321 PetscFunctionReturn(0); 322 } 323 324 // This is the RHS of the ODE, given as u_t = G(t,u) 325 // This function takes in a state vector Q and writes into G 326 static PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *userData) { 327 PetscErrorCode ierr; 328 User user = *(User *)userData; 329 PetscScalar *q, *g; 330 Vec Qloc, Gloc; 331 332 // Global-to-local 333 PetscFunctionBeginUser; 334 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 335 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 336 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 337 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 338 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 339 NULL, NULL, NULL); CHKERRQ(ierr); 340 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 341 342 // Ceed Vectors 343 ierr = VecGetArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 344 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 345 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); 346 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 347 348 // Apply CEED operator 349 CeedOperatorApply(user->op_rhs, user->qceed, user->gceed, 350 CEED_REQUEST_IMMEDIATE); 351 352 // Restore vectors 353 ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 354 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 355 356 ierr = VecZeroEntries(G); CHKERRQ(ierr); 357 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 358 359 // Inverse of the lumped mass matrix 360 ierr = VecPointwiseMult(G, G, user->M); // M is Minv 361 CHKERRQ(ierr); 362 363 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 364 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367 368 static PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Qdot, Vec G, 369 void *userData) { 370 PetscErrorCode ierr; 371 User user = *(User *)userData; 372 const PetscScalar *q, *qdot; 373 PetscScalar *g; 374 Vec Qloc, Qdotloc, Gloc; 375 376 // Global-to-local 377 PetscFunctionBeginUser; 378 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 379 ierr = DMGetLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 380 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 381 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 382 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 383 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 384 NULL, NULL, NULL); CHKERRQ(ierr); 385 ierr = VecZeroEntries(Qdotloc); CHKERRQ(ierr); 386 ierr = DMGlobalToLocal(user->dm, Qdot, INSERT_VALUES, Qdotloc); CHKERRQ(ierr); 387 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 388 389 // Ceed Vectors 390 ierr = VecGetArrayRead(Qloc, &q); CHKERRQ(ierr); 391 ierr = VecGetArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 392 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 393 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, 394 (PetscScalar *)q); 395 CeedVectorSetArray(user->qdotceed, CEED_MEM_HOST, CEED_USE_POINTER, 396 (PetscScalar *)qdot); 397 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 398 399 // Apply CEED operator 400 CeedOperatorApply(user->op_ifunction, user->qceed, user->gceed, 401 CEED_REQUEST_IMMEDIATE); 402 403 // Restore vectors 404 ierr = VecRestoreArrayRead(Qloc, &q); CHKERRQ(ierr); 405 ierr = VecRestoreArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 406 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 407 408 ierr = VecZeroEntries(G); CHKERRQ(ierr); 409 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 410 411 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 412 ierr = DMRestoreLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 413 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 414 PetscFunctionReturn(0); 415 } 416 417 // User provided TS Monitor 418 static PetscErrorCode TSMonitor_NS(TS ts, PetscInt stepno, PetscReal time, 419 Vec Q, void *ctx) { 420 User user = ctx; 421 Vec Qloc; 422 char filepath[PETSC_MAX_PATH_LEN]; 423 PetscViewer viewer; 424 PetscErrorCode ierr; 425 426 // Set up output 427 PetscFunctionBeginUser; 428 // Print every 'outputfreq' steps 429 if (stepno % user->outputfreq != 0) 430 PetscFunctionReturn(0); 431 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 432 ierr = PetscObjectSetName((PetscObject)Qloc, "StateVec"); CHKERRQ(ierr); 433 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 434 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 435 436 // Output 437 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-%03D.vtu", 438 user->outputfolder, stepno + user->contsteps); 439 CHKERRQ(ierr); 440 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), filepath, 441 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 442 ierr = VecView(Qloc, viewer); CHKERRQ(ierr); 443 if (user->dmviz) { 444 Vec Qrefined, Qrefined_loc; 445 char filepath_refined[PETSC_MAX_PATH_LEN]; 446 PetscViewer viewer_refined; 447 448 ierr = DMGetGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 449 ierr = DMGetLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 450 ierr = PetscObjectSetName((PetscObject)Qrefined_loc, "Refined"); 451 CHKERRQ(ierr); 452 ierr = MatInterpolate(user->interpviz, Q, Qrefined); CHKERRQ(ierr); 453 ierr = VecZeroEntries(Qrefined_loc); CHKERRQ(ierr); 454 ierr = DMGlobalToLocal(user->dmviz, Qrefined, INSERT_VALUES, Qrefined_loc); 455 CHKERRQ(ierr); 456 ierr = PetscSNPrintf(filepath_refined, sizeof filepath_refined, 457 "%s/nsrefined-%03D.vtu", 458 user->outputfolder, stepno + user->contsteps); 459 CHKERRQ(ierr); 460 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Qrefined), 461 filepath_refined, 462 FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 463 ierr = VecView(Qrefined_loc, viewer_refined); CHKERRQ(ierr); 464 ierr = DMRestoreLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 465 ierr = DMRestoreGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 466 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 467 } 468 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 469 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 470 471 // Save data in a binary file for continuation of simulations 472 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 473 user->outputfolder); CHKERRQ(ierr); 474 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 475 CHKERRQ(ierr); 476 ierr = VecView(Q, viewer); CHKERRQ(ierr); 477 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 478 479 // Save time stamp 480 // Dimensionalize time back 481 time /= user->units->second; 482 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 483 user->outputfolder); CHKERRQ(ierr); 484 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 485 CHKERRQ(ierr); 486 #if PETSC_VERSION_GE(3,13,0) 487 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 488 #else 489 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 490 #endif 491 CHKERRQ(ierr); 492 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 493 494 PetscFunctionReturn(0); 495 } 496 497 static PetscErrorCode ICs_PetscMultiplicity(CeedOperator op_ics, 498 CeedVector xcorners, CeedVector q0ceed, DM dm, Vec Qloc, Vec Q, 499 CeedElemRestriction restrictq, SetupContext ctxSetup, CeedScalar time) { 500 PetscErrorCode ierr; 501 CeedVector multlvec; 502 Vec Multiplicity, MultiplicityLoc; 503 504 ctxSetup->time = time; 505 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 506 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 507 CeedOperatorApply(op_ics, xcorners, q0ceed, CEED_REQUEST_IMMEDIATE); 508 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 509 ierr = DMLocalToGlobal(dm, Qloc, ADD_VALUES, Q); CHKERRQ(ierr); 510 511 // Fix multiplicity for output of ICs 512 ierr = DMGetLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 513 CeedElemRestrictionCreateVector(restrictq, &multlvec, NULL); 514 ierr = VectorPlacePetscVec(multlvec, MultiplicityLoc); CHKERRQ(ierr); 515 CeedElemRestrictionGetMultiplicity(restrictq, multlvec); 516 CeedVectorDestroy(&multlvec); 517 ierr = DMGetGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 518 ierr = VecZeroEntries(Multiplicity); CHKERRQ(ierr); 519 ierr = DMLocalToGlobal(dm, MultiplicityLoc, ADD_VALUES, Multiplicity); 520 CHKERRQ(ierr); 521 ierr = VecPointwiseDivide(Q, Q, Multiplicity); CHKERRQ(ierr); 522 ierr = VecPointwiseDivide(Qloc, Qloc, MultiplicityLoc); CHKERRQ(ierr); 523 ierr = DMRestoreLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 524 ierr = DMRestoreGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 525 526 PetscFunctionReturn(0); 527 } 528 529 static PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, 530 CeedElemRestriction restrictq, CeedBasis basisq, 531 CeedElemRestriction restrictqdi, CeedVector qdata, Vec M) { 532 PetscErrorCode ierr; 533 CeedQFunction qf_mass; 534 CeedOperator op_mass; 535 CeedVector mceed; 536 Vec Mloc; 537 CeedInt ncompq, qdatasize; 538 539 PetscFunctionBeginUser; 540 CeedElemRestrictionGetNumComponents(restrictq, &ncompq); 541 CeedElemRestrictionGetNumComponents(restrictqdi, &qdatasize); 542 // Create the Q-function that defines the action of the mass operator 543 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 544 CeedQFunctionAddInput(qf_mass, "q", ncompq, CEED_EVAL_INTERP); 545 CeedQFunctionAddInput(qf_mass, "qdata", qdatasize, CEED_EVAL_NONE); 546 CeedQFunctionAddOutput(qf_mass, "v", ncompq, CEED_EVAL_INTERP); 547 548 // Create the mass operator 549 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 550 CeedOperatorSetField(op_mass, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 551 CeedOperatorSetField(op_mass, "qdata", restrictqdi, 552 CEED_BASIS_COLLOCATED, qdata); 553 CeedOperatorSetField(op_mass, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 554 555 ierr = DMGetLocalVector(dm, &Mloc); CHKERRQ(ierr); 556 ierr = VecZeroEntries(Mloc); CHKERRQ(ierr); 557 CeedElemRestrictionCreateVector(restrictq, &mceed, NULL); 558 ierr = VectorPlacePetscVec(mceed, Mloc); CHKERRQ(ierr); 559 560 { 561 // Compute a lumped mass matrix 562 CeedVector onesvec; 563 CeedElemRestrictionCreateVector(restrictq, &onesvec, NULL); 564 CeedVectorSetValue(onesvec, 1.0); 565 CeedOperatorApply(op_mass, onesvec, mceed, CEED_REQUEST_IMMEDIATE); 566 CeedVectorDestroy(&onesvec); 567 CeedOperatorDestroy(&op_mass); 568 CeedVectorDestroy(&mceed); 569 } 570 CeedQFunctionDestroy(&qf_mass); 571 572 ierr = VecZeroEntries(M); CHKERRQ(ierr); 573 ierr = DMLocalToGlobal(dm, Mloc, ADD_VALUES, M); CHKERRQ(ierr); 574 ierr = DMRestoreLocalVector(dm, &Mloc); CHKERRQ(ierr); 575 576 // Invert diagonally lumped mass vector for RHS function 577 ierr = VecReciprocal(M); CHKERRQ(ierr); 578 PetscFunctionReturn(0); 579 } 580 581 PetscErrorCode SetUpDM(DM dm, problemData *problem, PetscInt degree, 582 SimpleBC bc, void *ctxSetup) { 583 PetscErrorCode ierr; 584 585 PetscFunctionBeginUser; 586 { 587 // Configure the finite element space and boundary conditions 588 PetscFE fe; 589 PetscSpace fespace; 590 PetscInt ncompq = 5; 591 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, ncompq, 592 PETSC_FALSE, degree, PETSC_DECIDE, 593 &fe); CHKERRQ(ierr); 594 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 595 ierr = DMAddField(dm,NULL,(PetscObject)fe); CHKERRQ(ierr); 596 ierr = DMCreateDS(dm); CHKERRQ(ierr); 597 /* Wall boundary conditions are zero velocity and zero flux for density and energy */ 598 { 599 PetscInt comps[3] = {1, 2, 3}; 600 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 601 3, comps, (void(*)(void))problem->bc, 602 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 603 } 604 { 605 PetscInt comps[1] = {1}; 606 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 607 1, comps, (void(*)(void))NULL, bc->nslip[0], 608 bc->slips[0], ctxSetup); CHKERRQ(ierr); 609 comps[0] = 2; 610 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 611 1, comps, (void(*)(void))NULL, bc->nslip[1], 612 bc->slips[1], ctxSetup); CHKERRQ(ierr); 613 comps[0] = 3; 614 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 615 1, comps, (void(*)(void))NULL, bc->nslip[2], 616 bc->slips[2], ctxSetup); CHKERRQ(ierr); 617 } 618 ierr = DMPlexSetClosurePermutationTensor(dm,PETSC_DETERMINE,NULL); 619 CHKERRQ(ierr); 620 ierr = PetscFEGetBasisSpace(fe, &fespace); CHKERRQ(ierr); 621 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 622 } 623 { 624 // Empty name for conserved field (because there is only one field) 625 PetscSection section; 626 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 627 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 628 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 629 CHKERRQ(ierr); 630 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 631 CHKERRQ(ierr); 632 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 633 CHKERRQ(ierr); 634 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 635 CHKERRQ(ierr); 636 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 637 CHKERRQ(ierr); 638 } 639 PetscFunctionReturn(0); 640 } 641 642 int main(int argc, char **argv) { 643 PetscInt ierr; 644 MPI_Comm comm; 645 DM dm, dmcoord, dmviz, dmvizfine; 646 Mat interpviz; 647 TS ts; 648 TSAdapt adapt; 649 User user; 650 Units units; 651 char ceedresource[4096] = "/cpu/self"; 652 PetscInt cStart, cEnd, localNelem, lnodes, steps; 653 const PetscInt ncompq = 5; 654 PetscMPIInt rank; 655 PetscScalar ftime; 656 Vec Q, Qloc, Xloc; 657 Ceed ceed; 658 CeedInt numP_Vol, numP_Sur, numQ_Vol, numQ_Sur; 659 CeedVector xcorners, qdataVol, qdataSur, q0ceed; 660 CeedBasis basisxVol, basisxcVol, basisqVol, basisxSur, basisxcSur, basisqSur; 661 CeedElemRestriction restrictxVol, restrictxcoordVol, restrictqVol, restrictqdiVol, 662 restrictxSur, restrictxcoordSur, restrictqSur, restrictqdiSur; 663 CeedQFunction qf_setupVol, qf_setupSur, qf_ics, qf_rhsVol, qf_rhsSur, 664 qf_ifunctionVol, qf_ifunctionSur; 665 CeedOperator op_setupVol, op_setupSur, op_ics; 666 CeedScalar Rd; 667 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 668 kgpersquaredms, Joulepercubicm; 669 problemType problemChoice; 670 problemData *problem = NULL; 671 StabilizationType stab; 672 PetscBool test, implicit; 673 PetscInt viz_refine = 0; 674 struct SimpleBC_ bc = { 675 .nwall = 6, 676 .walls = {1,2,3,4,5,6}, 677 }; 678 double start, cpu_time_used; 679 680 // Create the libCEED contexts 681 PetscScalar meter = 1e-2; // 1 meter in scaled length units 682 PetscScalar second = 1e-2; // 1 second in scaled time units 683 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 684 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 685 CeedScalar theta0 = 300.; // K 686 CeedScalar thetaC = -15.; // K 687 CeedScalar P0 = 1.e5; // Pa 688 CeedScalar N = 0.01; // 1/s 689 CeedScalar cv = 717.; // J/(kg K) 690 CeedScalar cp = 1004.; // J/(kg K) 691 CeedScalar g = 9.81; // m/s^2 692 CeedScalar lambda = -2./3.; // - 693 CeedScalar mu = 75.; // Pa s, dynamic viscosity 694 // mu = 75 is not physical for air, but is good for numerical stability 695 CeedScalar k = 0.02638; // W/(m K) 696 CeedScalar CtauS = 0.; // dimensionless 697 CeedScalar strong_form = 0.; // [0,1] 698 PetscScalar lx = 8000.; // m 699 PetscScalar ly = 8000.; // m 700 PetscScalar lz = 4000.; // m 701 CeedScalar rc = 1000.; // m (Radius of bubble) 702 PetscScalar resx = 1000.; // m (resolution in x) 703 PetscScalar resy = 1000.; // m (resolution in y) 704 PetscScalar resz = 1000.; // m (resolution in z) 705 PetscInt outputfreq = 10; // - 706 PetscInt contsteps = 0; // - 707 PetscInt degreeVol = 1; // - 708 PetscInt degreeSur = 1; // - 709 PetscInt qextraVol = 2; // - 710 PetscInt qextraSur = 2; // - 711 DMBoundaryType periodicity[] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 712 DM_BOUNDARY_NONE 713 }; 714 PetscReal center[3], dc_axis[3] = {0, 0, 0}; 715 716 ierr = PetscInitialize(&argc, &argv, NULL, help); 717 if (ierr) return ierr; 718 719 // Allocate PETSc context 720 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 721 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 722 723 // Parse command line options 724 comm = PETSC_COMM_WORLD; 725 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 726 NULL); CHKERRQ(ierr); 727 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 728 NULL, ceedresource, ceedresource, 729 sizeof(ceedresource), NULL); CHKERRQ(ierr); 730 ierr = PetscOptionsBool("-test", "Run in test mode", 731 NULL, test=PETSC_FALSE, &test, NULL); CHKERRQ(ierr); 732 problemChoice = NS_DENSITY_CURRENT; 733 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 734 problemTypes, (PetscEnum)problemChoice, 735 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 736 problem = &problemOptions[problemChoice]; 737 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 738 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 739 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 740 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 741 NULL, implicit=PETSC_FALSE, &implicit, NULL); 742 CHKERRQ(ierr); 743 { 744 PetscInt len; 745 PetscBool flg; 746 ierr = PetscOptionsIntArray("-bc_wall", 747 "Use wall boundary conditions on this list of faces", 748 NULL, bc.walls, 749 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 750 &len), &flg); CHKERRQ(ierr); 751 if (flg) bc.nwall = len; 752 for (PetscInt j=0; j<3; j++) { 753 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 754 ierr = PetscOptionsIntArray(flags[j], 755 "Use slip boundary conditions on this list of faces", 756 NULL, bc.slips[j], 757 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 758 &len), &flg); 759 CHKERRQ(ierr); 760 if (flg) bc.nslip[j] = len; 761 } 762 } 763 ierr = PetscOptionsInt("-viz_refine", 764 "Regular refinement levels for visualization", 765 NULL, viz_refine, &viz_refine, NULL); 766 CHKERRQ(ierr); 767 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 768 NULL, meter, &meter, NULL); CHKERRQ(ierr); 769 meter = fabs(meter); 770 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 771 NULL, second, &second, NULL); CHKERRQ(ierr); 772 second = fabs(second); 773 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 774 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 775 kilogram = fabs(kilogram); 776 ierr = PetscOptionsScalar("-units_Kelvin", 777 "1 Kelvin in scaled temperature units", 778 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 779 Kelvin = fabs(Kelvin); 780 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 781 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 782 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 783 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 784 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 785 NULL, P0, &P0, NULL); CHKERRQ(ierr); 786 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 787 NULL, N, &N, NULL); CHKERRQ(ierr); 788 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 789 NULL, cv, &cv, NULL); CHKERRQ(ierr); 790 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 791 NULL, cp, &cp, NULL); CHKERRQ(ierr); 792 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 793 NULL, g, &g, NULL); CHKERRQ(ierr); 794 ierr = PetscOptionsScalar("-lambda", 795 "Stokes hypothesis second viscosity coefficient", 796 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 797 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 798 NULL, mu, &mu, NULL); CHKERRQ(ierr); 799 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 800 NULL, k, &k, NULL); CHKERRQ(ierr); 801 ierr = PetscOptionsScalar("-CtauS", 802 "Scale coefficient for tau (nondimensional)", 803 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 804 ierr = PetscOptionsScalar("-strong_form", 805 "Strong (1) or weak/integrated by parts (0) advection residual", 806 NULL, strong_form, &strong_form, NULL); 807 CHKERRQ(ierr); 808 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 809 NULL, lx, &lx, NULL); CHKERRQ(ierr); 810 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 811 NULL, ly, &ly, NULL); CHKERRQ(ierr); 812 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 813 NULL, lz, &lz, NULL); CHKERRQ(ierr); 814 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 815 NULL, rc, &rc, NULL); CHKERRQ(ierr); 816 ierr = PetscOptionsScalar("-resx","Target resolution in x", 817 NULL, resx, &resx, NULL); CHKERRQ(ierr); 818 ierr = PetscOptionsScalar("-resy","Target resolution in y", 819 NULL, resy, &resy, NULL); CHKERRQ(ierr); 820 ierr = PetscOptionsScalar("-resz","Target resolution in z", 821 NULL, resz, &resz, NULL); CHKERRQ(ierr); 822 PetscInt n = problem->dim; 823 ierr = PetscOptionsEnumArray("-periodicity", "Periodicity per direction", 824 NULL, DMBoundaryTypes, (PetscEnum *)periodicity, 825 &n, NULL); CHKERRQ(ierr); 826 n = problem->dim; 827 center[0] = 0.5 * lx; 828 center[1] = 0.5 * ly; 829 center[2] = 0.5 * lz; 830 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 831 NULL, center, &n, NULL); CHKERRQ(ierr); 832 n = problem->dim; 833 ierr = PetscOptionsRealArray("-dc_axis", 834 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 835 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 836 { 837 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 838 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 839 if (norm > 0) { 840 for (int i=0; i<3; i++) dc_axis[i] /= norm; 841 } 842 } 843 ierr = PetscOptionsInt("-output_freq", 844 "Frequency of output, in number of steps", 845 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 846 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 847 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 848 ierr = PetscOptionsInt("-degree_volume", "Polynomial degree of finite elements for the Volume", 849 NULL, degreeVol, °reeVol, NULL); CHKERRQ(ierr); 850 ierr = PetscOptionsInt("-qextra_volume", "Number of extra quadrature points for the Volume", 851 NULL, qextraVol, &qextraVol, NULL); CHKERRQ(ierr); 852 PetscStrncpy(user->outputfolder, ".", 2); 853 ierr = PetscOptionsString("-of", "Output folder", 854 NULL, user->outputfolder, user->outputfolder, 855 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 856 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 857 858 // Define derived units 859 Pascal = kilogram / (meter * PetscSqr(second)); 860 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 861 mpersquareds = meter / PetscSqr(second); 862 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 863 kgpercubicm = kilogram / pow(meter,3); 864 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 865 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 866 867 // Scale variables to desired units 868 theta0 *= Kelvin; 869 thetaC *= Kelvin; 870 P0 *= Pascal; 871 N *= (1./second); 872 cv *= JperkgK; 873 cp *= JperkgK; 874 Rd = cp - cv; 875 g *= mpersquareds; 876 mu *= Pascal * second; 877 k *= WpermK; 878 lx = fabs(lx) * meter; 879 ly = fabs(ly) * meter; 880 lz = fabs(lz) * meter; 881 rc = fabs(rc) * meter; 882 resx = fabs(resx) * meter; 883 resy = fabs(resy) * meter; 884 resz = fabs(resz) * meter; 885 for (int i=0; i<3; i++) center[i] *= meter; 886 887 const CeedInt dim = problem->dim, ncompx = problem->dim, 888 qdatasizeVol = problem->qdatasizeVol; 889 // Set up the libCEED context 890 struct SetupContext_ ctxSetup = { 891 .theta0 = theta0, 892 .thetaC = thetaC, 893 .P0 = P0, 894 .N = N, 895 .cv = cv, 896 .cp = cp, 897 .Rd = Rd, 898 .g = g, 899 .rc = rc, 900 .lx = lx, 901 .ly = ly, 902 .lz = lz, 903 .periodicity0 = periodicity[0], 904 .periodicity1 = periodicity[1], 905 .periodicity2 = periodicity[2], 906 .center[0] = center[0], 907 .center[1] = center[1], 908 .center[2] = center[2], 909 .dc_axis[0] = dc_axis[0], 910 .dc_axis[1] = dc_axis[1], 911 .dc_axis[2] = dc_axis[2], 912 .time = 0, 913 }; 914 915 { 916 const PetscReal scale[3] = {lx, ly, lz}; 917 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 918 periodicity, PETSC_TRUE, &dm); 919 CHKERRQ(ierr); 920 } 921 if (1) { 922 DM dmDist = NULL; 923 PetscPartitioner part; 924 925 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 926 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 927 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 928 if (dmDist) { 929 ierr = DMDestroy(&dm); CHKERRQ(ierr); 930 dm = dmDist; 931 } 932 } 933 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 934 935 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 936 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 937 ierr = SetUpDM(dm, problem, degreeVol, &bc, &ctxSetup); CHKERRQ(ierr); 938 if (!test) { 939 ierr = PetscPrintf(PETSC_COMM_WORLD, 940 "Degree of Volumetric FEM Space: %D\n", 941 degreeVol); CHKERRQ(ierr); 942 } 943 dmviz = NULL; 944 interpviz = NULL; 945 if (viz_refine) { 946 DM dmhierarchy[viz_refine+1]; 947 948 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 949 dmhierarchy[0] = dm; 950 for (PetscInt i = 0, d = degreeVol; i < viz_refine; i++) { 951 Mat interp_next; 952 953 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 954 CHKERRQ(ierr); 955 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 956 d = (d + 1) / 2; 957 if (i + 1 == viz_refine) d = 1; 958 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 959 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 960 &interp_next, NULL); CHKERRQ(ierr); 961 if (!i) interpviz = interp_next; 962 else { 963 Mat C; 964 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 965 PETSC_DECIDE, &C); CHKERRQ(ierr); 966 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 967 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 968 interpviz = C; 969 } 970 } 971 for (PetscInt i=1; i<viz_refine; i++) { 972 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 973 } 974 dmviz = dmhierarchy[viz_refine]; 975 } 976 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 977 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 978 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 979 lnodes /= ncompq; 980 981 { 982 // Print grid information 983 CeedInt gdofs, odofs; 984 int comm_size; 985 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 986 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 987 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 988 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 989 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 990 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 991 if (!test) { 992 ierr = PetscPrintf(comm, "Global FEM dofs: %D (%D owned) on %d ranks\n", 993 gdofs, odofs, comm_size); CHKERRQ(ierr); 994 ierr = PetscPrintf(comm, "Local FEM nodes: %D\n", lnodes); CHKERRQ(ierr); 995 ierr = PetscPrintf(comm, "dm_plex_box_faces: %s\n", box_faces_str); 996 CHKERRQ(ierr); 997 } 998 999 } 1000 1001 // Set up global mass vector 1002 ierr = VecDuplicate(Q,&user->M); CHKERRQ(ierr); 1003 1004 // Set up CEED 1005 // CEED Bases 1006 CeedInit(ceedresource, &ceed); 1007 numP_Vol = degreeVol + 1; 1008 numQ_Vol = numP_Vol + qextraVol; 1009 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP_Vol, numQ_Vol, CEED_GAUSS, 1010 &basisqVol); 1011 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ_Vol, CEED_GAUSS, 1012 &basisxVol); 1013 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP_Vol, 1014 CEED_GAUSS_LOBATTO, &basisxcVol); 1015 1016 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1017 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1018 CHKERRQ(ierr); 1019 1020 // CEED Restrictions 1021 ierr = CreateRestrictionFromPlex(ceed, dm, degreeVol+1, 0, 0, 0, &restrictqVol); 1022 CHKERRQ(ierr); 1023 ierr = CreateRestrictionFromPlex(ceed, dmcoord, 2, 0, 0, 0, &restrictxVol); 1024 CHKERRQ(ierr); 1025 DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr); 1026 localNelem = cEnd - cStart; 1027 CeedInt numQdim = CeedIntPow(numQ_Vol, dim); 1028 CeedElemRestrictionCreateStrided(ceed, localNelem, numQdim, 1029 qdatasize, qdatasize*localNelem*numQdim, 1030 CEED_STRIDES_BACKEND, &restrictqdi); 1031 1032 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1033 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1034 1035 // Create the CEED vectors that will be needed in setup 1036 CeedInt Nqpts; 1037 CeedBasisGetNumQuadraturePoints(basisqVol, &Nqpts); 1038 CeedInt Ndofs = 1; 1039 for (int d=0; d<3; d++) Ndofs *= numP_Vol; 1040 CeedVectorCreate(ceed, qdatasizeVol*localNelem*Nqpts, &qdataVol); 1041 CeedElemRestrictionCreateVector(restrictqVol, &q0ceed, NULL); 1042 1043 // Create the Q-function that builds the quadrature data for the NS operator 1044 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1045 &qf_setupVol); 1046 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1047 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1048 CeedQFunctionAddOutput(qf_setupVol, "qdataVol", qdatasizeVol, CEED_EVAL_NONE); 1049 1050 // Create the Q-function that sets the ICs of the operator 1051 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1052 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1053 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1054 1055 qf_rhsVol = NULL; 1056 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1057 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1058 problem->applyVol_rhs_loc, &qf_rhsVol); 1059 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1060 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1061 CeedQFunctionAddInput(qf_rhsVol, "qdataVol", qdatasizeVol, CEED_EVAL_NONE); 1062 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1063 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1064 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1065 } 1066 1067 qf_ifunctionVol = NULL; 1068 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1069 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1070 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1071 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1072 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1073 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1074 CeedQFunctionAddInput(qf_ifunctionVol, "qdataVol", qdatasizeVol, CEED_EVAL_NONE); 1075 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1076 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1077 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1078 } 1079 1080 // Create the operator that builds the quadrature data for the NS operator 1081 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1082 CeedOperatorSetField(op_setupVol, "dx", restrictxVol, basisxVol, CEED_VECTOR_ACTIVE); 1083 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1084 basisxVol, CEED_VECTOR_NONE); 1085 CeedOperatorSetField(op_setupVol, "qdataVol", restrictqdiVol, 1086 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1087 1088 // Create the operator that sets the ICs 1089 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1090 CeedOperatorSetField(op_ics, "x", restrictxVol, basisxcVol, CEED_VECTOR_ACTIVE); 1091 CeedOperatorSetField(op_ics, "q0", restrictqVol, 1092 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1093 1094 CeedElemRestrictionCreateVector(restrictqVol, &user->qceed, NULL); 1095 CeedElemRestrictionCreateVector(restrictqVol, &user->qdotceed, NULL); 1096 CeedElemRestrictionCreateVector(restrictqVol, &user->gceed, NULL); 1097 1098 if (qf_rhsVol) { // Create the RHS physics operator 1099 CeedOperator op; 1100 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1101 CeedOperatorSetField(op, "q", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1102 CeedOperatorSetField(op, "dq", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1103 CeedOperatorSetField(op, "qdataVol", restrictqdiVol, 1104 CEED_BASIS_COLLOCATED, qdataVol); 1105 CeedOperatorSetField(op, "x", restrictxVol, basisxVol, xcorners); 1106 CeedOperatorSetField(op, "v", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1107 CeedOperatorSetField(op, "dv", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1108 user->op_rhs = op; 1109 } 1110 1111 if (qf_ifunctionVol) { // Create the IFunction operator 1112 CeedOperator op; 1113 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1114 CeedOperatorSetField(op, "q", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1115 CeedOperatorSetField(op, "dq", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1116 CeedOperatorSetField(op, "qdot", restrictqVol, basisqVol, user->qdotceed); 1117 CeedOperatorSetField(op, "qdataVol", restrictqdiVol, 1118 CEED_BASIS_COLLOCATED, qdataVol); 1119 CeedOperatorSetField(op, "x", restrictxVol, basisxVol, xcorners); 1120 CeedOperatorSetField(op, "v", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1121 CeedOperatorSetField(op, "dv", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1122 user->op_ifunction = op; 1123 } 1124 1125 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1126 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1127 struct Advection2dContext_ ctxAdvection2d = { 1128 .CtauS = CtauS, 1129 .strong_form = strong_form, 1130 .stabilization = stab, 1131 }; 1132 switch (problemChoice) { 1133 case NS_DENSITY_CURRENT: 1134 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1135 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1136 sizeof ctxNS); 1137 break; 1138 case NS_ADVECTION: 1139 case NS_ADVECTION2D: 1140 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1141 sizeof ctxAdvection2d); 1142 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1143 sizeof ctxAdvection2d); 1144 } 1145 1146 // Set up PETSc context 1147 // Set up units structure 1148 units->meter = meter; 1149 units->kilogram = kilogram; 1150 units->second = second; 1151 units->Kelvin = Kelvin; 1152 units->Pascal = Pascal; 1153 units->JperkgK = JperkgK; 1154 units->mpersquareds = mpersquareds; 1155 units->WpermK = WpermK; 1156 units->kgpercubicm = kgpercubicm; 1157 units->kgpersquaredms = kgpersquaredms; 1158 units->Joulepercubicm = Joulepercubicm; 1159 1160 // Set up user structure 1161 user->comm = comm; 1162 user->outputfreq = outputfreq; 1163 user->contsteps = contsteps; 1164 user->units = units; 1165 user->dm = dm; 1166 user->dmviz = dmviz; 1167 user->interpviz = interpviz; 1168 user->ceed = ceed; 1169 1170 // Calculate qdataVol and ICs 1171 // Set up state global and local vectors 1172 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1173 1174 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1175 1176 // Apply Setup Ceed Operators 1177 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1178 CeedOperatorApply(op_setupVol, xcorners, qdataVol, CEED_REQUEST_IMMEDIATE); 1179 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictqVol, basisqVol, restrictqdiVol, qdataVol, 1180 user->M); CHKERRQ(ierr); 1181 1182 ierr = ICs_PetscMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictqVol, 1183 &ctxSetup, 0.0); 1184 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1185 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1186 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1187 // slower execution. 1188 Vec Qbc; 1189 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1190 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1191 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1192 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1193 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1194 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1195 ierr = PetscObjectComposeFunction((PetscObject)dm, 1196 "DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_NS); CHKERRQ(ierr); 1197 } 1198 1199 MPI_Comm_rank(comm, &rank); 1200 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1201 // Gather initial Q values 1202 // In case of continuation of simulation, set up initial values from binary file 1203 if (contsteps) { // continue from existent solution 1204 PetscViewer viewer; 1205 char filepath[PETSC_MAX_PATH_LEN]; 1206 // Read input 1207 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1208 user->outputfolder); 1209 CHKERRQ(ierr); 1210 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1211 CHKERRQ(ierr); 1212 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1213 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1214 } else { 1215 //ierr = DMLocalToGlobal(dm, Qloc, INSERT_VALUES, Q);CHKERRQ(ierr); 1216 } 1217 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1218 1219 // Create and setup TS 1220 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1221 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1222 if (implicit) { 1223 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1224 if (user->op_ifunction) { 1225 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1226 } else { // Implicit integrators can fall back to using an RHSFunction 1227 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1228 } 1229 } else { 1230 if (!user->op_rhs) SETERRQ(comm,PETSC_ERR_ARG_NULL, 1231 "Problem does not provide RHSFunction"); 1232 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1233 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1234 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1235 } 1236 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1237 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1238 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1239 if (test) {ierr = TSSetMaxSteps(ts, 1); CHKERRQ(ierr);} 1240 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1241 ierr = TSAdaptSetStepLimits(adapt, 1242 1.e-12 * units->second, 1243 1.e2 * units->second); CHKERRQ(ierr); 1244 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1245 if (!contsteps) { // print initial condition 1246 if (!test) { 1247 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1248 } 1249 } else { // continue from time of last output 1250 PetscReal time; 1251 PetscInt count; 1252 PetscViewer viewer; 1253 char filepath[PETSC_MAX_PATH_LEN]; 1254 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1255 user->outputfolder); CHKERRQ(ierr); 1256 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1257 CHKERRQ(ierr); 1258 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1259 CHKERRQ(ierr); 1260 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1261 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1262 } 1263 if (!test) { 1264 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1265 } 1266 1267 // Solve 1268 start = MPI_Wtime(); 1269 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1270 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1271 cpu_time_used = MPI_Wtime() - start; 1272 ierr = TSGetSolveTime(ts,&ftime); CHKERRQ(ierr); 1273 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1274 comm); CHKERRQ(ierr); 1275 if (!test) { 1276 ierr = PetscPrintf(PETSC_COMM_WORLD, 1277 "Time taken for solution: %g\n", 1278 (double)cpu_time_used); CHKERRQ(ierr); 1279 } 1280 1281 // Get error 1282 if (problem->non_zero_time && !test) { 1283 Vec Qexact, Qexactloc; 1284 PetscReal norm; 1285 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1286 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1287 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1288 1289 ierr = ICs_PetscMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1290 restrictqVol, &ctxSetup, ftime); CHKERRQ(ierr); 1291 1292 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1293 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1294 CeedVectorDestroy(&q0ceed); 1295 ierr = PetscPrintf(PETSC_COMM_WORLD, 1296 "Max Error: %g\n", 1297 (double)norm); CHKERRQ(ierr); 1298 } 1299 1300 // Output Statistics 1301 ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr); 1302 if (!test) { 1303 ierr = PetscPrintf(PETSC_COMM_WORLD, 1304 "Time integrator took %D time steps to reach final time %g\n", 1305 steps,(double)ftime); CHKERRQ(ierr); 1306 } 1307 1308 // Clean up libCEED 1309 CeedVectorDestroy(&qdataVol); 1310 CeedVectorDestroy(&user->qceed); 1311 CeedVectorDestroy(&user->qdotceed); 1312 CeedVectorDestroy(&user->gceed); 1313 CeedVectorDestroy(&xcorners); 1314 CeedBasisDestroy(&basisqVol); 1315 CeedBasisDestroy(&basisxVol); 1316 CeedBasisDestroy(&basisxcVol); 1317 CeedElemRestrictionDestroy(&restrictqVol); 1318 CeedElemRestrictionDestroy(&restrictxVol); 1319 CeedElemRestrictionDestroy(&restrictqdiVol); 1320 CeedElemRestrictionDestroy(&restrictxcoordVol); 1321 CeedQFunctionDestroy(&qf_setupVol); 1322 CeedQFunctionDestroy(&qf_ics); 1323 CeedQFunctionDestroy(&qf_rhsVol); 1324 CeedQFunctionDestroy(&qf_ifunctionVol); 1325 CeedOperatorDestroy(&op_setupVol); 1326 CeedOperatorDestroy(&op_ics); 1327 CeedOperatorDestroy(&user->op_rhs); 1328 CeedOperatorDestroy(&user->op_ifunction); 1329 CeedDestroy(&ceed); 1330 1331 // Clean up PETSc 1332 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1333 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1334 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1335 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1336 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1337 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1338 ierr = PetscFree(units); CHKERRQ(ierr); 1339 ierr = PetscFree(user); CHKERRQ(ierr); 1340 return PetscFinalize(); 1341 } 1342 1343