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