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