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