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 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 961 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 962 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 963 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 964 NULL, implicit=PETSC_FALSE, &implicit, NULL); 965 CHKERRQ(ierr); 966 if (!implicit && stab != STAB_NONE) { 967 ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); 968 CHKERRQ(ierr); 969 } 970 { 971 PetscInt len; 972 PetscBool flg; 973 ierr = PetscOptionsIntArray("-bc_wall", 974 "Use wall boundary conditions on this list of faces", 975 NULL, bc.walls, 976 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 977 &len), &flg); CHKERRQ(ierr); 978 if (flg) { 979 bc.nwall = len; 980 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 981 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 982 } 983 for (PetscInt j=0; j<3; j++) { 984 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 985 ierr = PetscOptionsIntArray(flags[j], 986 "Use slip boundary conditions on this list of faces", 987 NULL, bc.slips[j], 988 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 989 &len), &flg); 990 CHKERRQ(ierr); 991 if (flg) { 992 bc.nslip[j] = len; 993 bc.userbc = PETSC_TRUE; 994 } 995 } 996 } 997 ierr = PetscOptionsInt("-viz_refine", 998 "Regular refinement levels for visualization", 999 NULL, viz_refine, &viz_refine, NULL); 1000 CHKERRQ(ierr); 1001 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1002 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1003 meter = fabs(meter); 1004 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1005 NULL, second, &second, NULL); CHKERRQ(ierr); 1006 second = fabs(second); 1007 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1008 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1009 kilogram = fabs(kilogram); 1010 ierr = PetscOptionsScalar("-units_Kelvin", 1011 "1 Kelvin in scaled temperature units", 1012 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1013 Kelvin = fabs(Kelvin); 1014 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1015 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1016 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1017 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1018 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1019 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1020 ierr = PetscOptionsScalar("-P_wind", "Inflow wind pressure", 1021 NULL, P_wind, &P_wind, NULL); CHKERRQ(ierr); 1022 ierr = PetscOptionsScalar("-rho_wind", "Inflow wind density", 1023 NULL, rho_wind, &rho_wind, NULL); CHKERRQ(ierr); 1024 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1025 NULL, N, &N, NULL); CHKERRQ(ierr); 1026 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1027 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1028 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1029 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1030 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1031 NULL, g, &g, NULL); CHKERRQ(ierr); 1032 ierr = PetscOptionsScalar("-lambda", 1033 "Stokes hypothesis second viscosity coefficient", 1034 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1035 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1036 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1037 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1038 NULL, k, &k, NULL); CHKERRQ(ierr); 1039 ierr = PetscOptionsScalar("-CtauS", 1040 "Scale coefficient for tau (nondimensional)", 1041 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1042 if (stab == STAB_NONE && CtauS != 0) { 1043 ierr = PetscPrintf(comm, 1044 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1045 CHKERRQ(ierr); 1046 } 1047 ierr = PetscOptionsScalar("-strong_form", 1048 "Strong (1) or weak/integrated by parts (0) advection residual", 1049 NULL, strong_form, &strong_form, NULL); 1050 CHKERRQ(ierr); 1051 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1052 ierr = PetscPrintf(comm, 1053 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1054 CHKERRQ(ierr); 1055 } 1056 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1057 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1058 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1059 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1060 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1061 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1062 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1063 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1064 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1065 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1066 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1067 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1068 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1069 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1070 n = problem->dim; 1071 center[0] = 0.5 * lx; 1072 center[1] = 0.5 * ly; 1073 center[2] = 0.5 * lz; 1074 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1075 NULL, center, &n, NULL); CHKERRQ(ierr); 1076 n = problem->dim; 1077 ierr = PetscOptionsRealArray("-dc_axis", 1078 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1079 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1080 { 1081 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1082 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1083 if (norm > 0) { 1084 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1085 } 1086 } 1087 ierr = PetscOptionsInt("-output_freq", 1088 "Frequency of output, in number of steps", 1089 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1090 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1091 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1092 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1093 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1094 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1095 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1096 PetscBool userQextraSur; 1097 ierr = PetscOptionsInt("-qextra_boundary", "Number of extra quadrature points on in/outflow faces", 1098 NULL, qextraSur, &qextraSur, &userQextraSur); CHKERRQ(ierr); 1099 if ((wind_type == ADVECTION_WIND_ROTATION || problemChoice == NS_DENSITY_CURRENT) && userQextraSur) { 1100 ierr = PetscPrintf(comm, "Warning! Use -qextra_boundary only with -problem_advection_wind translation\n"); 1101 CHKERRQ(ierr); 1102 } 1103 ierr = PetscStrncpy(user->outputfolder, ".", 2); CHKERRQ(ierr); 1104 ierr = PetscOptionsString("-of", "Output folder", 1105 NULL, user->outputfolder, user->outputfolder, 1106 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 1107 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1108 ierr = PetscOptionsEnum("-memtype", 1109 "CEED MemType requested", NULL, 1110 memTypes, (PetscEnum)memtyperequested, 1111 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1112 CHKERRQ(ierr); 1113 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1114 1115 // Define derived units 1116 Pascal = kilogram / (meter * PetscSqr(second)); 1117 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1118 mpersquareds = meter / PetscSqr(second); 1119 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1120 kgpercubicm = kilogram / pow(meter,3); 1121 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1122 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1123 1124 // Scale variables to desired units 1125 theta0 *= Kelvin; 1126 thetaC *= Kelvin; 1127 P0 *= Pascal; 1128 P_wind *= Pascal; 1129 rho_wind *= kgpercubicm; 1130 N *= (1./second); 1131 cv *= JperkgK; 1132 cp *= JperkgK; 1133 Rd = cp - cv; 1134 g *= mpersquareds; 1135 mu *= Pascal * second; 1136 k *= WpermK; 1137 lx = fabs(lx) * meter; 1138 ly = fabs(ly) * meter; 1139 lz = fabs(lz) * meter; 1140 rc = fabs(rc) * meter; 1141 resx = fabs(resx) * meter; 1142 resy = fabs(resy) * meter; 1143 resz = fabs(resz) * meter; 1144 for (int i=0; i<3; i++) center[i] *= meter; 1145 1146 const CeedInt dim = problem->dim, ncompx = problem->dim, 1147 qdatasizeVol = problem->qdatasizeVol; 1148 // Set up the libCEED context 1149 struct SetupContext_ ctxSetup = { 1150 .theta0 = theta0, 1151 .thetaC = thetaC, 1152 .P0 = P0, 1153 .N = N, 1154 .cv = cv, 1155 .cp = cp, 1156 .Rd = Rd, 1157 .g = g, 1158 .rc = rc, 1159 .lx = lx, 1160 .ly = ly, 1161 .lz = lz, 1162 .center[0] = center[0], 1163 .center[1] = center[1], 1164 .center[2] = center[2], 1165 .dc_axis[0] = dc_axis[0], 1166 .dc_axis[1] = dc_axis[1], 1167 .dc_axis[2] = dc_axis[2], 1168 .wind[0] = wind[0], 1169 .wind[1] = wind[1], 1170 .wind[2] = wind[2], 1171 .time = 0, 1172 .wind_type = wind_type, 1173 }; 1174 1175 // Create the mesh 1176 { 1177 const PetscReal scale[3] = {lx, ly, lz}; 1178 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1179 NULL, PETSC_TRUE, &dm); 1180 CHKERRQ(ierr); 1181 } 1182 1183 // Distribute the mesh over processes 1184 { 1185 DM dmDist = NULL; 1186 PetscPartitioner part; 1187 1188 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1189 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1190 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1191 if (dmDist) { 1192 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1193 dm = dmDist; 1194 } 1195 } 1196 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1197 1198 // Setup DM 1199 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1200 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1201 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetup); CHKERRQ(ierr); 1202 1203 // Refine DM for high-order viz 1204 dmviz = NULL; 1205 interpviz = NULL; 1206 if (viz_refine) { 1207 DM dmhierarchy[viz_refine+1]; 1208 1209 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1210 dmhierarchy[0] = dm; 1211 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1212 Mat interp_next; 1213 1214 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1215 CHKERRQ(ierr); 1216 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1217 d = (d + 1) / 2; 1218 if (i + 1 == viz_refine) d = 1; 1219 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 1220 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1221 &interp_next, NULL); CHKERRQ(ierr); 1222 if (!i) interpviz = interp_next; 1223 else { 1224 Mat C; 1225 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1226 PETSC_DECIDE, &C); CHKERRQ(ierr); 1227 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1228 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1229 interpviz = C; 1230 } 1231 } 1232 for (PetscInt i=1; i<viz_refine; i++) { 1233 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1234 } 1235 dmviz = dmhierarchy[viz_refine]; 1236 } 1237 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1238 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1239 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1240 lnodes /= ncompq; 1241 1242 // Initialize CEED 1243 CeedInit(ceedresource, &ceed); 1244 // Set memtype 1245 CeedMemType memtypebackend; 1246 CeedGetPreferredMemType(ceed, &memtypebackend); 1247 // Check memtype compatibility 1248 if (!setmemtyperequest) 1249 memtyperequested = memtypebackend; 1250 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1251 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1252 "PETSc was not built with CUDA. " 1253 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1254 1255 // Set number of 1D nodes and quadrature points 1256 numP = degree + 1; 1257 numQ = numP + qextra; 1258 1259 // Print summary 1260 if (testChoice == TEST_NONE) { 1261 CeedInt gdofs, odofs; 1262 int comm_size; 1263 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1264 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1265 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1266 gnodes = gdofs/ncompq; 1267 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1268 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1269 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1270 const char *usedresource; 1271 CeedGetResource(ceed, &usedresource); 1272 1273 ierr = PetscPrintf(comm, 1274 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1275 " rank(s) : %d\n" 1276 " Problem:\n" 1277 " Problem Name : %s\n" 1278 " Stabilization : %s\n" 1279 " PETSc:\n" 1280 " Box Faces : %s\n" 1281 " libCEED:\n" 1282 " libCEED Backend : %s\n" 1283 " libCEED Backend MemType : %s\n" 1284 " libCEED User Requested MemType : %s\n" 1285 " Mesh:\n" 1286 " Number of 1D Basis Nodes (P) : %d\n" 1287 " Number of 1D Quadrature Points (Q) : %d\n" 1288 " Global DoFs : %D\n" 1289 " Owned DoFs : %D\n" 1290 " DoFs per node : %D\n" 1291 " Global nodes : %D\n" 1292 " Owned nodes : %D\n", 1293 comm_size, problemTypes[problemChoice], 1294 StabilizationTypes[stab], box_faces_str, usedresource, 1295 CeedMemTypes[memtypebackend], 1296 (setmemtyperequest) ? 1297 CeedMemTypes[memtyperequested] : "none", 1298 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1299 CHKERRQ(ierr); 1300 } 1301 1302 // Set up global mass vector 1303 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1304 1305 // Set up libCEED 1306 // CEED Bases 1307 CeedInit(ceedresource, &ceed); 1308 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1309 &basisq); 1310 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1311 &basisx); 1312 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1313 CEED_GAUSS_LOBATTO, &basisxc); 1314 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1315 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1316 CHKERRQ(ierr); 1317 1318 // CEED Restrictions 1319 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1320 qdatasizeVol, &restrictq, &restrictx, 1321 &restrictqdi); CHKERRQ(ierr); 1322 1323 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1324 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1325 1326 // Create the CEED vectors that will be needed in setup 1327 CeedInt NqptsVol; 1328 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1329 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1330 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1331 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1332 1333 // Create the Q-function that builds the quadrature data for the NS operator 1334 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1335 &qf_setupVol); 1336 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1337 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1338 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1339 1340 // Create the Q-function that sets the ICs of the operator 1341 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1342 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1343 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1344 1345 qf_rhsVol = NULL; 1346 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1347 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1348 problem->applyVol_rhs_loc, &qf_rhsVol); 1349 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1350 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1351 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1352 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1353 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1354 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1355 } 1356 1357 qf_ifunctionVol = NULL; 1358 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1359 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1360 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1361 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1362 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1363 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1364 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1365 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1366 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1367 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1368 } 1369 1370 // Create the operator that builds the quadrature data for the NS operator 1371 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1372 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1373 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1374 basisx, CEED_VECTOR_NONE); 1375 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1376 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1377 1378 // Create the operator that sets the ICs 1379 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1380 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1381 CeedOperatorSetField(op_ics, "q0", restrictq, 1382 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1383 1384 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1385 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1386 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1387 1388 if (qf_rhsVol) { // Create the RHS physics operator 1389 CeedOperator op; 1390 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1391 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1392 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1393 CeedOperatorSetField(op, "qdata", restrictqdi, 1394 CEED_BASIS_COLLOCATED, qdata); 1395 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1396 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1397 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1398 user->op_rhs_vol = op; 1399 } 1400 1401 if (qf_ifunctionVol) { // Create the IFunction operator 1402 CeedOperator op; 1403 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1404 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1405 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1406 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1407 CeedOperatorSetField(op, "qdata", restrictqdi, 1408 CEED_BASIS_COLLOCATED, qdata); 1409 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1410 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1411 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1412 user->op_ifunction_vol = op; 1413 } 1414 1415 // Set up CEED for the boundaries 1416 CeedInt height = 1; 1417 CeedInt dimSur = dim - height; 1418 CeedInt numP_Sur = degree + 1; 1419 CeedInt numQ_Sur = numP_Sur + qextraSur; 1420 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1421 CeedBasis basisxSur, basisxcSur, basisqSur; 1422 CeedInt NqptsSur; 1423 CeedQFunction qf_setupSur, qf_Sur; 1424 1425 // CEED bases for the boundaries 1426 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, CEED_GAUSS, 1427 &basisqSur); 1428 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1429 &basisxSur); 1430 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1431 CEED_GAUSS_LOBATTO, &basisxcSur); 1432 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1433 1434 // Create the Q-function that builds the quadrature data for the Surface operator 1435 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1436 &qf_setupSur); 1437 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1438 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1439 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1440 1441 // Creat Q-Function for Boundaries 1442 qf_Sur = NULL; 1443 if (problem->applySur) { 1444 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1445 problem->applySur_loc, &qf_Sur); 1446 CeedQFunctionAddInput(qf_Sur, "q", ncompq, CEED_EVAL_INTERP); 1447 CeedQFunctionAddInput(qf_Sur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1448 CeedQFunctionAddInput(qf_Sur, "x", ncompx, CEED_EVAL_INTERP); 1449 CeedQFunctionAddOutput(qf_Sur, "v", ncompq, CEED_EVAL_INTERP); 1450 } 1451 1452 // Create CEED Operator for the whole domain 1453 if (!implicit) 1454 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, user->op_rhs_vol, qf_Sur, qf_setupSur, 1455 height, numP_Sur, numQ_Sur, qdatasizeSur, NqptsSur, basisxSur, 1456 basisqSur, &user->op_rhs); CHKERRQ(ierr); 1457 if (implicit) 1458 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, user->op_ifunction_vol, qf_Sur, qf_setupSur, 1459 height, numP_Sur, numQ_Sur, qdatasizeSur, NqptsSur, basisxSur, 1460 basisqSur, &user->op_ifunction); CHKERRQ(ierr); 1461 // Set up contex for QFunctions 1462 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1463 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1464 struct Advection2dContext_ ctxAdvection2d = { 1465 .CtauS = CtauS, 1466 .strong_form = strong_form, 1467 .stabilization = stab, 1468 }; 1469 struct SurfaceContext_ ctxSurface = { 1470 .cv = cv, 1471 .cp = cp, 1472 .Rd = Rd, 1473 .P_wind = P_wind, 1474 .rho_wind = rho_wind, 1475 .strong_form = strong_form, 1476 .wind[0] = wind[0], 1477 .wind[1] = wind[1], 1478 .wind[2] = wind[2], 1479 .wind_type = wind_type, 1480 .implicit = implicit, 1481 }; 1482 switch (problemChoice) { 1483 case NS_DENSITY_CURRENT: 1484 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1485 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1486 sizeof ctxNS); 1487 break; 1488 case NS_ADVECTION: 1489 case NS_ADVECTION2D: 1490 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1491 sizeof ctxAdvection2d); 1492 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1493 sizeof ctxAdvection2d); 1494 if (qf_Sur) CeedQFunctionSetContext(qf_Sur, &ctxSurface, sizeof ctxSurface); 1495 } 1496 1497 // Set up PETSc context 1498 // Set up units structure 1499 units->meter = meter; 1500 units->kilogram = kilogram; 1501 units->second = second; 1502 units->Kelvin = Kelvin; 1503 units->Pascal = Pascal; 1504 units->JperkgK = JperkgK; 1505 units->mpersquareds = mpersquareds; 1506 units->WpermK = WpermK; 1507 units->kgpercubicm = kgpercubicm; 1508 units->kgpersquaredms = kgpersquaredms; 1509 units->Joulepercubicm = Joulepercubicm; 1510 1511 // Set up user structure 1512 user->comm = comm; 1513 user->outputfreq = outputfreq; 1514 user->contsteps = contsteps; 1515 user->units = units; 1516 user->dm = dm; 1517 user->dmviz = dmviz; 1518 user->interpviz = interpviz; 1519 user->ceed = ceed; 1520 1521 // Calculate qdata and ICs 1522 // Set up state global and local vectors 1523 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1524 1525 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1526 1527 // Apply Setup Ceed Operators 1528 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1529 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1530 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1531 user->M); CHKERRQ(ierr); 1532 1533 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1534 &ctxSetup, 0.0); CHKERRQ(ierr); 1535 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1536 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1537 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1538 // slower execution. 1539 Vec Qbc; 1540 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1541 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1542 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1543 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1544 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1545 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1546 ierr = PetscObjectComposeFunction((PetscObject)dm, 1547 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1548 CHKERRQ(ierr); 1549 } 1550 1551 MPI_Comm_rank(comm, &rank); 1552 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1553 // Gather initial Q values 1554 // In case of continuation of simulation, set up initial values from binary file 1555 if (contsteps) { // continue from existent solution 1556 PetscViewer viewer; 1557 char filepath[PETSC_MAX_PATH_LEN]; 1558 // Read input 1559 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1560 user->outputfolder); 1561 CHKERRQ(ierr); 1562 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1563 CHKERRQ(ierr); 1564 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1565 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1566 } 1567 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1568 1569 // Create and setup TS 1570 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1571 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1572 if (implicit) { 1573 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1574 if (user->op_ifunction) { 1575 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1576 } else { // Implicit integrators can fall back to using an RHSFunction 1577 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1578 } 1579 } else { 1580 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1581 "Problem does not provide RHSFunction"); 1582 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1583 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1584 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1585 } 1586 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1587 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1588 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1589 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1590 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1591 ierr = TSAdaptSetStepLimits(adapt, 1592 1.e-12 * units->second, 1593 1.e2 * units->second); CHKERRQ(ierr); 1594 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1595 if (!contsteps) { // print initial condition 1596 if (testChoice == TEST_NONE) { 1597 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1598 } 1599 } else { // continue from time of last output 1600 PetscReal time; 1601 PetscInt count; 1602 PetscViewer viewer; 1603 char filepath[PETSC_MAX_PATH_LEN]; 1604 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1605 user->outputfolder); CHKERRQ(ierr); 1606 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1607 CHKERRQ(ierr); 1608 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1609 CHKERRQ(ierr); 1610 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1611 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1612 } 1613 if (testChoice == TEST_NONE) { 1614 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1615 } 1616 1617 // Solve 1618 start = MPI_Wtime(); 1619 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1620 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1621 cpu_time_used = MPI_Wtime() - start; 1622 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1623 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1624 comm); CHKERRQ(ierr); 1625 if (testChoice == TEST_NONE) { 1626 ierr = PetscPrintf(PETSC_COMM_WORLD, 1627 "Time taken for solution (sec): %g\n", 1628 (double)cpu_time_used); CHKERRQ(ierr); 1629 } 1630 1631 // Get error 1632 if (problem->non_zero_time && testChoice == TEST_NONE) { 1633 Vec Qexact, Qexactloc; 1634 PetscReal norm; 1635 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1636 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1637 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1638 1639 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1640 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1641 1642 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1643 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1644 CeedVectorDestroy(&q0ceed); 1645 ierr = PetscPrintf(PETSC_COMM_WORLD, 1646 "Max Error: %g\n", 1647 (double)norm); CHKERRQ(ierr); 1648 // Clean up vectors 1649 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1650 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1651 } 1652 1653 // Output Statistics 1654 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1655 if (testChoice == TEST_NONE) { 1656 ierr = PetscPrintf(PETSC_COMM_WORLD, 1657 "Time integrator took %D time steps to reach final time %g\n", 1658 steps, (double)ftime); CHKERRQ(ierr); 1659 } 1660 // Output numerical values from command line 1661 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1662 1663 // compare reference solution values with current run 1664 if (testChoice != TEST_NONE) { 1665 PetscViewer viewer; 1666 // Read reference file 1667 Vec Qref; 1668 PetscReal error, Qrefnorm; 1669 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1670 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1671 CHKERRQ(ierr); 1672 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1673 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1674 1675 // Compute error with respect to reference solution 1676 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1677 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1678 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1679 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1680 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1681 // Check error 1682 if (error > test->testtol) { 1683 ierr = PetscPrintf(PETSC_COMM_WORLD, 1684 "Test failed with error norm %g\n", 1685 (double)error); CHKERRQ(ierr); 1686 } 1687 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1688 } 1689 1690 // Clean up libCEED 1691 CeedVectorDestroy(&qdata); 1692 CeedVectorDestroy(&user->qceed); 1693 CeedVectorDestroy(&user->qdotceed); 1694 CeedVectorDestroy(&user->gceed); 1695 CeedVectorDestroy(&xcorners); 1696 CeedBasisDestroy(&basisq); 1697 CeedBasisDestroy(&basisx); 1698 CeedBasisDestroy(&basisxc); 1699 CeedElemRestrictionDestroy(&restrictq); 1700 CeedElemRestrictionDestroy(&restrictx); 1701 CeedElemRestrictionDestroy(&restrictqdi); 1702 CeedQFunctionDestroy(&qf_setupVol); 1703 CeedQFunctionDestroy(&qf_ics); 1704 CeedQFunctionDestroy(&qf_rhsVol); 1705 CeedQFunctionDestroy(&qf_ifunctionVol); 1706 CeedOperatorDestroy(&op_setupVol); 1707 CeedOperatorDestroy(&op_ics); 1708 CeedOperatorDestroy(&user->op_rhs_vol); 1709 CeedOperatorDestroy(&user->op_ifunction_vol); 1710 CeedDestroy(&ceed); 1711 CeedBasisDestroy(&basisqSur); 1712 CeedBasisDestroy(&basisxSur); 1713 CeedBasisDestroy(&basisxcSur); 1714 CeedQFunctionDestroy(&qf_setupSur); 1715 CeedQFunctionDestroy(&qf_Sur); 1716 CeedOperatorDestroy(&user->op_rhs); 1717 CeedOperatorDestroy(&user->op_ifunction); 1718 1719 // Clean up PETSc 1720 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1721 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1722 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1723 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1724 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1725 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1726 ierr = PetscFree(units); CHKERRQ(ierr); 1727 ierr = PetscFree(user); CHKERRQ(ierr); 1728 return PetscFinalize(); 1729 } 1730 1731