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 1 32 // ./navierstokes -ceed /gpu/cuda -problem advection -degree 1 33 // 34 //TESTARGS -ceed {ceed_resource} -test explicit -degree 3 -dm_plex_box_faces 1,1,2 -units_kilogram 1e-9 -lx 125 -ly 125 -lz 250 -center 62.5,62.5,187.5 -rc 100. -thetaC -35. -ksp_atol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 35 //TESTARGS -ceed {ceed_resource} -test implicit_stab_none -degree 3 -dm_plex_box_faces 1,1,2 -units_kilogram 1e-9 -lx 125 -ly 125 -lz 250 -center 62.5,62.5,187.5 -rc 100. -thetaC -35. -ksp_atol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 -implicit -ts_type alpha 36 //TESTARGS -ceed {ceed_resource} -test implicit_stab_supg -degree 3 -dm_plex_box_faces 1,1,2 -units_kilogram 1e-9 -lx 125 -ly 125 -lz 250 -center 62.5,62.5,187.5 -rc 100. -thetaC -35. -ksp_atol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 -implicit -ts_type alpha -stab supg 37 38 /// @file 39 /// Navier-Stokes example using PETSc 40 41 const char help[] = "Solve Navier-Stokes using PETSc and libCEED\n"; 42 43 #include <ceed.h> 44 #include <petscdmplex.h> 45 #include <petscsys.h> 46 #include <petscts.h> 47 #include <stdbool.h> 48 #include "advection.h" 49 #include "advection2d.h" 50 #include "common.h" 51 #include "densitycurrent.h" 52 #include "setup-boundary.h" 53 54 #if PETSC_VERSION_LT(3,14,0) 55 # define DMPlexGetClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexGetClosureIndices(a,b,c,d,f,g,i) 56 # define DMPlexRestoreClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexRestoreClosureIndices(a,b,c,d,f,g,i) 57 #endif 58 59 #if PETSC_VERSION_LT(3,14,0) 60 # define DMAddBoundary(a,b,c,d,e,f,g,h,i,j,k,l) DMAddBoundary(a,b,c,d,e,f,g,h,j,k,l) 61 #endif 62 63 // MemType Options 64 static const char *const memTypes[] = { 65 "host", 66 "device", 67 "memType", "CEED_MEM_", NULL 68 }; 69 70 // Problem Options 71 typedef enum { 72 NS_DENSITY_CURRENT = 0, 73 NS_ADVECTION = 1, 74 NS_ADVECTION2D = 2, 75 } problemType; 76 static const char *const problemTypes[] = { 77 "density_current", 78 "advection", 79 "advection2d", 80 "problemType", "NS_", NULL 81 }; 82 83 // Wind Options for Advection 84 typedef enum { 85 ADVECTION_WIND_ROTATION = 0, 86 ADVECTION_WIND_TRANSLATION = 1, 87 } WindType; 88 static const char *const WindTypes[] = { 89 "rotation", 90 "translation", 91 "WindType", "ADVECTION_WIND_", NULL 92 }; 93 94 typedef enum { 95 STAB_NONE = 0, 96 STAB_SU = 1, // Streamline Upwind 97 STAB_SUPG = 2, // Streamline Upwind Petrov-Galerkin 98 } StabilizationType; 99 static const char *const StabilizationTypes[] = { 100 "none", 101 "SU", 102 "SUPG", 103 "StabilizationType", "STAB_", NULL 104 }; 105 106 // Test Options 107 typedef enum { 108 TEST_NONE = 0, // Non test mode 109 TEST_EXPLICIT = 1, // Explicit test 110 TEST_IMPLICIT_STAB_NONE = 2, // Implicit test no stab 111 TEST_IMPLICIT_STAB_SUPG = 3, // Implicit test supg stab 112 } testType; 113 static const char *const testTypes[] = { 114 "none", 115 "explicit", 116 "implicit_stab_none", 117 "implicit_stab_supg", 118 "testType", "TEST_", NULL 119 }; 120 121 // Tests specific data 122 typedef struct { 123 PetscScalar testtol; 124 const char *filepath; 125 } testData; 126 127 testData testOptions[] = { 128 [TEST_NONE] = { 129 .testtol = 0., 130 .filepath = NULL 131 }, 132 [TEST_EXPLICIT] = { 133 .testtol = 1E-5, 134 .filepath = "examples/fluids/tests-output/fluids-navierstokes-explicit.bin" 135 }, 136 [TEST_IMPLICIT_STAB_NONE] = { 137 .testtol = 5E-4, 138 .filepath = "examples/fluids/tests-output/fluids-navierstokes-implicit-stab-none.bin" 139 }, 140 [TEST_IMPLICIT_STAB_SUPG] = { 141 .testtol = 5E-4, 142 .filepath = "examples/fluids/tests-output/fluids-navierstokes-implicit-stab-supg.bin" 143 } 144 }; 145 146 // Problem specific data 147 typedef struct { 148 CeedInt dim, qdatasizeVol, qdatasizeSur; 149 CeedQFunctionUser setupVol, setupSur, ics, applyVol_rhs, applyVol_ifunction, 150 applySur; 151 PetscErrorCode (*bc)(PetscInt, PetscReal, const PetscReal[], PetscInt, 152 PetscScalar[], void *); 153 const char *setupVol_loc, *setupSur_loc, *ics_loc, *applyVol_rhs_loc, 154 *applyVol_ifunction_loc, *applySur_loc; 155 const bool non_zero_time; 156 } problemData; 157 158 problemData problemOptions[] = { 159 [NS_DENSITY_CURRENT] = { 160 .dim = 3, 161 .qdatasizeVol = 10, 162 .qdatasizeSur = 4, 163 .setupVol = Setup, 164 .setupVol_loc = Setup_loc, 165 .setupSur = SetupBoundary, 166 .setupSur_loc = SetupBoundary_loc, 167 .ics = ICsDC, 168 .ics_loc = ICsDC_loc, 169 .applyVol_rhs = DC, 170 .applyVol_rhs_loc = DC_loc, 171 .applyVol_ifunction = IFunction_DC, 172 .applyVol_ifunction_loc = IFunction_DC_loc, 173 .bc = Exact_DC, 174 .non_zero_time = PETSC_FALSE, 175 }, 176 [NS_ADVECTION] = { 177 .dim = 3, 178 .qdatasizeVol = 10, 179 .qdatasizeSur = 4, 180 .setupVol = Setup, 181 .setupVol_loc = Setup_loc, 182 .setupSur = SetupBoundary, 183 .setupSur_loc = SetupBoundary_loc, 184 .ics = ICsAdvection, 185 .ics_loc = ICsAdvection_loc, 186 .applyVol_rhs = Advection, 187 .applyVol_rhs_loc = Advection_loc, 188 .applyVol_ifunction = IFunction_Advection, 189 .applyVol_ifunction_loc = IFunction_Advection_loc, 190 .applySur = Advection_Sur, 191 .applySur_loc = Advection_Sur_loc, 192 .bc = Exact_Advection, 193 .non_zero_time = PETSC_FALSE, 194 }, 195 [NS_ADVECTION2D] = { 196 .dim = 2, 197 .qdatasizeVol = 5, 198 .qdatasizeSur = 3, 199 .setupVol = Setup2d, 200 .setupVol_loc = Setup2d_loc, 201 .setupSur = SetupBoundary2d, 202 .setupSur_loc = SetupBoundary2d_loc, 203 .ics = ICsAdvection2d, 204 .ics_loc = ICsAdvection2d_loc, 205 .applyVol_rhs = Advection2d, 206 .applyVol_rhs_loc = Advection2d_loc, 207 .applyVol_ifunction = IFunction_Advection2d, 208 .applyVol_ifunction_loc = IFunction_Advection2d_loc, 209 .applySur = Advection2d_Sur, 210 .applySur_loc = Advection2d_Sur_loc, 211 .bc = Exact_Advection2d, 212 .non_zero_time = PETSC_TRUE, 213 }, 214 }; 215 216 // PETSc user data 217 typedef struct User_ *User; 218 typedef struct Units_ *Units; 219 220 struct User_ { 221 MPI_Comm comm; 222 PetscInt outputfreq; 223 DM dm; 224 DM dmviz; 225 Mat interpviz; 226 Ceed ceed; 227 Units units; 228 CeedVector qceed, qdotceed, gceed; 229 CeedOperator op_rhs_vol, op_rhs, op_ifunction_vol, op_ifunction; 230 Vec M; 231 char outputdir[PETSC_MAX_PATH_LEN]; 232 PetscInt contsteps; 233 }; 234 235 struct Units_ { 236 // fundamental units 237 PetscScalar meter; 238 PetscScalar kilogram; 239 PetscScalar second; 240 PetscScalar Kelvin; 241 // derived units 242 PetscScalar Pascal; 243 PetscScalar JperkgK; 244 PetscScalar mpersquareds; 245 PetscScalar WpermK; 246 PetscScalar kgpercubicm; 247 PetscScalar kgpersquaredms; 248 PetscScalar Joulepercubicm; 249 PetscScalar Joule; 250 }; 251 252 typedef struct SimpleBC_ *SimpleBC; 253 struct SimpleBC_ { 254 PetscInt nwall, nslip[3]; 255 PetscInt walls[6], slips[3][6]; 256 PetscBool userbc; 257 }; 258 259 // Essential BC dofs are encoded in closure indices as -(i+1). 260 static PetscInt Involute(PetscInt i) { 261 return i >= 0 ? i : -(i+1); 262 } 263 264 // Utility function to create local CEED restriction 265 static PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 266 CeedInt height, DMLabel domainLabel, CeedInt value, 267 CeedElemRestriction *Erestrict) { 268 269 PetscSection section; 270 PetscInt p, Nelem, Ndof, *erestrict, eoffset, nfields, dim, depth; 271 DMLabel depthLabel; 272 IS depthIS, iterIS; 273 Vec Uloc; 274 const PetscInt *iterIndices; 275 PetscErrorCode ierr; 276 277 PetscFunctionBeginUser; 278 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 279 dim -= height; 280 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 281 ierr = PetscSectionGetNumFields(section, &nfields); CHKERRQ(ierr); 282 PetscInt ncomp[nfields], fieldoff[nfields+1]; 283 fieldoff[0] = 0; 284 for (PetscInt f=0; f<nfields; f++) { 285 ierr = PetscSectionGetFieldComponents(section, f, &ncomp[f]); CHKERRQ(ierr); 286 fieldoff[f+1] = fieldoff[f] + ncomp[f]; 287 } 288 289 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 290 ierr = DMPlexGetDepthLabel(dm, &depthLabel); CHKERRQ(ierr); 291 ierr = DMLabelGetStratumIS(depthLabel, depth - height, &depthIS); CHKERRQ(ierr); 292 if (domainLabel) { 293 IS domainIS; 294 ierr = DMLabelGetStratumIS(domainLabel, value, &domainIS); CHKERRQ(ierr); 295 if (domainIS) { // domainIS is non-empty 296 ierr = ISIntersect(depthIS, domainIS, &iterIS); CHKERRQ(ierr); 297 ierr = ISDestroy(&domainIS); CHKERRQ(ierr); 298 } else { // domainIS is NULL (empty) 299 iterIS = NULL; 300 } 301 ierr = ISDestroy(&depthIS); CHKERRQ(ierr); 302 } else { 303 iterIS = depthIS; 304 } 305 if (iterIS) { 306 ierr = ISGetLocalSize(iterIS, &Nelem); CHKERRQ(ierr); 307 ierr = ISGetIndices(iterIS, &iterIndices); CHKERRQ(ierr); 308 } else { 309 Nelem = 0; 310 iterIndices = NULL; 311 } 312 ierr = PetscMalloc1(Nelem*PetscPowInt(P, dim), &erestrict); CHKERRQ(ierr); 313 for (p=0,eoffset=0; p<Nelem; p++) { 314 PetscInt c = iterIndices[p]; 315 PetscInt numindices, *indices, nnodes; 316 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, 317 &numindices, &indices, NULL, NULL); 318 CHKERRQ(ierr); 319 bool flip = false; 320 if (height > 0) { 321 PetscInt numCells, numFaces, start = -1; 322 const PetscInt *orients, *faces, *cells; 323 ierr = DMPlexGetSupport(dm, c, &cells); CHKERRQ(ierr); 324 ierr = DMPlexGetSupportSize(dm, c, &numCells); CHKERRQ(ierr); 325 if (numCells != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 326 "Expected one cell in support of exterior face, but got %D cells", 327 numCells); 328 ierr = DMPlexGetCone(dm, cells[0], &faces); CHKERRQ(ierr); 329 ierr = DMPlexGetConeSize(dm, cells[0], &numFaces); CHKERRQ(ierr); 330 for (PetscInt i=0; i<numFaces; i++) {if (faces[i] == c) start = i;} 331 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, 332 "Could not find face %D in cone of its support", 333 c); 334 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients); CHKERRQ(ierr); 335 if (orients[start] < 0) flip = true; 336 } 337 if (numindices % fieldoff[nfields]) SETERRQ1(PETSC_COMM_SELF, 338 PETSC_ERR_ARG_INCOMP, "Number of closure indices not compatible with Cell %D", 339 c); 340 nnodes = numindices / fieldoff[nfields]; 341 for (PetscInt i=0; i<nnodes; i++) { 342 PetscInt ii = i; 343 if (flip) { 344 if (P == nnodes) ii = nnodes - 1 - i; 345 else if (P*P == nnodes) { 346 PetscInt row = i / P, col = i % P; 347 ii = row + col * P; 348 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, 349 "No support for flipping point with %D nodes != P (%D) or P^2", 350 nnodes, P); 351 } 352 // Check that indices are blocked by node and thus can be coalesced as a single field with 353 // fieldoff[nfields] = sum(ncomp) components. 354 for (PetscInt f=0; f<nfields; f++) { 355 for (PetscInt j=0; j<ncomp[f]; j++) { 356 if (Involute(indices[fieldoff[f]*nnodes + ii*ncomp[f] + j]) 357 != Involute(indices[ii*ncomp[0]]) + fieldoff[f] + j) 358 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 359 "Cell %D closure indices not interlaced for node %D field %D component %D", 360 c, ii, f, j); 361 } 362 } 363 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 364 PetscInt loc = Involute(indices[ii*ncomp[0]]); 365 erestrict[eoffset++] = loc; 366 } 367 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, 368 &numindices, &indices, NULL, NULL); 369 CHKERRQ(ierr); 370 } 371 if (eoffset != Nelem*PetscPowInt(P, dim)) 372 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 373 "ElemRestriction of size (%D,%D) initialized %D nodes", Nelem, 374 PetscPowInt(P, dim),eoffset); 375 if (iterIS) { 376 ierr = ISRestoreIndices(iterIS, &iterIndices); CHKERRQ(ierr); 377 } 378 ierr = ISDestroy(&iterIS); CHKERRQ(ierr); 379 380 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 381 ierr = VecGetLocalSize(Uloc, &Ndof); CHKERRQ(ierr); 382 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 383 CeedElemRestrictionCreate(ceed, Nelem, PetscPowInt(P, dim), fieldoff[nfields], 384 1, Ndof, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 385 Erestrict); 386 ierr = PetscFree(erestrict); CHKERRQ(ierr); 387 PetscFunctionReturn(0); 388 } 389 390 // Utility function to get Ceed Restriction for each domain 391 static PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 392 DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, CeedInt qdatasize, 393 CeedElemRestriction *restrictq, CeedElemRestriction *restrictx, 394 CeedElemRestriction *restrictqdi) { 395 396 DM dmcoord; 397 CeedInt dim, localNelem; 398 CeedInt Qdim; 399 PetscErrorCode ierr; 400 401 PetscFunctionBeginUser; 402 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 403 dim -= height; 404 Qdim = CeedIntPow(Q, dim); 405 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 406 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 407 CHKERRQ(ierr); 408 ierr = CreateRestrictionFromPlex(ceed, dm, P, height, domainLabel, value, 409 restrictq); 410 CHKERRQ(ierr); 411 ierr = CreateRestrictionFromPlex(ceed, dmcoord, 2, height, domainLabel, value, 412 restrictx); 413 CHKERRQ(ierr); 414 CeedElemRestrictionGetNumElements(*restrictq, &localNelem); 415 CeedElemRestrictionCreateStrided(ceed, localNelem, Qdim, 416 qdatasize, qdatasize*localNelem*Qdim, 417 CEED_STRIDES_BACKEND, restrictqdi); 418 PetscFunctionReturn(0); 419 } 420 421 // Utility function to create CEED Composite Operator for the entire domain 422 static PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, 423 WindType wind_type, CeedOperator op_applyVol, CeedQFunction qf_applySur, 424 CeedQFunction qf_setupSur, CeedInt height, CeedInt numP_Sur, CeedInt numQ_Sur, 425 CeedInt qdatasizeSur, CeedInt NqptsSur, CeedBasis basisxSur, 426 CeedBasis basisqSur, CeedOperator *op_apply) { 427 428 CeedInt dim, nFace; 429 PetscInt lsize, localNelemSur[6]; 430 Vec Xloc; 431 CeedVector xcorners, qdataSur[6]; 432 CeedOperator op_setupSur[6], op_applySur[6]; 433 CeedElemRestriction restrictxSur[6], restrictqSur[6], restrictqdiSur[6]; 434 DMLabel domainLabel; 435 PetscScalar *x; 436 PetscErrorCode ierr; 437 438 PetscFunctionBeginUser; 439 // Composite Operaters 440 CeedCompositeOperatorCreate(ceed, op_apply); 441 CeedCompositeOperatorAddSub(*op_apply, 442 op_applyVol); // Apply a Sub-Operator for the volume 443 444 if (wind_type == ADVECTION_WIND_TRANSLATION) { 445 bc->nwall = 0; 446 bc->nslip[0] = bc->nslip[1] = bc->nslip[2] = 0; 447 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 448 ierr = VecGetLocalSize(Xloc, &lsize); CHKERRQ(ierr); 449 ierr = CeedVectorCreate(ceed, lsize, &xcorners); CHKERRQ(ierr); 450 ierr = VecGetArray(Xloc, &x); CHKERRQ(ierr); 451 CeedVectorSetArray(xcorners, CEED_MEM_HOST, CEED_USE_POINTER, x); 452 ierr = DMGetLabel(dm, "Face Sets", &domainLabel); CHKERRQ(ierr); 453 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 454 if (dim == 2) nFace = 4; 455 if (dim == 3) nFace = 6; 456 457 // Create CEED Operator for each boundary face 458 for (CeedInt i=0; i<nFace; i++) { 459 ierr = GetRestrictionForDomain(ceed, dm, height, domainLabel, i+1, numP_Sur, 460 numQ_Sur, qdatasizeSur, &restrictqSur[i], 461 &restrictxSur[i], &restrictqdiSur[i]); 462 CHKERRQ(ierr); 463 // Create the CEED vectors that will be needed in Boundary setup 464 CeedElemRestrictionGetNumElements(restrictqSur[i], &localNelemSur[i]); 465 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, 466 &qdataSur[i]); 467 // Create the operator that builds the quadrature data for the Boundary operator 468 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur[i]); 469 CeedOperatorSetField(op_setupSur[i], "dx", restrictxSur[i], basisxSur, 470 CEED_VECTOR_ACTIVE); 471 CeedOperatorSetField(op_setupSur[i], "weight", CEED_ELEMRESTRICTION_NONE, 472 basisxSur, CEED_VECTOR_NONE); 473 CeedOperatorSetField(op_setupSur[i], "qdataSur", restrictqdiSur[i], 474 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 475 // Create Boundary operator 476 CeedOperatorCreate(ceed, qf_applySur, NULL, NULL, &op_applySur[i]); 477 CeedOperatorSetField(op_applySur[i], "q", restrictqSur[i], basisqSur, 478 CEED_VECTOR_ACTIVE); 479 CeedOperatorSetField(op_applySur[i], "qdataSur", restrictqdiSur[i], 480 CEED_BASIS_COLLOCATED, qdataSur[i]); 481 CeedOperatorSetField(op_applySur[i], "x", restrictxSur[i], basisxSur, 482 xcorners); 483 CeedOperatorSetField(op_applySur[i], "v", restrictqSur[i], basisqSur, 484 CEED_VECTOR_ACTIVE); 485 // Apply CEED operator for Boundary setup 486 CeedOperatorApply(op_setupSur[i], xcorners, qdataSur[i], 487 CEED_REQUEST_IMMEDIATE); 488 // Apply Sub-Operator for the Boundary 489 CeedCompositeOperatorAddSub(*op_apply, op_applySur[i]); 490 } 491 CeedVectorDestroy(&xcorners); 492 } 493 PetscFunctionReturn(0); 494 } 495 496 static int CreateVectorFromPetscVec(Ceed ceed, Vec p, CeedVector *v) { 497 PetscErrorCode ierr; 498 PetscInt m; 499 500 PetscFunctionBeginUser; 501 ierr = VecGetLocalSize(p, &m); CHKERRQ(ierr); 502 ierr = CeedVectorCreate(ceed, m, v); CHKERRQ(ierr); 503 PetscFunctionReturn(0); 504 } 505 506 static int VectorPlacePetscVec(CeedVector c, Vec p) { 507 PetscErrorCode ierr; 508 PetscInt mceed, mpetsc; 509 PetscScalar *a; 510 511 PetscFunctionBeginUser; 512 ierr = CeedVectorGetLength(c, &mceed); CHKERRQ(ierr); 513 ierr = VecGetLocalSize(p, &mpetsc); CHKERRQ(ierr); 514 if (mceed != mpetsc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 515 "Cannot place PETSc Vec of length %D in CeedVector of length %D", 516 mpetsc, mceed); 517 ierr = VecGetArray(p, &a); CHKERRQ(ierr); 518 CeedVectorSetArray(c, CEED_MEM_HOST, CEED_USE_POINTER, a); 519 PetscFunctionReturn(0); 520 } 521 522 static PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 523 PetscBool insertEssential, Vec Qloc, PetscReal time, Vec faceGeomFVM, 524 Vec cellGeomFVM, Vec gradFVM) { 525 PetscErrorCode ierr; 526 Vec Qbc; 527 528 PetscFunctionBegin; 529 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 530 ierr = VecAXPY(Qloc, 1., Qbc); CHKERRQ(ierr); 531 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 532 PetscFunctionReturn(0); 533 } 534 535 // This is the RHS of the ODE, given as u_t = G(t,u) 536 // This function takes in a state vector Q and writes into G 537 static PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *userData) { 538 PetscErrorCode ierr; 539 User user = *(User *)userData; 540 PetscScalar *q, *g; 541 Vec Qloc, Gloc; 542 543 // Global-to-local 544 PetscFunctionBeginUser; 545 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 546 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 547 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 548 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 549 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 550 NULL, NULL, NULL); CHKERRQ(ierr); 551 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 552 553 // Ceed Vectors 554 ierr = VecGetArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 555 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 556 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); 557 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 558 559 // Apply CEED operator 560 CeedOperatorApply(user->op_rhs, user->qceed, user->gceed, 561 CEED_REQUEST_IMMEDIATE); 562 563 // Restore vectors 564 ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 565 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 566 567 ierr = VecZeroEntries(G); CHKERRQ(ierr); 568 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 569 570 // Inverse of the lumped mass matrix 571 ierr = VecPointwiseMult(G, G, user->M); // M is Minv 572 CHKERRQ(ierr); 573 574 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 575 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 576 PetscFunctionReturn(0); 577 } 578 579 static PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Qdot, Vec G, 580 void *userData) { 581 PetscErrorCode ierr; 582 User user = *(User *)userData; 583 const PetscScalar *q, *qdot; 584 PetscScalar *g; 585 Vec Qloc, Qdotloc, Gloc; 586 587 // Global-to-local 588 PetscFunctionBeginUser; 589 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 590 ierr = DMGetLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 591 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 592 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 593 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 594 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 595 NULL, NULL, NULL); CHKERRQ(ierr); 596 ierr = VecZeroEntries(Qdotloc); CHKERRQ(ierr); 597 ierr = DMGlobalToLocal(user->dm, Qdot, INSERT_VALUES, Qdotloc); CHKERRQ(ierr); 598 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 599 600 // Ceed Vectors 601 ierr = VecGetArrayRead(Qloc, &q); CHKERRQ(ierr); 602 ierr = VecGetArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 603 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 604 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, 605 (PetscScalar *)q); 606 CeedVectorSetArray(user->qdotceed, CEED_MEM_HOST, CEED_USE_POINTER, 607 (PetscScalar *)qdot); 608 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 609 610 // Apply CEED operator 611 CeedOperatorApply(user->op_ifunction, user->qceed, user->gceed, 612 CEED_REQUEST_IMMEDIATE); 613 614 // Restore vectors 615 ierr = VecRestoreArrayRead(Qloc, &q); CHKERRQ(ierr); 616 ierr = VecRestoreArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 617 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 618 619 ierr = VecZeroEntries(G); CHKERRQ(ierr); 620 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 621 622 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 623 ierr = DMRestoreLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 624 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 625 PetscFunctionReturn(0); 626 } 627 628 // User provided TS Monitor 629 static PetscErrorCode TSMonitor_NS(TS ts, PetscInt stepno, PetscReal time, 630 Vec Q, void *ctx) { 631 User user = ctx; 632 Vec Qloc; 633 char filepath[PETSC_MAX_PATH_LEN]; 634 PetscViewer viewer; 635 PetscErrorCode ierr; 636 637 // Set up output 638 PetscFunctionBeginUser; 639 // Print every 'outputfreq' steps 640 if (stepno % user->outputfreq != 0) 641 PetscFunctionReturn(0); 642 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 643 ierr = PetscObjectSetName((PetscObject)Qloc, "StateVec"); CHKERRQ(ierr); 644 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 645 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 646 647 // Output 648 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-%03D.vtu", 649 user->outputdir, stepno + user->contsteps); 650 CHKERRQ(ierr); 651 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), filepath, 652 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 653 ierr = VecView(Qloc, viewer); CHKERRQ(ierr); 654 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 655 if (user->dmviz) { 656 Vec Qrefined, Qrefined_loc; 657 char filepath_refined[PETSC_MAX_PATH_LEN]; 658 PetscViewer viewer_refined; 659 660 ierr = DMGetGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 661 ierr = DMGetLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 662 ierr = PetscObjectSetName((PetscObject)Qrefined_loc, "Refined"); 663 CHKERRQ(ierr); 664 ierr = MatInterpolate(user->interpviz, Q, Qrefined); CHKERRQ(ierr); 665 ierr = VecZeroEntries(Qrefined_loc); CHKERRQ(ierr); 666 ierr = DMGlobalToLocal(user->dmviz, Qrefined, INSERT_VALUES, Qrefined_loc); 667 CHKERRQ(ierr); 668 ierr = PetscSNPrintf(filepath_refined, sizeof filepath_refined, 669 "%s/nsrefined-%03D.vtu", 670 user->outputdir, stepno + user->contsteps); 671 CHKERRQ(ierr); 672 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Qrefined), 673 filepath_refined, 674 FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 675 ierr = VecView(Qrefined_loc, viewer_refined); CHKERRQ(ierr); 676 ierr = DMRestoreLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 677 ierr = DMRestoreGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 678 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 679 } 680 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 681 682 // Save data in a binary file for continuation of simulations 683 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 684 user->outputdir); CHKERRQ(ierr); 685 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 686 CHKERRQ(ierr); 687 ierr = VecView(Q, viewer); CHKERRQ(ierr); 688 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 689 690 // Save time stamp 691 // Dimensionalize time back 692 time /= user->units->second; 693 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 694 user->outputdir); CHKERRQ(ierr); 695 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 696 CHKERRQ(ierr); 697 #if PETSC_VERSION_GE(3,13,0) 698 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 699 #else 700 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 701 #endif 702 CHKERRQ(ierr); 703 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 704 705 PetscFunctionReturn(0); 706 } 707 708 static PetscErrorCode ICs_FixMultiplicity(CeedOperator op_ics, 709 CeedVector xcorners, CeedVector q0ceed, DM dm, Vec Qloc, Vec Q, 710 CeedElemRestriction restrictq, CeedQFunctionContext ctxSetup, CeedScalar time) { 711 PetscErrorCode ierr; 712 CeedVector multlvec; 713 Vec Multiplicity, MultiplicityLoc; 714 715 SetupContext ctxSetupData; 716 CeedQFunctionContextGetData(ctxSetup, CEED_MEM_HOST, (void **)&ctxSetupData); 717 ctxSetupData->time = time; 718 CeedQFunctionContextRestoreData(ctxSetup, (void **)&ctxSetupData); 719 720 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 721 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 722 CeedOperatorApply(op_ics, xcorners, q0ceed, CEED_REQUEST_IMMEDIATE); 723 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 724 ierr = DMLocalToGlobal(dm, Qloc, ADD_VALUES, Q); CHKERRQ(ierr); 725 726 // Fix multiplicity for output of ICs 727 ierr = DMGetLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 728 CeedElemRestrictionCreateVector(restrictq, &multlvec, NULL); 729 ierr = VectorPlacePetscVec(multlvec, MultiplicityLoc); CHKERRQ(ierr); 730 CeedElemRestrictionGetMultiplicity(restrictq, multlvec); 731 CeedVectorDestroy(&multlvec); 732 ierr = DMGetGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 733 ierr = VecZeroEntries(Multiplicity); CHKERRQ(ierr); 734 ierr = DMLocalToGlobal(dm, MultiplicityLoc, ADD_VALUES, Multiplicity); 735 CHKERRQ(ierr); 736 ierr = VecPointwiseDivide(Q, Q, Multiplicity); CHKERRQ(ierr); 737 ierr = VecPointwiseDivide(Qloc, Qloc, MultiplicityLoc); CHKERRQ(ierr); 738 ierr = DMRestoreLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 739 ierr = DMRestoreGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 740 741 PetscFunctionReturn(0); 742 } 743 744 static PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, 745 CeedElemRestriction restrictq, CeedBasis basisq, 746 CeedElemRestriction restrictqdi, CeedVector qdata, Vec M) { 747 PetscErrorCode ierr; 748 CeedQFunction qf_mass; 749 CeedOperator op_mass; 750 CeedVector mceed; 751 Vec Mloc; 752 CeedInt ncompq, qdatasize; 753 754 PetscFunctionBeginUser; 755 CeedElemRestrictionGetNumComponents(restrictq, &ncompq); 756 CeedElemRestrictionGetNumComponents(restrictqdi, &qdatasize); 757 // Create the Q-function that defines the action of the mass operator 758 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 759 CeedQFunctionAddInput(qf_mass, "q", ncompq, CEED_EVAL_INTERP); 760 CeedQFunctionAddInput(qf_mass, "qdata", qdatasize, CEED_EVAL_NONE); 761 CeedQFunctionAddOutput(qf_mass, "v", ncompq, CEED_EVAL_INTERP); 762 763 // Create the mass operator 764 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 765 CeedOperatorSetField(op_mass, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 766 CeedOperatorSetField(op_mass, "qdata", restrictqdi, 767 CEED_BASIS_COLLOCATED, qdata); 768 CeedOperatorSetField(op_mass, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 769 770 ierr = DMGetLocalVector(dm, &Mloc); CHKERRQ(ierr); 771 ierr = VecZeroEntries(Mloc); CHKERRQ(ierr); 772 CeedElemRestrictionCreateVector(restrictq, &mceed, NULL); 773 ierr = VectorPlacePetscVec(mceed, Mloc); CHKERRQ(ierr); 774 775 { 776 // Compute a lumped mass matrix 777 CeedVector onesvec; 778 CeedElemRestrictionCreateVector(restrictq, &onesvec, NULL); 779 CeedVectorSetValue(onesvec, 1.0); 780 CeedOperatorApply(op_mass, onesvec, mceed, CEED_REQUEST_IMMEDIATE); 781 CeedVectorDestroy(&onesvec); 782 CeedOperatorDestroy(&op_mass); 783 CeedVectorDestroy(&mceed); 784 } 785 CeedQFunctionDestroy(&qf_mass); 786 787 ierr = VecZeroEntries(M); CHKERRQ(ierr); 788 ierr = DMLocalToGlobal(dm, Mloc, ADD_VALUES, M); CHKERRQ(ierr); 789 ierr = DMRestoreLocalVector(dm, &Mloc); CHKERRQ(ierr); 790 791 // Invert diagonally lumped mass vector for RHS function 792 ierr = VecReciprocal(M); CHKERRQ(ierr); 793 PetscFunctionReturn(0); 794 } 795 796 static PetscErrorCode SetUpDM(DM dm, problemData *problem, PetscInt degree, 797 SimpleBC bc, void *ctxSetupData) { 798 PetscErrorCode ierr; 799 800 PetscFunctionBeginUser; 801 { 802 // Configure the finite element space and boundary conditions 803 PetscFE fe; 804 PetscInt ncompq = 5; 805 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, ncompq, 806 PETSC_FALSE, degree, PETSC_DECIDE, 807 &fe); CHKERRQ(ierr); 808 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 809 ierr = DMAddField(dm, NULL,(PetscObject)fe); CHKERRQ(ierr); 810 ierr = DMCreateDS(dm); CHKERRQ(ierr); 811 { 812 PetscInt comps[1] = {1}; 813 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 814 1, comps, (void(*)(void))NULL, NULL, bc->nslip[0], 815 bc->slips[0], ctxSetupData); CHKERRQ(ierr); 816 comps[0] = 2; 817 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 818 1, comps, (void(*)(void))NULL, NULL, bc->nslip[1], 819 bc->slips[1], ctxSetupData); CHKERRQ(ierr); 820 comps[0] = 3; 821 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 822 1, comps, (void(*)(void))NULL, NULL, bc->nslip[2], 823 bc->slips[2], ctxSetupData); CHKERRQ(ierr); 824 } 825 if (bc->userbc == PETSC_TRUE) { 826 for (PetscInt c = 0; c < 3; c++) { 827 for (PetscInt s = 0; s < bc->nslip[c]; s++) { 828 for (PetscInt w = 0; w < bc->nwall; w++) { 829 if (bc->slips[c][s] == bc->walls[w]) 830 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, 831 "Boundary condition already set on face %D!\n", 832 bc->walls[w]); 833 834 } 835 } 836 } 837 } 838 // Wall boundary conditions are zero energy density and zero flux for 839 // velocity in advection/advection2d, and zero velocity and zero flux 840 // for mass density and energy density in density_current 841 { 842 if (problem->bc == Exact_Advection || problem->bc == Exact_Advection2d) { 843 PetscInt comps[1] = {4}; 844 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 845 1, comps, (void(*)(void))problem->bc, NULL, 846 bc->nwall, bc->walls, ctxSetupData); CHKERRQ(ierr); 847 } else if (problem->bc == Exact_DC) { 848 PetscInt comps[3] = {1, 2, 3}; 849 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 850 3, comps, (void(*)(void))problem->bc, NULL, 851 bc->nwall, bc->walls, ctxSetupData); CHKERRQ(ierr); 852 } else 853 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, 854 "Undefined boundary conditions for this problem"); 855 } 856 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 857 CHKERRQ(ierr); 858 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 859 } 860 { 861 // Empty name for conserved field (because there is only one field) 862 PetscSection section; 863 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 864 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 865 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 866 CHKERRQ(ierr); 867 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 868 CHKERRQ(ierr); 869 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 870 CHKERRQ(ierr); 871 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 872 CHKERRQ(ierr); 873 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 874 CHKERRQ(ierr); 875 } 876 PetscFunctionReturn(0); 877 } 878 879 int main(int argc, char **argv) { 880 PetscInt ierr; 881 MPI_Comm comm; 882 DM dm, dmcoord, dmviz; 883 Mat interpviz; 884 TS ts; 885 TSAdapt adapt; 886 User user; 887 Units units; 888 char ceedresource[4096] = "/cpu/self"; 889 PetscInt localNelemVol, lnodes, gnodes, steps; 890 const PetscInt ncompq = 5; 891 PetscMPIInt rank; 892 PetscScalar ftime; 893 Vec Q, Qloc, Xloc; 894 Ceed ceed; 895 CeedInt numP, numQ; 896 CeedVector xcorners, qdata, q0ceed; 897 CeedBasis basisx, basisxc, basisq; 898 CeedElemRestriction restrictx, restrictq, restrictqdi; 899 CeedQFunction qf_setupVol, qf_ics, qf_rhsVol, qf_ifunctionVol; 900 CeedQFunctionContext ctxSetup, ctxNS, ctxAdvection2d, ctxSurface; 901 CeedOperator op_setupVol, op_ics; 902 CeedScalar Rd; 903 CeedMemType memtyperequested; 904 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 905 kgpersquaredms, Joulepercubicm, Joule; 906 problemType problemChoice; 907 problemData *problem = NULL; 908 WindType wind_type; 909 StabilizationType stab; 910 testType testChoice; 911 testData *test = NULL; 912 PetscBool implicit; 913 PetscInt viz_refine = 0; 914 struct SimpleBC_ bc = { 915 .nslip = {2, 2, 2}, 916 .slips = {{5, 6}, {3, 4}, {1, 2}} 917 }; 918 double start, cpu_time_used; 919 // Check PETSc CUDA support 920 PetscBool petschavecuda, setmemtyperequest = PETSC_FALSE; 921 // *INDENT-OFF* 922 #ifdef PETSC_HAVE_CUDA 923 petschavecuda = PETSC_TRUE; 924 #else 925 petschavecuda = PETSC_FALSE; 926 #endif 927 // *INDENT-ON* 928 929 // Create the libCEED contexts 930 PetscScalar meter = 1e-2; // 1 meter in scaled length units 931 PetscScalar second = 1e-2; // 1 second in scaled time units 932 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 933 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 934 CeedScalar theta0 = 300.; // K 935 CeedScalar thetaC = -15.; // K 936 CeedScalar P0 = 1.e5; // Pa 937 CeedScalar E_wind = 1.e6; // J 938 CeedScalar N = 0.01; // 1/s 939 CeedScalar cv = 717.; // J/(kg K) 940 CeedScalar cp = 1004.; // J/(kg K) 941 CeedScalar g = 9.81; // m/s^2 942 CeedScalar lambda = -2./3.; // - 943 CeedScalar mu = 75.; // Pa s, dynamic viscosity 944 // mu = 75 is not physical for air, but is good for numerical stability 945 CeedScalar k = 0.02638; // W/(m K) 946 CeedScalar CtauS = 0.; // dimensionless 947 CeedScalar strong_form = 0.; // [0,1] 948 PetscScalar lx = 8000.; // m 949 PetscScalar ly = 8000.; // m 950 PetscScalar lz = 4000.; // m 951 CeedScalar rc = 1000.; // m (Radius of bubble) 952 PetscScalar resx = 1000.; // m (resolution in x) 953 PetscScalar resy = 1000.; // m (resolution in y) 954 PetscScalar resz = 1000.; // m (resolution in z) 955 PetscInt outputfreq = 10; // - 956 PetscInt contsteps = 0; // - 957 PetscInt degree = 1; // - 958 PetscInt qextra = 2; // - 959 PetscInt qextraSur = 2; // - 960 PetscReal center[3], dc_axis[3] = {0, 0, 0}, wind[3] = {1., 0, 0}; 961 962 ierr = PetscInitialize(&argc, &argv, NULL, help); 963 if (ierr) return ierr; 964 965 // Allocate PETSc context 966 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 967 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 968 969 // Parse command line options 970 comm = PETSC_COMM_WORLD; 971 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 972 NULL); CHKERRQ(ierr); 973 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 974 NULL, ceedresource, ceedresource, 975 sizeof(ceedresource), NULL); CHKERRQ(ierr); 976 testChoice = TEST_NONE; 977 ierr = PetscOptionsEnum("-test", "Run tests", NULL, 978 testTypes, (PetscEnum)testChoice, 979 (PetscEnum *)&testChoice, 980 NULL); CHKERRQ(ierr); 981 test = &testOptions[testChoice]; 982 problemChoice = NS_DENSITY_CURRENT; 983 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 984 problemTypes, (PetscEnum)problemChoice, 985 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 986 problem = &problemOptions[problemChoice]; 987 ierr = PetscOptionsEnum("-problem_advection_wind", "Wind type in Advection", 988 NULL, WindTypes, 989 (PetscEnum)(wind_type = ADVECTION_WIND_ROTATION), 990 (PetscEnum *)&wind_type, NULL); CHKERRQ(ierr); 991 PetscInt n = problem->dim; 992 PetscBool userWind; 993 ierr = PetscOptionsRealArray("-problem_advection_wind_translation", 994 "Constant wind vector", 995 NULL, wind, &n, &userWind); CHKERRQ(ierr); 996 if (wind_type == ADVECTION_WIND_ROTATION && userWind) { 997 ierr = PetscPrintf(comm, 998 "Warning! Use -problem_advection_wind_translation only with -problem_advection_wind translation\n"); 999 CHKERRQ(ierr); 1000 } 1001 if (wind_type == ADVECTION_WIND_TRANSLATION 1002 && problemChoice == NS_DENSITY_CURRENT) { 1003 SETERRQ(comm, PETSC_ERR_ARG_INCOMP, 1004 "-problem_advection_wind translation is not defined for -problem density_current"); 1005 } 1006 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 1007 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 1008 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 1009 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 1010 NULL, implicit=PETSC_FALSE, &implicit, NULL); 1011 CHKERRQ(ierr); 1012 if (!implicit && stab != STAB_NONE) { 1013 ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); 1014 CHKERRQ(ierr); 1015 } 1016 { 1017 PetscInt len; 1018 PetscBool flg; 1019 ierr = PetscOptionsIntArray("-bc_wall", 1020 "Use wall boundary conditions on this list of faces", 1021 NULL, bc.walls, 1022 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 1023 &len), &flg); CHKERRQ(ierr); 1024 if (flg) { 1025 bc.nwall = len; 1026 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 1027 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1028 } 1029 for (PetscInt j=0; j<3; j++) { 1030 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 1031 ierr = PetscOptionsIntArray(flags[j], 1032 "Use slip boundary conditions on this list of faces", 1033 NULL, bc.slips[j], 1034 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 1035 &len), &flg); 1036 CHKERRQ(ierr); 1037 if (flg) { 1038 bc.nslip[j] = len; 1039 bc.userbc = PETSC_TRUE; 1040 } 1041 } 1042 } 1043 ierr = PetscOptionsInt("-viz_refine", 1044 "Regular refinement levels for visualization", 1045 NULL, viz_refine, &viz_refine, NULL); 1046 CHKERRQ(ierr); 1047 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1048 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1049 meter = fabs(meter); 1050 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1051 NULL, second, &second, NULL); CHKERRQ(ierr); 1052 second = fabs(second); 1053 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1054 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1055 kilogram = fabs(kilogram); 1056 ierr = PetscOptionsScalar("-units_Kelvin", 1057 "1 Kelvin in scaled temperature units", 1058 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1059 Kelvin = fabs(Kelvin); 1060 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1061 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1062 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1063 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1064 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1065 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1066 ierr = PetscOptionsScalar("-E_wind", "Total energy of inflow wind", 1067 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 1068 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1069 NULL, N, &N, NULL); CHKERRQ(ierr); 1070 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1071 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1072 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1073 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1074 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1075 NULL, g, &g, NULL); CHKERRQ(ierr); 1076 ierr = PetscOptionsScalar("-lambda", 1077 "Stokes hypothesis second viscosity coefficient", 1078 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1079 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1080 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1081 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1082 NULL, k, &k, NULL); CHKERRQ(ierr); 1083 ierr = PetscOptionsScalar("-CtauS", 1084 "Scale coefficient for tau (nondimensional)", 1085 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1086 if (stab == STAB_NONE && CtauS != 0) { 1087 ierr = PetscPrintf(comm, 1088 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1089 CHKERRQ(ierr); 1090 } 1091 ierr = PetscOptionsScalar("-strong_form", 1092 "Strong (1) or weak/integrated by parts (0) advection residual", 1093 NULL, strong_form, &strong_form, NULL); 1094 CHKERRQ(ierr); 1095 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1096 ierr = PetscPrintf(comm, 1097 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1098 CHKERRQ(ierr); 1099 } 1100 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1101 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1102 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1103 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1104 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1105 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1106 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1107 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1108 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1109 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1110 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1111 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1112 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1113 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1114 n = problem->dim; 1115 center[0] = 0.5 * lx; 1116 center[1] = 0.5 * ly; 1117 center[2] = 0.5 * lz; 1118 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1119 NULL, center, &n, NULL); CHKERRQ(ierr); 1120 n = problem->dim; 1121 ierr = PetscOptionsRealArray("-dc_axis", 1122 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1123 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1124 { 1125 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1126 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1127 if (norm > 0) { 1128 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1129 } 1130 } 1131 ierr = PetscOptionsInt("-output_freq", 1132 "Frequency of output, in number of steps", 1133 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1134 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1135 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1136 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1137 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1138 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1139 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1140 PetscBool userQextraSur; 1141 ierr = PetscOptionsInt("-qextra_boundary", 1142 "Number of extra quadrature points on in/outflow faces", 1143 NULL, qextraSur, &qextraSur, &userQextraSur); 1144 CHKERRQ(ierr); 1145 if ((wind_type == ADVECTION_WIND_ROTATION 1146 || problemChoice == NS_DENSITY_CURRENT) && userQextraSur) { 1147 ierr = PetscPrintf(comm, 1148 "Warning! Use -qextra_boundary only with -problem_advection_wind translation\n"); 1149 CHKERRQ(ierr); 1150 } 1151 ierr = PetscStrncpy(user->outputdir, ".", 2); CHKERRQ(ierr); 1152 ierr = PetscOptionsString("-output_dir", "Output directory", 1153 NULL, user->outputdir, user->outputdir, 1154 sizeof(user->outputdir), NULL); CHKERRQ(ierr); 1155 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1156 ierr = PetscOptionsEnum("-memtype", 1157 "CEED MemType requested", NULL, 1158 memTypes, (PetscEnum)memtyperequested, 1159 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1160 CHKERRQ(ierr); 1161 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1162 1163 // Define derived units 1164 Pascal = kilogram / (meter * PetscSqr(second)); 1165 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1166 mpersquareds = meter / PetscSqr(second); 1167 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1168 kgpercubicm = kilogram / pow(meter,3); 1169 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1170 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1171 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 1172 1173 // Scale variables to desired units 1174 theta0 *= Kelvin; 1175 thetaC *= Kelvin; 1176 P0 *= Pascal; 1177 E_wind *= Joule; 1178 N *= (1./second); 1179 cv *= JperkgK; 1180 cp *= JperkgK; 1181 Rd = cp - cv; 1182 g *= mpersquareds; 1183 mu *= Pascal * second; 1184 k *= WpermK; 1185 lx = fabs(lx) * meter; 1186 ly = fabs(ly) * meter; 1187 lz = fabs(lz) * meter; 1188 rc = fabs(rc) * meter; 1189 resx = fabs(resx) * meter; 1190 resy = fabs(resy) * meter; 1191 resz = fabs(resz) * meter; 1192 for (int i=0; i<3; i++) center[i] *= meter; 1193 1194 const CeedInt dim = problem->dim, ncompx = problem->dim, 1195 qdatasizeVol = problem->qdatasizeVol; 1196 // Set up the libCEED context 1197 struct SetupContext_ ctxSetupData = { 1198 .theta0 = theta0, 1199 .thetaC = thetaC, 1200 .P0 = P0, 1201 .N = N, 1202 .cv = cv, 1203 .cp = cp, 1204 .Rd = Rd, 1205 .g = g, 1206 .rc = rc, 1207 .lx = lx, 1208 .ly = ly, 1209 .lz = lz, 1210 .center[0] = center[0], 1211 .center[1] = center[1], 1212 .center[2] = center[2], 1213 .dc_axis[0] = dc_axis[0], 1214 .dc_axis[1] = dc_axis[1], 1215 .dc_axis[2] = dc_axis[2], 1216 .wind[0] = wind[0], 1217 .wind[1] = wind[1], 1218 .wind[2] = wind[2], 1219 .time = 0, 1220 .wind_type = wind_type, 1221 }; 1222 1223 // Create the mesh 1224 { 1225 const PetscReal scale[3] = {lx, ly, lz}; 1226 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1227 NULL, PETSC_TRUE, &dm); 1228 CHKERRQ(ierr); 1229 } 1230 1231 // Distribute the mesh over processes 1232 { 1233 DM dmDist = NULL; 1234 PetscPartitioner part; 1235 1236 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1237 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1238 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1239 if (dmDist) { 1240 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1241 dm = dmDist; 1242 } 1243 } 1244 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1245 1246 // Setup DM 1247 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1248 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1249 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetupData); CHKERRQ(ierr); 1250 1251 // Refine DM for high-order viz 1252 dmviz = NULL; 1253 interpviz = NULL; 1254 if (viz_refine) { 1255 DM dmhierarchy[viz_refine+1]; 1256 1257 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1258 dmhierarchy[0] = dm; 1259 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1260 Mat interp_next; 1261 1262 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1263 CHKERRQ(ierr); 1264 ierr = DMClearDS(dmhierarchy[i+1]); CHKERRQ(ierr); 1265 ierr = DMClearFields(dmhierarchy[i+1]); CHKERRQ(ierr); 1266 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1267 d = (d + 1) / 2; 1268 if (i + 1 == viz_refine) d = 1; 1269 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetupData); 1270 CHKERRQ(ierr); 1271 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1272 &interp_next, NULL); CHKERRQ(ierr); 1273 if (!i) interpviz = interp_next; 1274 else { 1275 Mat C; 1276 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1277 PETSC_DECIDE, &C); CHKERRQ(ierr); 1278 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1279 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1280 interpviz = C; 1281 } 1282 } 1283 for (PetscInt i=1; i<viz_refine; i++) { 1284 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1285 } 1286 dmviz = dmhierarchy[viz_refine]; 1287 } 1288 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1289 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1290 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1291 lnodes /= ncompq; 1292 1293 // Initialize CEED 1294 CeedInit(ceedresource, &ceed); 1295 // Set memtype 1296 CeedMemType memtypebackend; 1297 CeedGetPreferredMemType(ceed, &memtypebackend); 1298 // Check memtype compatibility 1299 if (!setmemtyperequest) 1300 memtyperequested = memtypebackend; 1301 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1302 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1303 "PETSc was not built with CUDA. " 1304 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1305 1306 // Set number of 1D nodes and quadrature points 1307 numP = degree + 1; 1308 numQ = numP + qextra; 1309 1310 // Print summary 1311 if (testChoice == TEST_NONE) { 1312 CeedInt gdofs, odofs; 1313 int comm_size; 1314 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1315 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1316 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1317 gnodes = gdofs/ncompq; 1318 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1319 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1320 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1321 const char *usedresource; 1322 CeedGetResource(ceed, &usedresource); 1323 1324 ierr = PetscPrintf(comm, 1325 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1326 " rank(s) : %d\n" 1327 " Problem:\n" 1328 " Problem Name : %s\n" 1329 " Stabilization : %s\n" 1330 " PETSc:\n" 1331 " Box Faces : %s\n" 1332 " libCEED:\n" 1333 " libCEED Backend : %s\n" 1334 " libCEED Backend MemType : %s\n" 1335 " libCEED User Requested MemType : %s\n" 1336 " Mesh:\n" 1337 " Number of 1D Basis Nodes (P) : %d\n" 1338 " Number of 1D Quadrature Points (Q) : %d\n" 1339 " Global DoFs : %D\n" 1340 " Owned DoFs : %D\n" 1341 " DoFs per node : %D\n" 1342 " Global nodes : %D\n" 1343 " Owned nodes : %D\n", 1344 comm_size, problemTypes[problemChoice], 1345 StabilizationTypes[stab], box_faces_str, usedresource, 1346 CeedMemTypes[memtypebackend], 1347 (setmemtyperequest) ? 1348 CeedMemTypes[memtyperequested] : "none", 1349 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1350 CHKERRQ(ierr); 1351 } 1352 1353 // Set up global mass vector 1354 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1355 1356 // Set up libCEED 1357 // CEED Bases 1358 CeedInit(ceedresource, &ceed); 1359 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1360 &basisq); 1361 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1362 &basisx); 1363 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1364 CEED_GAUSS_LOBATTO, &basisxc); 1365 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1366 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1367 CHKERRQ(ierr); 1368 1369 // CEED Restrictions 1370 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1371 qdatasizeVol, &restrictq, &restrictx, 1372 &restrictqdi); CHKERRQ(ierr); 1373 1374 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1375 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1376 1377 // Create the CEED vectors that will be needed in setup 1378 CeedInt NqptsVol; 1379 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1380 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1381 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1382 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1383 1384 // Create the Q-function that builds the quadrature data for the NS operator 1385 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1386 &qf_setupVol); 1387 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1388 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1389 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1390 1391 // Create the Q-function that sets the ICs of the operator 1392 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1393 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1394 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1395 1396 qf_rhsVol = NULL; 1397 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1398 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1399 problem->applyVol_rhs_loc, &qf_rhsVol); 1400 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1401 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1402 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1403 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1404 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1405 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1406 } 1407 1408 qf_ifunctionVol = NULL; 1409 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1410 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1411 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1412 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1413 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1414 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1415 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1416 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1417 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1418 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1419 } 1420 1421 // Create the operator that builds the quadrature data for the NS operator 1422 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1423 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1424 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1425 basisx, CEED_VECTOR_NONE); 1426 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1427 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1428 1429 // Create the operator that sets the ICs 1430 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1431 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1432 CeedOperatorSetField(op_ics, "q0", restrictq, 1433 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1434 1435 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1436 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1437 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1438 1439 if (qf_rhsVol) { // Create the RHS physics operator 1440 CeedOperator op; 1441 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1442 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1443 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1444 CeedOperatorSetField(op, "qdata", restrictqdi, 1445 CEED_BASIS_COLLOCATED, qdata); 1446 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1447 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1448 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1449 user->op_rhs_vol = op; 1450 } 1451 1452 if (qf_ifunctionVol) { // Create the IFunction operator 1453 CeedOperator op; 1454 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1455 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1456 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1457 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1458 CeedOperatorSetField(op, "qdata", restrictqdi, 1459 CEED_BASIS_COLLOCATED, qdata); 1460 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1461 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1462 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1463 user->op_ifunction_vol = op; 1464 } 1465 1466 // Set up CEED for the boundaries 1467 CeedInt height = 1; 1468 CeedInt dimSur = dim - height; 1469 CeedInt numP_Sur = degree + 1; 1470 CeedInt numQ_Sur = numP_Sur + qextraSur; 1471 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1472 CeedBasis basisxSur, basisxcSur, basisqSur; 1473 CeedInt NqptsSur; 1474 CeedQFunction qf_setupSur, qf_applySur; 1475 1476 // CEED bases for the boundaries 1477 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, 1478 CEED_GAUSS, 1479 &basisqSur); 1480 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1481 &basisxSur); 1482 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1483 CEED_GAUSS_LOBATTO, &basisxcSur); 1484 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1485 1486 // Create the Q-function that builds the quadrature data for the Surface operator 1487 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1488 &qf_setupSur); 1489 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1490 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1491 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1492 1493 // Creat Q-Function for Boundaries 1494 qf_applySur = NULL; 1495 if (problem->applySur) { 1496 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1497 problem->applySur_loc, &qf_applySur); 1498 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1499 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1500 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1501 CeedQFunctionAddOutput(qf_applySur, "v", ncompq, CEED_EVAL_INTERP); 1502 } 1503 1504 // Create CEED Operator for the whole domain 1505 if (!implicit) 1506 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, user->op_rhs_vol, 1507 qf_applySur, qf_setupSur, 1508 height, numP_Sur, numQ_Sur, qdatasizeSur, 1509 NqptsSur, basisxSur, basisqSur, 1510 &user->op_rhs); CHKERRQ(ierr); 1511 if (implicit) 1512 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, 1513 user->op_ifunction_vol, 1514 qf_applySur, qf_setupSur, 1515 height, numP_Sur, numQ_Sur, qdatasizeSur, 1516 NqptsSur, basisxSur, basisqSur, 1517 &user->op_ifunction); CHKERRQ(ierr); 1518 // Set up contex for QFunctions 1519 CeedQFunctionContextCreate(ceed, &ctxSetup); 1520 CeedQFunctionContextSetData(ctxSetup, CEED_MEM_HOST, CEED_USE_POINTER, 1521 sizeof ctxSetupData, &ctxSetupData); 1522 CeedQFunctionSetContext(qf_ics, ctxSetup); 1523 1524 CeedScalar ctxNSData[8] = {lambda, mu, k, cv, cp, g, Rd}; 1525 CeedQFunctionContextCreate(ceed, &ctxNS); 1526 CeedQFunctionContextSetData(ctxNS, CEED_MEM_HOST, CEED_USE_POINTER, 1527 sizeof ctxNSData, &ctxNSData); 1528 1529 struct Advection2dContext_ ctxAdvection2dData = { 1530 .CtauS = CtauS, 1531 .strong_form = strong_form, 1532 .stabilization = stab, 1533 }; 1534 CeedQFunctionContextCreate(ceed, &ctxAdvection2d); 1535 CeedQFunctionContextSetData(ctxAdvection2d, CEED_MEM_HOST, CEED_USE_POINTER, 1536 sizeof ctxAdvection2dData, &ctxAdvection2dData); 1537 1538 struct SurfaceContext_ ctxSurfaceData = { 1539 .E_wind = E_wind, 1540 .strong_form = strong_form, 1541 .implicit = implicit, 1542 }; 1543 CeedQFunctionContextCreate(ceed, &ctxSurface); 1544 CeedQFunctionContextSetData(ctxSurface, CEED_MEM_HOST, CEED_USE_POINTER, 1545 sizeof ctxSurfaceData, &ctxSurfaceData); 1546 1547 switch (problemChoice) { 1548 case NS_DENSITY_CURRENT: 1549 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxNS); 1550 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxNS); 1551 break; 1552 case NS_ADVECTION: 1553 case NS_ADVECTION2D: 1554 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxAdvection2d); 1555 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxAdvection2d); 1556 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, ctxSurface); 1557 } 1558 1559 // Set up PETSc context 1560 // Set up units structure 1561 units->meter = meter; 1562 units->kilogram = kilogram; 1563 units->second = second; 1564 units->Kelvin = Kelvin; 1565 units->Pascal = Pascal; 1566 units->JperkgK = JperkgK; 1567 units->mpersquareds = mpersquareds; 1568 units->WpermK = WpermK; 1569 units->kgpercubicm = kgpercubicm; 1570 units->kgpersquaredms = kgpersquaredms; 1571 units->Joulepercubicm = Joulepercubicm; 1572 units->Joule = Joule; 1573 1574 // Set up user structure 1575 user->comm = comm; 1576 user->outputfreq = outputfreq; 1577 user->contsteps = contsteps; 1578 user->units = units; 1579 user->dm = dm; 1580 user->dmviz = dmviz; 1581 user->interpviz = interpviz; 1582 user->ceed = ceed; 1583 1584 // Calculate qdata and ICs 1585 // Set up state global and local vectors 1586 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1587 1588 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1589 1590 // Apply Setup Ceed Operators 1591 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1592 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1593 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1594 user->M); CHKERRQ(ierr); 1595 1596 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1597 ctxSetup, 0.0); CHKERRQ(ierr); 1598 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1599 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1600 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1601 // slower execution. 1602 Vec Qbc; 1603 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1604 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1605 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1606 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1607 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1608 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1609 ierr = PetscObjectComposeFunction((PetscObject)dm, 1610 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1611 CHKERRQ(ierr); 1612 } 1613 1614 MPI_Comm_rank(comm, &rank); 1615 if (!rank) {ierr = PetscMkdir(user->outputdir); CHKERRQ(ierr);} 1616 // Gather initial Q values 1617 // In case of continuation of simulation, set up initial values from binary file 1618 if (contsteps) { // continue from existent solution 1619 PetscViewer viewer; 1620 char filepath[PETSC_MAX_PATH_LEN]; 1621 // Read input 1622 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1623 user->outputdir); 1624 CHKERRQ(ierr); 1625 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1626 CHKERRQ(ierr); 1627 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1628 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1629 } 1630 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1631 1632 // Create and setup TS 1633 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1634 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1635 if (implicit) { 1636 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1637 if (user->op_ifunction) { 1638 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1639 } else { // Implicit integrators can fall back to using an RHSFunction 1640 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1641 } 1642 } else { 1643 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1644 "Problem does not provide RHSFunction"); 1645 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1646 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1647 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1648 } 1649 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1650 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1651 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1652 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1653 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1654 ierr = TSAdaptSetStepLimits(adapt, 1655 1.e-12 * units->second, 1656 1.e2 * units->second); CHKERRQ(ierr); 1657 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1658 if (!contsteps) { // print initial condition 1659 if (testChoice == TEST_NONE) { 1660 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1661 } 1662 } else { // continue from time of last output 1663 PetscReal time; 1664 PetscInt count; 1665 PetscViewer viewer; 1666 char filepath[PETSC_MAX_PATH_LEN]; 1667 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1668 user->outputdir); CHKERRQ(ierr); 1669 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1670 CHKERRQ(ierr); 1671 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1672 CHKERRQ(ierr); 1673 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1674 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1675 } 1676 if (testChoice == TEST_NONE) { 1677 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1678 } 1679 1680 // Solve 1681 start = MPI_Wtime(); 1682 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1683 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1684 cpu_time_used = MPI_Wtime() - start; 1685 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1686 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1687 comm); CHKERRQ(ierr); 1688 if (testChoice == TEST_NONE) { 1689 ierr = PetscPrintf(PETSC_COMM_WORLD, 1690 "Time taken for solution (sec): %g\n", 1691 (double)cpu_time_used); CHKERRQ(ierr); 1692 } 1693 1694 // Get error 1695 if (problem->non_zero_time && testChoice == TEST_NONE) { 1696 Vec Qexact, Qexactloc; 1697 PetscReal norm; 1698 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1699 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1700 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1701 1702 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1703 restrictq, ctxSetup, ftime); CHKERRQ(ierr); 1704 1705 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1706 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1707 CeedVectorDestroy(&q0ceed); 1708 ierr = PetscPrintf(PETSC_COMM_WORLD, 1709 "Max Error: %g\n", 1710 (double)norm); CHKERRQ(ierr); 1711 // Clean up vectors 1712 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1713 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1714 } 1715 1716 // Output Statistics 1717 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1718 if (testChoice == TEST_NONE) { 1719 ierr = PetscPrintf(PETSC_COMM_WORLD, 1720 "Time integrator took %D time steps to reach final time %g\n", 1721 steps, (double)ftime); CHKERRQ(ierr); 1722 } 1723 // Output numerical values from command line 1724 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1725 1726 // Compare reference solution values with current test run for CI 1727 if (testChoice != TEST_NONE) { 1728 PetscViewer viewer; 1729 // Read reference file 1730 Vec Qref; 1731 PetscReal error, Qrefnorm; 1732 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1733 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1734 CHKERRQ(ierr); 1735 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1736 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1737 1738 // Compute error with respect to reference solution 1739 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1740 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1741 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1742 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1743 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1744 // Check error 1745 if (error > test->testtol) { 1746 ierr = PetscPrintf(PETSC_COMM_WORLD, 1747 "Test failed with error norm %g\n", 1748 (double)error); CHKERRQ(ierr); 1749 } 1750 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1751 } 1752 1753 // Clean up libCEED 1754 CeedVectorDestroy(&qdata); 1755 CeedVectorDestroy(&user->qceed); 1756 CeedVectorDestroy(&user->qdotceed); 1757 CeedVectorDestroy(&user->gceed); 1758 CeedVectorDestroy(&xcorners); 1759 CeedBasisDestroy(&basisq); 1760 CeedBasisDestroy(&basisx); 1761 CeedBasisDestroy(&basisxc); 1762 CeedElemRestrictionDestroy(&restrictq); 1763 CeedElemRestrictionDestroy(&restrictx); 1764 CeedElemRestrictionDestroy(&restrictqdi); 1765 CeedQFunctionDestroy(&qf_setupVol); 1766 CeedQFunctionDestroy(&qf_ics); 1767 CeedQFunctionDestroy(&qf_rhsVol); 1768 CeedQFunctionDestroy(&qf_ifunctionVol); 1769 CeedQFunctionContextDestroy(&ctxSetup); 1770 CeedQFunctionContextDestroy(&ctxNS); 1771 CeedQFunctionContextDestroy(&ctxAdvection2d); 1772 CeedQFunctionContextDestroy(&ctxSurface); 1773 CeedOperatorDestroy(&op_setupVol); 1774 CeedOperatorDestroy(&op_ics); 1775 CeedOperatorDestroy(&user->op_rhs_vol); 1776 CeedOperatorDestroy(&user->op_ifunction_vol); 1777 CeedDestroy(&ceed); 1778 CeedBasisDestroy(&basisqSur); 1779 CeedBasisDestroy(&basisxSur); 1780 CeedBasisDestroy(&basisxcSur); 1781 CeedQFunctionDestroy(&qf_setupSur); 1782 CeedQFunctionDestroy(&qf_applySur); 1783 CeedOperatorDestroy(&user->op_rhs); 1784 CeedOperatorDestroy(&user->op_ifunction); 1785 1786 // Clean up PETSc 1787 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1788 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1789 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1790 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1791 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1792 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1793 ierr = PetscFree(units); CHKERRQ(ierr); 1794 ierr = PetscFree(user); CHKERRQ(ierr); 1795 return PetscFinalize(); 1796 } 1797