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