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/occa -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 <petscts.h> 44 #include <petscdmplex.h> 45 #include <ceed.h> 46 #include <stdbool.h> 47 #include <petscsys.h> 48 #include "common.h" 49 #include "setup-boundary.h" 50 #include "advection.h" 51 #include "advection2d.h" 52 #include "densitycurrent.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 outputfolder[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->outputfolder, 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->outputfolder, 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->outputfolder); 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->outputfolder); 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, SetupContext ctxSetup, CeedScalar time) { 711 PetscErrorCode ierr; 712 CeedVector multlvec; 713 Vec Multiplicity, MultiplicityLoc; 714 715 ctxSetup->time = time; 716 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 717 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 718 CeedOperatorApply(op_ics, xcorners, q0ceed, CEED_REQUEST_IMMEDIATE); 719 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 720 ierr = DMLocalToGlobal(dm, Qloc, ADD_VALUES, Q); CHKERRQ(ierr); 721 722 // Fix multiplicity for output of ICs 723 ierr = DMGetLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 724 CeedElemRestrictionCreateVector(restrictq, &multlvec, NULL); 725 ierr = VectorPlacePetscVec(multlvec, MultiplicityLoc); CHKERRQ(ierr); 726 CeedElemRestrictionGetMultiplicity(restrictq, multlvec); 727 CeedVectorDestroy(&multlvec); 728 ierr = DMGetGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 729 ierr = VecZeroEntries(Multiplicity); CHKERRQ(ierr); 730 ierr = DMLocalToGlobal(dm, MultiplicityLoc, ADD_VALUES, Multiplicity); 731 CHKERRQ(ierr); 732 ierr = VecPointwiseDivide(Q, Q, Multiplicity); CHKERRQ(ierr); 733 ierr = VecPointwiseDivide(Qloc, Qloc, MultiplicityLoc); CHKERRQ(ierr); 734 ierr = DMRestoreLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 735 ierr = DMRestoreGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 736 737 PetscFunctionReturn(0); 738 } 739 740 static PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, 741 CeedElemRestriction restrictq, CeedBasis basisq, 742 CeedElemRestriction restrictqdi, CeedVector qdata, Vec M) { 743 PetscErrorCode ierr; 744 CeedQFunction qf_mass; 745 CeedOperator op_mass; 746 CeedVector mceed; 747 Vec Mloc; 748 CeedInt ncompq, qdatasize; 749 750 PetscFunctionBeginUser; 751 CeedElemRestrictionGetNumComponents(restrictq, &ncompq); 752 CeedElemRestrictionGetNumComponents(restrictqdi, &qdatasize); 753 // Create the Q-function that defines the action of the mass operator 754 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 755 CeedQFunctionAddInput(qf_mass, "q", ncompq, CEED_EVAL_INTERP); 756 CeedQFunctionAddInput(qf_mass, "qdata", qdatasize, CEED_EVAL_NONE); 757 CeedQFunctionAddOutput(qf_mass, "v", ncompq, CEED_EVAL_INTERP); 758 759 // Create the mass operator 760 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 761 CeedOperatorSetField(op_mass, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 762 CeedOperatorSetField(op_mass, "qdata", restrictqdi, 763 CEED_BASIS_COLLOCATED, qdata); 764 CeedOperatorSetField(op_mass, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 765 766 ierr = DMGetLocalVector(dm, &Mloc); CHKERRQ(ierr); 767 ierr = VecZeroEntries(Mloc); CHKERRQ(ierr); 768 CeedElemRestrictionCreateVector(restrictq, &mceed, NULL); 769 ierr = VectorPlacePetscVec(mceed, Mloc); CHKERRQ(ierr); 770 771 { 772 // Compute a lumped mass matrix 773 CeedVector onesvec; 774 CeedElemRestrictionCreateVector(restrictq, &onesvec, NULL); 775 CeedVectorSetValue(onesvec, 1.0); 776 CeedOperatorApply(op_mass, onesvec, mceed, CEED_REQUEST_IMMEDIATE); 777 CeedVectorDestroy(&onesvec); 778 CeedOperatorDestroy(&op_mass); 779 CeedVectorDestroy(&mceed); 780 } 781 CeedQFunctionDestroy(&qf_mass); 782 783 ierr = VecZeroEntries(M); CHKERRQ(ierr); 784 ierr = DMLocalToGlobal(dm, Mloc, ADD_VALUES, M); CHKERRQ(ierr); 785 ierr = DMRestoreLocalVector(dm, &Mloc); CHKERRQ(ierr); 786 787 // Invert diagonally lumped mass vector for RHS function 788 ierr = VecReciprocal(M); CHKERRQ(ierr); 789 PetscFunctionReturn(0); 790 } 791 792 static PetscErrorCode SetUpDM(DM dm, problemData *problem, PetscInt degree, 793 SimpleBC bc, void *ctxSetup) { 794 PetscErrorCode ierr; 795 796 PetscFunctionBeginUser; 797 { 798 // Configure the finite element space and boundary conditions 799 PetscFE fe; 800 PetscInt ncompq = 5; 801 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, ncompq, 802 PETSC_FALSE, degree, PETSC_DECIDE, 803 &fe); CHKERRQ(ierr); 804 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 805 ierr = DMAddField(dm, NULL,(PetscObject)fe); CHKERRQ(ierr); 806 ierr = DMCreateDS(dm); CHKERRQ(ierr); 807 { 808 PetscInt comps[1] = {1}; 809 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 810 1, comps, (void(*)(void))NULL, NULL, bc->nslip[0], 811 bc->slips[0], ctxSetup); CHKERRQ(ierr); 812 comps[0] = 2; 813 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 814 1, comps, (void(*)(void))NULL, NULL, bc->nslip[1], 815 bc->slips[1], ctxSetup); CHKERRQ(ierr); 816 comps[0] = 3; 817 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 818 1, comps, (void(*)(void))NULL, NULL, bc->nslip[2], 819 bc->slips[2], ctxSetup); CHKERRQ(ierr); 820 } 821 if (bc->userbc == PETSC_TRUE) { 822 for (PetscInt c = 0; c < 3; c++) { 823 for (PetscInt s = 0; s < bc->nslip[c]; s++) { 824 for (PetscInt w = 0; w < bc->nwall; w++) { 825 if (bc->slips[c][s] == bc->walls[w]) 826 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, 827 "Boundary condition already set on face %D!\n", 828 bc->walls[w]); 829 830 } 831 } 832 } 833 } 834 // Wall boundary conditions are zero energy density and zero flux for 835 // velocity in advection/advection2d, and zero velocity and zero flux 836 // for mass density and energy density in density_current 837 { 838 if (problem->bc == Exact_Advection || problem->bc == Exact_Advection2d) { 839 PetscInt comps[1] = {4}; 840 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 841 1, comps, (void(*)(void))problem->bc, NULL, 842 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 843 } else if (problem->bc == Exact_DC) { 844 PetscInt comps[3] = {1, 2, 3}; 845 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 846 3, comps, (void(*)(void))problem->bc, NULL, 847 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 848 } else 849 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, 850 "Undefined boundary conditions for this problem"); 851 } 852 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 853 CHKERRQ(ierr); 854 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 855 } 856 { 857 // Empty name for conserved field (because there is only one field) 858 PetscSection section; 859 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 860 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 861 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 862 CHKERRQ(ierr); 863 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 864 CHKERRQ(ierr); 865 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 866 CHKERRQ(ierr); 867 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 868 CHKERRQ(ierr); 869 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 870 CHKERRQ(ierr); 871 } 872 PetscFunctionReturn(0); 873 } 874 875 int main(int argc, char **argv) { 876 PetscInt ierr; 877 MPI_Comm comm; 878 DM dm, dmcoord, dmviz; 879 Mat interpviz; 880 TS ts; 881 TSAdapt adapt; 882 User user; 883 Units units; 884 char ceedresource[4096] = "/cpu/self"; 885 PetscInt localNelemVol, lnodes, gnodes, steps; 886 const PetscInt ncompq = 5; 887 PetscMPIInt rank; 888 PetscScalar ftime; 889 Vec Q, Qloc, Xloc; 890 Ceed ceed; 891 CeedInt numP, numQ; 892 CeedVector xcorners, qdata, q0ceed; 893 CeedBasis basisx, basisxc, basisq; 894 CeedElemRestriction restrictx, restrictq, restrictqdi; 895 CeedQFunction qf_setupVol, qf_ics, qf_rhsVol, qf_ifunctionVol; 896 CeedOperator op_setupVol, op_ics; 897 CeedScalar Rd; 898 CeedMemType memtyperequested; 899 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 900 kgpersquaredms, Joulepercubicm, Joule; 901 problemType problemChoice; 902 problemData *problem = NULL; 903 WindType wind_type; 904 StabilizationType stab; 905 testType testChoice; 906 testData *test = NULL; 907 PetscBool implicit; 908 PetscInt viz_refine = 0; 909 struct SimpleBC_ bc = { 910 .nslip = {2, 2, 2}, 911 .slips = {{5, 6}, {3, 4}, {1, 2}} 912 }; 913 double start, cpu_time_used; 914 // Check PETSc CUDA support 915 PetscBool petschavecuda, setmemtyperequest = PETSC_FALSE; 916 // *INDENT-OFF* 917 #ifdef PETSC_HAVE_CUDA 918 petschavecuda = PETSC_TRUE; 919 #else 920 petschavecuda = PETSC_FALSE; 921 #endif 922 // *INDENT-ON* 923 924 // Create the libCEED contexts 925 PetscScalar meter = 1e-2; // 1 meter in scaled length units 926 PetscScalar second = 1e-2; // 1 second in scaled time units 927 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 928 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 929 CeedScalar theta0 = 300.; // K 930 CeedScalar thetaC = -15.; // K 931 CeedScalar P0 = 1.e5; // Pa 932 CeedScalar E_wind = 1.e6; // J 933 CeedScalar N = 0.01; // 1/s 934 CeedScalar cv = 717.; // J/(kg K) 935 CeedScalar cp = 1004.; // J/(kg K) 936 CeedScalar g = 9.81; // m/s^2 937 CeedScalar lambda = -2./3.; // - 938 CeedScalar mu = 75.; // Pa s, dynamic viscosity 939 // mu = 75 is not physical for air, but is good for numerical stability 940 CeedScalar k = 0.02638; // W/(m K) 941 CeedScalar CtauS = 0.; // dimensionless 942 CeedScalar strong_form = 0.; // [0,1] 943 PetscScalar lx = 8000.; // m 944 PetscScalar ly = 8000.; // m 945 PetscScalar lz = 4000.; // m 946 CeedScalar rc = 1000.; // m (Radius of bubble) 947 PetscScalar resx = 1000.; // m (resolution in x) 948 PetscScalar resy = 1000.; // m (resolution in y) 949 PetscScalar resz = 1000.; // m (resolution in z) 950 PetscInt outputfreq = 10; // - 951 PetscInt contsteps = 0; // - 952 PetscInt degree = 1; // - 953 PetscInt qextra = 2; // - 954 PetscInt qextraSur = 2; // - 955 PetscReal center[3], dc_axis[3] = {0, 0, 0}, wind[3] = {1., 0, 0}; 956 957 ierr = PetscInitialize(&argc, &argv, NULL, help); 958 if (ierr) return ierr; 959 960 // Allocate PETSc context 961 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 962 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 963 964 // Parse command line options 965 comm = PETSC_COMM_WORLD; 966 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 967 NULL); CHKERRQ(ierr); 968 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 969 NULL, ceedresource, ceedresource, 970 sizeof(ceedresource), NULL); CHKERRQ(ierr); 971 testChoice = TEST_NONE; 972 ierr = PetscOptionsEnum("-test", "Run tests", NULL, 973 testTypes, (PetscEnum)testChoice, 974 (PetscEnum *)&testChoice, 975 NULL); CHKERRQ(ierr); 976 test = &testOptions[testChoice]; 977 problemChoice = NS_DENSITY_CURRENT; 978 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 979 problemTypes, (PetscEnum)problemChoice, 980 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 981 problem = &problemOptions[problemChoice]; 982 ierr = PetscOptionsEnum("-problem_advection_wind", "Wind type in Advection", 983 NULL, WindTypes, 984 (PetscEnum)(wind_type = ADVECTION_WIND_ROTATION), 985 (PetscEnum *)&wind_type, NULL); CHKERRQ(ierr); 986 PetscInt n = problem->dim; 987 PetscBool userWind; 988 ierr = PetscOptionsRealArray("-problem_advection_wind_translation", 989 "Constant wind vector", 990 NULL, wind, &n, &userWind); CHKERRQ(ierr); 991 if (wind_type == ADVECTION_WIND_ROTATION && userWind) { 992 ierr = PetscPrintf(comm, 993 "Warning! Use -problem_advection_wind_translation only with -problem_advection_wind translation\n"); 994 CHKERRQ(ierr); 995 } 996 if (wind_type == ADVECTION_WIND_TRANSLATION 997 && problemChoice == NS_DENSITY_CURRENT) { 998 SETERRQ(comm, PETSC_ERR_ARG_INCOMP, 999 "-problem_advection_wind translation is not defined for -problem density_current"); 1000 } 1001 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 1002 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 1003 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 1004 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 1005 NULL, implicit=PETSC_FALSE, &implicit, NULL); 1006 CHKERRQ(ierr); 1007 if (!implicit && stab != STAB_NONE) { 1008 ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); 1009 CHKERRQ(ierr); 1010 } 1011 { 1012 PetscInt len; 1013 PetscBool flg; 1014 ierr = PetscOptionsIntArray("-bc_wall", 1015 "Use wall boundary conditions on this list of faces", 1016 NULL, bc.walls, 1017 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 1018 &len), &flg); CHKERRQ(ierr); 1019 if (flg) { 1020 bc.nwall = len; 1021 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 1022 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1023 } 1024 for (PetscInt j=0; j<3; j++) { 1025 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 1026 ierr = PetscOptionsIntArray(flags[j], 1027 "Use slip boundary conditions on this list of faces", 1028 NULL, bc.slips[j], 1029 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 1030 &len), &flg); 1031 CHKERRQ(ierr); 1032 if (flg) { 1033 bc.nslip[j] = len; 1034 bc.userbc = PETSC_TRUE; 1035 } 1036 } 1037 } 1038 ierr = PetscOptionsInt("-viz_refine", 1039 "Regular refinement levels for visualization", 1040 NULL, viz_refine, &viz_refine, NULL); 1041 CHKERRQ(ierr); 1042 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1043 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1044 meter = fabs(meter); 1045 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1046 NULL, second, &second, NULL); CHKERRQ(ierr); 1047 second = fabs(second); 1048 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1049 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1050 kilogram = fabs(kilogram); 1051 ierr = PetscOptionsScalar("-units_Kelvin", 1052 "1 Kelvin in scaled temperature units", 1053 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1054 Kelvin = fabs(Kelvin); 1055 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1056 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1057 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1058 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1059 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1060 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1061 ierr = PetscOptionsScalar("-E_wind", "Total energy of inflow wind", 1062 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 1063 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1064 NULL, N, &N, NULL); CHKERRQ(ierr); 1065 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1066 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1067 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1068 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1069 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1070 NULL, g, &g, NULL); CHKERRQ(ierr); 1071 ierr = PetscOptionsScalar("-lambda", 1072 "Stokes hypothesis second viscosity coefficient", 1073 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1074 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1075 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1076 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1077 NULL, k, &k, NULL); CHKERRQ(ierr); 1078 ierr = PetscOptionsScalar("-CtauS", 1079 "Scale coefficient for tau (nondimensional)", 1080 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1081 if (stab == STAB_NONE && CtauS != 0) { 1082 ierr = PetscPrintf(comm, 1083 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1084 CHKERRQ(ierr); 1085 } 1086 ierr = PetscOptionsScalar("-strong_form", 1087 "Strong (1) or weak/integrated by parts (0) advection residual", 1088 NULL, strong_form, &strong_form, NULL); 1089 CHKERRQ(ierr); 1090 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1091 ierr = PetscPrintf(comm, 1092 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1093 CHKERRQ(ierr); 1094 } 1095 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1096 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1097 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1098 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1099 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1100 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1101 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1102 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1103 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1104 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1105 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1106 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1107 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1108 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1109 n = problem->dim; 1110 center[0] = 0.5 * lx; 1111 center[1] = 0.5 * ly; 1112 center[2] = 0.5 * lz; 1113 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1114 NULL, center, &n, NULL); CHKERRQ(ierr); 1115 n = problem->dim; 1116 ierr = PetscOptionsRealArray("-dc_axis", 1117 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1118 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1119 { 1120 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1121 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1122 if (norm > 0) { 1123 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1124 } 1125 } 1126 ierr = PetscOptionsInt("-output_freq", 1127 "Frequency of output, in number of steps", 1128 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1129 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1130 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1131 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1132 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1133 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1134 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1135 PetscBool userQextraSur; 1136 ierr = PetscOptionsInt("-qextra_boundary", 1137 "Number of extra quadrature points on in/outflow faces", 1138 NULL, qextraSur, &qextraSur, &userQextraSur); 1139 CHKERRQ(ierr); 1140 if ((wind_type == ADVECTION_WIND_ROTATION 1141 || problemChoice == NS_DENSITY_CURRENT) && userQextraSur) { 1142 ierr = PetscPrintf(comm, 1143 "Warning! Use -qextra_boundary only with -problem_advection_wind translation\n"); 1144 CHKERRQ(ierr); 1145 } 1146 ierr = PetscStrncpy(user->outputfolder, ".", 2); CHKERRQ(ierr); 1147 ierr = PetscOptionsString("-of", "Output folder", 1148 NULL, user->outputfolder, user->outputfolder, 1149 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 1150 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1151 ierr = PetscOptionsEnum("-memtype", 1152 "CEED MemType requested", NULL, 1153 memTypes, (PetscEnum)memtyperequested, 1154 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1155 CHKERRQ(ierr); 1156 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1157 1158 // Define derived units 1159 Pascal = kilogram / (meter * PetscSqr(second)); 1160 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1161 mpersquareds = meter / PetscSqr(second); 1162 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1163 kgpercubicm = kilogram / pow(meter,3); 1164 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1165 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1166 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 1167 1168 // Scale variables to desired units 1169 theta0 *= Kelvin; 1170 thetaC *= Kelvin; 1171 P0 *= Pascal; 1172 E_wind *= Joule; 1173 N *= (1./second); 1174 cv *= JperkgK; 1175 cp *= JperkgK; 1176 Rd = cp - cv; 1177 g *= mpersquareds; 1178 mu *= Pascal * second; 1179 k *= WpermK; 1180 lx = fabs(lx) * meter; 1181 ly = fabs(ly) * meter; 1182 lz = fabs(lz) * meter; 1183 rc = fabs(rc) * meter; 1184 resx = fabs(resx) * meter; 1185 resy = fabs(resy) * meter; 1186 resz = fabs(resz) * meter; 1187 for (int i=0; i<3; i++) center[i] *= meter; 1188 1189 const CeedInt dim = problem->dim, ncompx = problem->dim, 1190 qdatasizeVol = problem->qdatasizeVol; 1191 // Set up the libCEED context 1192 struct SetupContext_ ctxSetup = { 1193 .theta0 = theta0, 1194 .thetaC = thetaC, 1195 .P0 = P0, 1196 .N = N, 1197 .cv = cv, 1198 .cp = cp, 1199 .Rd = Rd, 1200 .g = g, 1201 .rc = rc, 1202 .lx = lx, 1203 .ly = ly, 1204 .lz = lz, 1205 .center[0] = center[0], 1206 .center[1] = center[1], 1207 .center[2] = center[2], 1208 .dc_axis[0] = dc_axis[0], 1209 .dc_axis[1] = dc_axis[1], 1210 .dc_axis[2] = dc_axis[2], 1211 .wind[0] = wind[0], 1212 .wind[1] = wind[1], 1213 .wind[2] = wind[2], 1214 .time = 0, 1215 .wind_type = wind_type, 1216 }; 1217 1218 // Create the mesh 1219 { 1220 const PetscReal scale[3] = {lx, ly, lz}; 1221 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1222 NULL, PETSC_TRUE, &dm); 1223 CHKERRQ(ierr); 1224 } 1225 1226 // Distribute the mesh over processes 1227 { 1228 DM dmDist = NULL; 1229 PetscPartitioner part; 1230 1231 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1232 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1233 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1234 if (dmDist) { 1235 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1236 dm = dmDist; 1237 } 1238 } 1239 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1240 1241 // Setup DM 1242 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1243 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1244 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetup); CHKERRQ(ierr); 1245 1246 // Refine DM for high-order viz 1247 dmviz = NULL; 1248 interpviz = NULL; 1249 if (viz_refine) { 1250 DM dmhierarchy[viz_refine+1]; 1251 1252 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1253 dmhierarchy[0] = dm; 1254 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1255 Mat interp_next; 1256 1257 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1258 CHKERRQ(ierr); 1259 ierr = DMClearDS(dmhierarchy[i+1]); CHKERRQ(ierr); 1260 ierr = DMClearFields(dmhierarchy[i+1]); CHKERRQ(ierr); 1261 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1262 d = (d + 1) / 2; 1263 if (i + 1 == viz_refine) d = 1; 1264 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); 1265 CHKERRQ(ierr); 1266 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1267 &interp_next, NULL); CHKERRQ(ierr); 1268 if (!i) interpviz = interp_next; 1269 else { 1270 Mat C; 1271 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1272 PETSC_DECIDE, &C); CHKERRQ(ierr); 1273 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1274 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1275 interpviz = C; 1276 } 1277 } 1278 for (PetscInt i=1; i<viz_refine; i++) { 1279 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1280 } 1281 dmviz = dmhierarchy[viz_refine]; 1282 } 1283 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1284 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1285 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1286 lnodes /= ncompq; 1287 1288 // Initialize CEED 1289 CeedInit(ceedresource, &ceed); 1290 // Set memtype 1291 CeedMemType memtypebackend; 1292 CeedGetPreferredMemType(ceed, &memtypebackend); 1293 // Check memtype compatibility 1294 if (!setmemtyperequest) 1295 memtyperequested = memtypebackend; 1296 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1297 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1298 "PETSc was not built with CUDA. " 1299 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1300 1301 // Set number of 1D nodes and quadrature points 1302 numP = degree + 1; 1303 numQ = numP + qextra; 1304 1305 // Print summary 1306 if (testChoice == TEST_NONE) { 1307 CeedInt gdofs, odofs; 1308 int comm_size; 1309 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1310 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1311 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1312 gnodes = gdofs/ncompq; 1313 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1314 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1315 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1316 const char *usedresource; 1317 CeedGetResource(ceed, &usedresource); 1318 1319 ierr = PetscPrintf(comm, 1320 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1321 " rank(s) : %d\n" 1322 " Problem:\n" 1323 " Problem Name : %s\n" 1324 " Stabilization : %s\n" 1325 " PETSc:\n" 1326 " Box Faces : %s\n" 1327 " libCEED:\n" 1328 " libCEED Backend : %s\n" 1329 " libCEED Backend MemType : %s\n" 1330 " libCEED User Requested MemType : %s\n" 1331 " Mesh:\n" 1332 " Number of 1D Basis Nodes (P) : %d\n" 1333 " Number of 1D Quadrature Points (Q) : %d\n" 1334 " Global DoFs : %D\n" 1335 " Owned DoFs : %D\n" 1336 " DoFs per node : %D\n" 1337 " Global nodes : %D\n" 1338 " Owned nodes : %D\n", 1339 comm_size, problemTypes[problemChoice], 1340 StabilizationTypes[stab], box_faces_str, usedresource, 1341 CeedMemTypes[memtypebackend], 1342 (setmemtyperequest) ? 1343 CeedMemTypes[memtyperequested] : "none", 1344 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1345 CHKERRQ(ierr); 1346 } 1347 1348 // Set up global mass vector 1349 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1350 1351 // Set up libCEED 1352 // CEED Bases 1353 CeedInit(ceedresource, &ceed); 1354 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1355 &basisq); 1356 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1357 &basisx); 1358 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1359 CEED_GAUSS_LOBATTO, &basisxc); 1360 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1361 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1362 CHKERRQ(ierr); 1363 1364 // CEED Restrictions 1365 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1366 qdatasizeVol, &restrictq, &restrictx, 1367 &restrictqdi); CHKERRQ(ierr); 1368 1369 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1370 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1371 1372 // Create the CEED vectors that will be needed in setup 1373 CeedInt NqptsVol; 1374 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1375 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1376 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1377 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1378 1379 // Create the Q-function that builds the quadrature data for the NS operator 1380 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1381 &qf_setupVol); 1382 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1383 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1384 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1385 1386 // Create the Q-function that sets the ICs of the operator 1387 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1388 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1389 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1390 1391 qf_rhsVol = NULL; 1392 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1393 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1394 problem->applyVol_rhs_loc, &qf_rhsVol); 1395 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1396 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1397 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1398 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1399 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1400 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1401 } 1402 1403 qf_ifunctionVol = NULL; 1404 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1405 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1406 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1407 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1408 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1409 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1410 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1411 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1412 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1413 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1414 } 1415 1416 // Create the operator that builds the quadrature data for the NS operator 1417 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1418 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1419 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1420 basisx, CEED_VECTOR_NONE); 1421 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1422 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1423 1424 // Create the operator that sets the ICs 1425 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1426 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1427 CeedOperatorSetField(op_ics, "q0", restrictq, 1428 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1429 1430 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1431 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1432 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1433 1434 if (qf_rhsVol) { // Create the RHS physics operator 1435 CeedOperator op; 1436 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1437 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1438 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1439 CeedOperatorSetField(op, "qdata", restrictqdi, 1440 CEED_BASIS_COLLOCATED, qdata); 1441 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1442 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1443 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1444 user->op_rhs_vol = op; 1445 } 1446 1447 if (qf_ifunctionVol) { // Create the IFunction operator 1448 CeedOperator op; 1449 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1450 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1451 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1452 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1453 CeedOperatorSetField(op, "qdata", restrictqdi, 1454 CEED_BASIS_COLLOCATED, qdata); 1455 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1456 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1457 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1458 user->op_ifunction_vol = op; 1459 } 1460 1461 // Set up CEED for the boundaries 1462 CeedInt height = 1; 1463 CeedInt dimSur = dim - height; 1464 CeedInt numP_Sur = degree + 1; 1465 CeedInt numQ_Sur = numP_Sur + qextraSur; 1466 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1467 CeedBasis basisxSur, basisxcSur, basisqSur; 1468 CeedInt NqptsSur; 1469 CeedQFunction qf_setupSur, qf_applySur; 1470 1471 // CEED bases for the boundaries 1472 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, 1473 CEED_GAUSS, 1474 &basisqSur); 1475 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1476 &basisxSur); 1477 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1478 CEED_GAUSS_LOBATTO, &basisxcSur); 1479 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1480 1481 // Create the Q-function that builds the quadrature data for the Surface operator 1482 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1483 &qf_setupSur); 1484 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1485 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1486 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1487 1488 // Creat Q-Function for Boundaries 1489 qf_applySur = NULL; 1490 if (problem->applySur) { 1491 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1492 problem->applySur_loc, &qf_applySur); 1493 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1494 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1495 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1496 CeedQFunctionAddOutput(qf_applySur, "v", ncompq, CEED_EVAL_INTERP); 1497 } 1498 1499 // Create CEED Operator for the whole domain 1500 if (!implicit) 1501 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, user->op_rhs_vol, 1502 qf_applySur, qf_setupSur, 1503 height, numP_Sur, numQ_Sur, qdatasizeSur, 1504 NqptsSur, basisxSur, basisqSur, 1505 &user->op_rhs); CHKERRQ(ierr); 1506 if (implicit) 1507 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, 1508 user->op_ifunction_vol, 1509 qf_applySur, qf_setupSur, 1510 height, numP_Sur, numQ_Sur, qdatasizeSur, 1511 NqptsSur, basisxSur, basisqSur, 1512 &user->op_ifunction); CHKERRQ(ierr); 1513 // Set up contex for QFunctions 1514 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1515 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1516 struct Advection2dContext_ ctxAdvection2d = { 1517 .CtauS = CtauS, 1518 .strong_form = strong_form, 1519 .stabilization = stab, 1520 }; 1521 struct SurfaceContext_ ctxSurface = { 1522 .E_wind = E_wind, 1523 .strong_form = strong_form, 1524 .implicit = implicit, 1525 }; 1526 switch (problemChoice) { 1527 case NS_DENSITY_CURRENT: 1528 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1529 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1530 sizeof ctxNS); 1531 break; 1532 case NS_ADVECTION: 1533 case NS_ADVECTION2D: 1534 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1535 sizeof ctxAdvection2d); 1536 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1537 sizeof ctxAdvection2d); 1538 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &ctxSurface, 1539 sizeof ctxSurface); 1540 } 1541 1542 // Set up PETSc context 1543 // Set up units structure 1544 units->meter = meter; 1545 units->kilogram = kilogram; 1546 units->second = second; 1547 units->Kelvin = Kelvin; 1548 units->Pascal = Pascal; 1549 units->JperkgK = JperkgK; 1550 units->mpersquareds = mpersquareds; 1551 units->WpermK = WpermK; 1552 units->kgpercubicm = kgpercubicm; 1553 units->kgpersquaredms = kgpersquaredms; 1554 units->Joulepercubicm = Joulepercubicm; 1555 units->Joule = Joule; 1556 1557 // Set up user structure 1558 user->comm = comm; 1559 user->outputfreq = outputfreq; 1560 user->contsteps = contsteps; 1561 user->units = units; 1562 user->dm = dm; 1563 user->dmviz = dmviz; 1564 user->interpviz = interpviz; 1565 user->ceed = ceed; 1566 1567 // Calculate qdata and ICs 1568 // Set up state global and local vectors 1569 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1570 1571 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1572 1573 // Apply Setup Ceed Operators 1574 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1575 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1576 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1577 user->M); CHKERRQ(ierr); 1578 1579 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1580 &ctxSetup, 0.0); CHKERRQ(ierr); 1581 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1582 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1583 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1584 // slower execution. 1585 Vec Qbc; 1586 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1587 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1588 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1589 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1590 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1591 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1592 ierr = PetscObjectComposeFunction((PetscObject)dm, 1593 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1594 CHKERRQ(ierr); 1595 } 1596 1597 MPI_Comm_rank(comm, &rank); 1598 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1599 // Gather initial Q values 1600 // In case of continuation of simulation, set up initial values from binary file 1601 if (contsteps) { // continue from existent solution 1602 PetscViewer viewer; 1603 char filepath[PETSC_MAX_PATH_LEN]; 1604 // Read input 1605 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1606 user->outputfolder); 1607 CHKERRQ(ierr); 1608 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1609 CHKERRQ(ierr); 1610 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1611 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1612 } 1613 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1614 1615 // Create and setup TS 1616 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1617 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1618 if (implicit) { 1619 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1620 if (user->op_ifunction) { 1621 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1622 } else { // Implicit integrators can fall back to using an RHSFunction 1623 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1624 } 1625 } else { 1626 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1627 "Problem does not provide RHSFunction"); 1628 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1629 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1630 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1631 } 1632 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1633 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1634 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1635 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1636 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1637 ierr = TSAdaptSetStepLimits(adapt, 1638 1.e-12 * units->second, 1639 1.e2 * units->second); CHKERRQ(ierr); 1640 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1641 if (!contsteps) { // print initial condition 1642 if (testChoice == TEST_NONE) { 1643 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1644 } 1645 } else { // continue from time of last output 1646 PetscReal time; 1647 PetscInt count; 1648 PetscViewer viewer; 1649 char filepath[PETSC_MAX_PATH_LEN]; 1650 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1651 user->outputfolder); CHKERRQ(ierr); 1652 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1653 CHKERRQ(ierr); 1654 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1655 CHKERRQ(ierr); 1656 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1657 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1658 } 1659 if (testChoice == TEST_NONE) { 1660 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1661 } 1662 1663 // Solve 1664 start = MPI_Wtime(); 1665 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1666 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1667 cpu_time_used = MPI_Wtime() - start; 1668 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1669 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1670 comm); CHKERRQ(ierr); 1671 if (testChoice == TEST_NONE) { 1672 ierr = PetscPrintf(PETSC_COMM_WORLD, 1673 "Time taken for solution (sec): %g\n", 1674 (double)cpu_time_used); CHKERRQ(ierr); 1675 } 1676 1677 // Get error 1678 if (problem->non_zero_time && testChoice == TEST_NONE) { 1679 Vec Qexact, Qexactloc; 1680 PetscReal norm; 1681 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1682 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1683 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1684 1685 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1686 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1687 1688 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1689 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1690 CeedVectorDestroy(&q0ceed); 1691 ierr = PetscPrintf(PETSC_COMM_WORLD, 1692 "Max Error: %g\n", 1693 (double)norm); CHKERRQ(ierr); 1694 // Clean up vectors 1695 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1696 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1697 } 1698 1699 // Output Statistics 1700 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1701 if (testChoice == TEST_NONE) { 1702 ierr = PetscPrintf(PETSC_COMM_WORLD, 1703 "Time integrator took %D time steps to reach final time %g\n", 1704 steps, (double)ftime); CHKERRQ(ierr); 1705 } 1706 // Output numerical values from command line 1707 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1708 1709 // Compare reference solution values with current test run for CI 1710 if (testChoice != TEST_NONE) { 1711 PetscViewer viewer; 1712 // Read reference file 1713 Vec Qref; 1714 PetscReal error, Qrefnorm; 1715 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1716 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1717 CHKERRQ(ierr); 1718 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1719 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1720 1721 // Compute error with respect to reference solution 1722 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1723 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1724 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1725 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1726 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1727 // Check error 1728 if (error > test->testtol) { 1729 ierr = PetscPrintf(PETSC_COMM_WORLD, 1730 "Test failed with error norm %g\n", 1731 (double)error); CHKERRQ(ierr); 1732 } 1733 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1734 } 1735 1736 // Clean up libCEED 1737 CeedVectorDestroy(&qdata); 1738 CeedVectorDestroy(&user->qceed); 1739 CeedVectorDestroy(&user->qdotceed); 1740 CeedVectorDestroy(&user->gceed); 1741 CeedVectorDestroy(&xcorners); 1742 CeedBasisDestroy(&basisq); 1743 CeedBasisDestroy(&basisx); 1744 CeedBasisDestroy(&basisxc); 1745 CeedElemRestrictionDestroy(&restrictq); 1746 CeedElemRestrictionDestroy(&restrictx); 1747 CeedElemRestrictionDestroy(&restrictqdi); 1748 CeedQFunctionDestroy(&qf_setupVol); 1749 CeedQFunctionDestroy(&qf_ics); 1750 CeedQFunctionDestroy(&qf_rhsVol); 1751 CeedQFunctionDestroy(&qf_ifunctionVol); 1752 CeedOperatorDestroy(&op_setupVol); 1753 CeedOperatorDestroy(&op_ics); 1754 CeedOperatorDestroy(&user->op_rhs_vol); 1755 CeedOperatorDestroy(&user->op_ifunction_vol); 1756 CeedDestroy(&ceed); 1757 CeedBasisDestroy(&basisqSur); 1758 CeedBasisDestroy(&basisxSur); 1759 CeedBasisDestroy(&basisxcSur); 1760 CeedQFunctionDestroy(&qf_setupSur); 1761 CeedQFunctionDestroy(&qf_applySur); 1762 CeedOperatorDestroy(&user->op_rhs); 1763 CeedOperatorDestroy(&user->op_ifunction); 1764 1765 // Clean up PETSc 1766 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1767 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1768 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1769 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1770 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1771 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1772 ierr = PetscFree(units); CHKERRQ(ierr); 1773 ierr = PetscFree(user); CHKERRQ(ierr); 1774 return PetscFinalize(); 1775 } 1776 1777