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 PetscScalar Joule; 245 }; 246 247 typedef struct SimpleBC_ *SimpleBC; 248 struct SimpleBC_ { 249 PetscInt nwall, nslip[3]; 250 PetscInt walls[6], slips[3][6]; 251 PetscBool userbc; 252 }; 253 254 // Essential BC dofs are encoded in closure indices as -(i+1). 255 static PetscInt Involute(PetscInt i) { 256 return i >= 0 ? i : -(i+1); 257 } 258 259 // Utility function to create local CEED restriction 260 static PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 261 CeedInt height, DMLabel domainLabel, CeedInt value, 262 CeedElemRestriction *Erestrict) { 263 264 PetscSection section; 265 PetscInt p, Nelem, Ndof, *erestrict, eoffset, nfields, dim, depth; 266 DMLabel depthLabel; 267 IS depthIS, iterIS; 268 Vec Uloc; 269 const PetscInt *iterIndices; 270 PetscErrorCode ierr; 271 272 PetscFunctionBeginUser; 273 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 274 dim -= height; 275 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 276 ierr = PetscSectionGetNumFields(section, &nfields); CHKERRQ(ierr); 277 PetscInt ncomp[nfields], fieldoff[nfields+1]; 278 fieldoff[0] = 0; 279 for (PetscInt f=0; f<nfields; f++) { 280 ierr = PetscSectionGetFieldComponents(section, f, &ncomp[f]); CHKERRQ(ierr); 281 fieldoff[f+1] = fieldoff[f] + ncomp[f]; 282 } 283 284 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 285 ierr = DMPlexGetDepthLabel(dm, &depthLabel); CHKERRQ(ierr); 286 ierr = DMLabelGetStratumIS(depthLabel, depth - height, &depthIS); CHKERRQ(ierr); 287 if (domainLabel) { 288 IS domainIS; 289 ierr = DMLabelGetStratumIS(domainLabel, value, &domainIS); CHKERRQ(ierr); 290 ierr = ISIntersect(depthIS, domainIS, &iterIS); CHKERRQ(ierr); 291 ierr = ISDestroy(&domainIS); CHKERRQ(ierr); 292 ierr = ISDestroy(&depthIS); CHKERRQ(ierr); 293 } else { 294 iterIS = depthIS; 295 } 296 ierr = ISGetLocalSize(iterIS, &Nelem); CHKERRQ(ierr); 297 ierr = ISGetIndices(iterIS, &iterIndices); CHKERRQ(ierr); 298 ierr = PetscMalloc1(Nelem*PetscPowInt(P, dim), &erestrict); CHKERRQ(ierr); 299 for (p=0,eoffset=0; p<Nelem; p++) { 300 PetscInt c = iterIndices[p]; 301 PetscInt numindices, *indices, nnodes; 302 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, 303 &numindices, &indices, NULL, NULL); 304 CHKERRQ(ierr); 305 bool flip = false; 306 if (height > 0) { 307 PetscInt numCells, numFaces, start = -1; 308 const PetscInt *orients, *faces, *cells; 309 ierr = DMPlexGetSupport(dm, c, &cells); CHKERRQ(ierr); 310 ierr = DMPlexGetSupportSize(dm, c, &numCells); CHKERRQ(ierr); 311 if (numCells != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected one cell in support of exterior face, but got %D cells", numCells); 312 ierr = DMPlexGetCone(dm, cells[0], &faces); CHKERRQ(ierr); 313 ierr = DMPlexGetConeSize(dm, cells[0], &numFaces); CHKERRQ(ierr); 314 for (PetscInt i=0; i<numFaces; i++) {if (faces[i] == c) start = i;} 315 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %D in cone of its support", c); 316 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients); CHKERRQ(ierr); 317 if (orients[start] < 0) flip = true; 318 } 319 if (numindices % fieldoff[nfields]) SETERRQ1(PETSC_COMM_SELF, 320 PETSC_ERR_ARG_INCOMP, "Number of closure indices not compatible with Cell %D", 321 c); 322 nnodes = numindices / fieldoff[nfields]; 323 for (PetscInt i=0; i<nnodes; i++) { 324 PetscInt ii = i; 325 if (flip) { 326 if (P == nnodes) ii = nnodes - 1 - i; 327 else if (P*P == nnodes) { 328 PetscInt row = i / P, col = i % P; 329 ii = row + col * P; 330 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for flipping point with %D nodes != P (%D) or P^2", nnodes, P); 331 } 332 // Check that indices are blocked by node and thus can be coalesced as a single field with 333 // fieldoff[nfields] = sum(ncomp) components. 334 for (PetscInt f=0; f<nfields; f++) { 335 for (PetscInt j=0; j<ncomp[f]; j++) { 336 if (Involute(indices[fieldoff[f]*nnodes + ii*ncomp[f] + j]) 337 != Involute(indices[ii*ncomp[0]]) + fieldoff[f] + j) 338 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 339 "Cell %D closure indices not interlaced for node %D field %D component %D", 340 c, ii, f, j); 341 } 342 } 343 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 344 PetscInt loc = Involute(indices[ii*ncomp[0]]); 345 erestrict[eoffset++] = loc; 346 } 347 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, 348 &numindices, &indices, NULL, NULL); 349 CHKERRQ(ierr); 350 } 351 if (eoffset != Nelem*PetscPowInt(P, dim)) 352 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 353 "ElemRestriction of size (%D,%D) initialized %D nodes", Nelem, 354 PetscPowInt(P, dim),eoffset); 355 ierr = ISRestoreIndices(iterIS, &iterIndices); CHKERRQ(ierr); 356 ierr = ISDestroy(&iterIS); CHKERRQ(ierr); 357 358 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 359 ierr = VecGetLocalSize(Uloc, &Ndof); CHKERRQ(ierr); 360 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 361 CeedElemRestrictionCreate(ceed, Nelem, PetscPowInt(P, dim), fieldoff[nfields], 362 1, Ndof, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 363 Erestrict); 364 ierr = PetscFree(erestrict); CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367 368 // Utility function to get Ceed Restriction for each domain 369 static PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 370 DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, CeedInt qdatasize, 371 CeedElemRestriction *restrictq, CeedElemRestriction *restrictx, 372 CeedElemRestriction *restrictqdi) { 373 374 DM dmcoord; 375 CeedInt dim, localNelem; 376 CeedInt Qdim; 377 PetscErrorCode ierr; 378 379 PetscFunctionBeginUser; 380 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 381 dim -= height; 382 Qdim = CeedIntPow(Q, dim); 383 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 384 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 385 CHKERRQ(ierr); 386 ierr = CreateRestrictionFromPlex(ceed, dm, P, height, domainLabel, value, restrictq); 387 CHKERRQ(ierr); 388 ierr = CreateRestrictionFromPlex(ceed, dmcoord, 2, height, domainLabel, value, restrictx); 389 CHKERRQ(ierr); 390 CeedElemRestrictionGetNumElements(*restrictq, &localNelem); 391 CeedElemRestrictionCreateStrided(ceed, localNelem, Qdim, 392 qdatasize, qdatasize*localNelem*Qdim, 393 CEED_STRIDES_BACKEND, restrictqdi); 394 PetscFunctionReturn(0); 395 } 396 397 // Utility function to create CEED Composite Operator for the entire domain 398 static PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, 399 WindType wind_type, CeedOperator op_applyVol, CeedQFunction qf_applySur, 400 CeedQFunction qf_setupSur, CeedInt height, CeedInt numP_Sur, CeedInt numQ_Sur, 401 CeedInt qdatasizeSur, CeedInt NqptsSur, CeedBasis basisxSur, 402 CeedBasis basisqSur, CeedOperator *op_apply) { 403 404 CeedInt dim, nFace; 405 PetscInt lsize, localNelemSur[6]; 406 Vec Xloc; 407 CeedVector xcorners, qdataSur[6]; 408 CeedOperator op_setupSur[6], op_applySur[6]; 409 CeedElemRestriction restrictxSur[6], restrictqSur[6], restrictqdiSur[6]; 410 DMLabel domainLabel; 411 PetscScalar *x; 412 PetscErrorCode ierr; 413 414 PetscFunctionBeginUser; 415 // Composite Operaters 416 CeedCompositeOperatorCreate(ceed, op_apply); 417 CeedCompositeOperatorAddSub(*op_apply, op_applyVol); // Apply a Sub-Operator for the volume 418 419 if (wind_type == ADVECTION_WIND_TRANSLATION) { 420 bc->nwall = 0; 421 bc->nslip[0] = bc->nslip[1] = bc->nslip[2] = 0; 422 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 423 ierr = VecGetLocalSize(Xloc, &lsize); CHKERRQ(ierr); 424 ierr = CeedVectorCreate(ceed, lsize, &xcorners); CHKERRQ(ierr); 425 ierr = VecGetArray(Xloc, &x); CHKERRQ(ierr); 426 CeedVectorSetArray(xcorners, CEED_MEM_HOST, CEED_USE_POINTER, x); 427 ierr = DMGetLabel(dm, "Face Sets", &domainLabel); CHKERRQ(ierr); 428 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 429 if (dim == 2) nFace = 4; 430 if (dim == 3) nFace = 6; 431 432 // Create CEED Operator for each boundary face 433 for (CeedInt i=0; i<nFace; i++) { 434 ierr = GetRestrictionForDomain(ceed, dm, height, domainLabel, i+1, numP_Sur, 435 numQ_Sur, qdatasizeSur, &restrictqSur[i], &restrictxSur[i], 436 &restrictqdiSur[i]); CHKERRQ(ierr); 437 // Create the CEED vectors that will be needed in Boundary setup 438 CeedElemRestrictionGetNumElements(restrictqSur[i], &localNelemSur[i]); 439 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, &qdataSur[i]); 440 // Create the operator that builds the quadrature data for the Boundary operator 441 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur[i]); 442 CeedOperatorSetField(op_setupSur[i], "dx", restrictxSur[i], basisxSur, CEED_VECTOR_ACTIVE); 443 CeedOperatorSetField(op_setupSur[i], "weight", CEED_ELEMRESTRICTION_NONE, 444 basisxSur, CEED_VECTOR_NONE); 445 CeedOperatorSetField(op_setupSur[i], "qdataSur", restrictqdiSur[i], 446 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 447 // Create Boundary operator 448 CeedOperatorCreate(ceed, qf_applySur, NULL, NULL, &op_applySur[i]); 449 CeedOperatorSetField(op_applySur[i], "q", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 450 CeedOperatorSetField(op_applySur[i], "qdataSur", restrictqdiSur[i], 451 CEED_BASIS_COLLOCATED, qdataSur[i]); 452 CeedOperatorSetField(op_applySur[i], "x", restrictxSur[i], basisxSur, xcorners); 453 CeedOperatorSetField(op_applySur[i], "v", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 454 // Apply CEED operator for Boundary setup 455 CeedOperatorApply(op_setupSur[i], xcorners, qdataSur[i], CEED_REQUEST_IMMEDIATE); 456 // Apply Sub-Operator for the Boundary 457 CeedCompositeOperatorAddSub(*op_apply, op_applySur[i]); 458 } 459 CeedVectorDestroy(&xcorners); 460 } 461 PetscFunctionReturn(0); 462 } 463 464 static int CreateVectorFromPetscVec(Ceed ceed, Vec p, CeedVector *v) { 465 PetscErrorCode ierr; 466 PetscInt m; 467 468 PetscFunctionBeginUser; 469 ierr = VecGetLocalSize(p, &m); CHKERRQ(ierr); 470 ierr = CeedVectorCreate(ceed, m, v); CHKERRQ(ierr); 471 PetscFunctionReturn(0); 472 } 473 474 static int VectorPlacePetscVec(CeedVector c, Vec p) { 475 PetscErrorCode ierr; 476 PetscInt mceed, mpetsc; 477 PetscScalar *a; 478 479 PetscFunctionBeginUser; 480 ierr = CeedVectorGetLength(c, &mceed); CHKERRQ(ierr); 481 ierr = VecGetLocalSize(p, &mpetsc); CHKERRQ(ierr); 482 if (mceed != mpetsc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 483 "Cannot place PETSc Vec of length %D in CeedVector of length %D", 484 mpetsc, mceed); 485 ierr = VecGetArray(p, &a); CHKERRQ(ierr); 486 CeedVectorSetArray(c, CEED_MEM_HOST, CEED_USE_POINTER, a); 487 PetscFunctionReturn(0); 488 } 489 490 static PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 491 PetscBool insertEssential, Vec Qloc, PetscReal time, Vec faceGeomFVM, 492 Vec cellGeomFVM, Vec gradFVM) { 493 PetscErrorCode ierr; 494 Vec Qbc; 495 496 PetscFunctionBegin; 497 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 498 ierr = VecAXPY(Qloc, 1., Qbc); CHKERRQ(ierr); 499 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 500 PetscFunctionReturn(0); 501 } 502 503 // This is the RHS of the ODE, given as u_t = G(t,u) 504 // This function takes in a state vector Q and writes into G 505 static PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *userData) { 506 PetscErrorCode ierr; 507 User user = *(User *)userData; 508 PetscScalar *q, *g; 509 Vec Qloc, Gloc; 510 511 // Global-to-local 512 PetscFunctionBeginUser; 513 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 514 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 515 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 516 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 517 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 518 NULL, NULL, NULL); CHKERRQ(ierr); 519 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 520 521 // Ceed Vectors 522 ierr = VecGetArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 523 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 524 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); 525 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 526 527 // Apply CEED operator 528 CeedOperatorApply(user->op_rhs, user->qceed, user->gceed, 529 CEED_REQUEST_IMMEDIATE); 530 531 // Restore vectors 532 ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 533 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 534 535 ierr = VecZeroEntries(G); CHKERRQ(ierr); 536 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 537 538 // Inverse of the lumped mass matrix 539 ierr = VecPointwiseMult(G, G, user->M); // M is Minv 540 CHKERRQ(ierr); 541 542 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 543 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 544 PetscFunctionReturn(0); 545 } 546 547 static PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Qdot, Vec G, 548 void *userData) { 549 PetscErrorCode ierr; 550 User user = *(User *)userData; 551 const PetscScalar *q, *qdot; 552 PetscScalar *g; 553 Vec Qloc, Qdotloc, Gloc; 554 555 // Global-to-local 556 PetscFunctionBeginUser; 557 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 558 ierr = DMGetLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 559 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 560 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 561 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 562 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 563 NULL, NULL, NULL); CHKERRQ(ierr); 564 ierr = VecZeroEntries(Qdotloc); CHKERRQ(ierr); 565 ierr = DMGlobalToLocal(user->dm, Qdot, INSERT_VALUES, Qdotloc); CHKERRQ(ierr); 566 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 567 568 // Ceed Vectors 569 ierr = VecGetArrayRead(Qloc, &q); CHKERRQ(ierr); 570 ierr = VecGetArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 571 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 572 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, 573 (PetscScalar *)q); 574 CeedVectorSetArray(user->qdotceed, CEED_MEM_HOST, CEED_USE_POINTER, 575 (PetscScalar *)qdot); 576 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 577 578 // Apply CEED operator 579 CeedOperatorApply(user->op_ifunction, user->qceed, user->gceed, 580 CEED_REQUEST_IMMEDIATE); 581 582 // Restore vectors 583 ierr = VecRestoreArrayRead(Qloc, &q); CHKERRQ(ierr); 584 ierr = VecRestoreArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 585 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 586 587 ierr = VecZeroEntries(G); CHKERRQ(ierr); 588 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 589 590 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 591 ierr = DMRestoreLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 592 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 593 PetscFunctionReturn(0); 594 } 595 596 // User provided TS Monitor 597 static PetscErrorCode TSMonitor_NS(TS ts, PetscInt stepno, PetscReal time, 598 Vec Q, void *ctx) { 599 User user = ctx; 600 Vec Qloc; 601 char filepath[PETSC_MAX_PATH_LEN]; 602 PetscViewer viewer; 603 PetscErrorCode ierr; 604 605 // Set up output 606 PetscFunctionBeginUser; 607 // Print every 'outputfreq' steps 608 if (stepno % user->outputfreq != 0) 609 PetscFunctionReturn(0); 610 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 611 ierr = PetscObjectSetName((PetscObject)Qloc, "StateVec"); CHKERRQ(ierr); 612 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 613 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 614 615 // Output 616 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-%03D.vtu", 617 user->outputfolder, stepno + user->contsteps); 618 CHKERRQ(ierr); 619 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), filepath, 620 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 621 ierr = VecView(Qloc, viewer); CHKERRQ(ierr); 622 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 623 if (user->dmviz) { 624 Vec Qrefined, Qrefined_loc; 625 char filepath_refined[PETSC_MAX_PATH_LEN]; 626 PetscViewer viewer_refined; 627 628 ierr = DMGetGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 629 ierr = DMGetLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 630 ierr = PetscObjectSetName((PetscObject)Qrefined_loc, "Refined"); 631 CHKERRQ(ierr); 632 ierr = MatInterpolate(user->interpviz, Q, Qrefined); CHKERRQ(ierr); 633 ierr = VecZeroEntries(Qrefined_loc); CHKERRQ(ierr); 634 ierr = DMGlobalToLocal(user->dmviz, Qrefined, INSERT_VALUES, Qrefined_loc); 635 CHKERRQ(ierr); 636 ierr = PetscSNPrintf(filepath_refined, sizeof filepath_refined, 637 "%s/nsrefined-%03D.vtu", 638 user->outputfolder, stepno + user->contsteps); 639 CHKERRQ(ierr); 640 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Qrefined), 641 filepath_refined, 642 FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 643 ierr = VecView(Qrefined_loc, viewer_refined); CHKERRQ(ierr); 644 ierr = DMRestoreLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 645 ierr = DMRestoreGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 646 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 647 } 648 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 649 650 // Save data in a binary file for continuation of simulations 651 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 652 user->outputfolder); CHKERRQ(ierr); 653 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 654 CHKERRQ(ierr); 655 ierr = VecView(Q, viewer); CHKERRQ(ierr); 656 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 657 658 // Save time stamp 659 // Dimensionalize time back 660 time /= user->units->second; 661 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 662 user->outputfolder); CHKERRQ(ierr); 663 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 664 CHKERRQ(ierr); 665 #if PETSC_VERSION_GE(3,13,0) 666 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 667 #else 668 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 669 #endif 670 CHKERRQ(ierr); 671 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 672 673 PetscFunctionReturn(0); 674 } 675 676 static PetscErrorCode ICs_FixMultiplicity(CeedOperator op_ics, 677 CeedVector xcorners, CeedVector q0ceed, DM dm, Vec Qloc, Vec Q, 678 CeedElemRestriction restrictq, SetupContext ctxSetup, CeedScalar time) { 679 PetscErrorCode ierr; 680 CeedVector multlvec; 681 Vec Multiplicity, MultiplicityLoc; 682 683 ctxSetup->time = time; 684 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 685 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 686 CeedOperatorApply(op_ics, xcorners, q0ceed, CEED_REQUEST_IMMEDIATE); 687 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 688 ierr = DMLocalToGlobal(dm, Qloc, ADD_VALUES, Q); CHKERRQ(ierr); 689 690 // Fix multiplicity for output of ICs 691 ierr = DMGetLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 692 CeedElemRestrictionCreateVector(restrictq, &multlvec, NULL); 693 ierr = VectorPlacePetscVec(multlvec, MultiplicityLoc); CHKERRQ(ierr); 694 CeedElemRestrictionGetMultiplicity(restrictq, multlvec); 695 CeedVectorDestroy(&multlvec); 696 ierr = DMGetGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 697 ierr = VecZeroEntries(Multiplicity); CHKERRQ(ierr); 698 ierr = DMLocalToGlobal(dm, MultiplicityLoc, ADD_VALUES, Multiplicity); 699 CHKERRQ(ierr); 700 ierr = VecPointwiseDivide(Q, Q, Multiplicity); CHKERRQ(ierr); 701 ierr = VecPointwiseDivide(Qloc, Qloc, MultiplicityLoc); CHKERRQ(ierr); 702 ierr = DMRestoreLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 703 ierr = DMRestoreGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 704 705 PetscFunctionReturn(0); 706 } 707 708 static PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, 709 CeedElemRestriction restrictq, CeedBasis basisq, 710 CeedElemRestriction restrictqdi, CeedVector qdata, Vec M) { 711 PetscErrorCode ierr; 712 CeedQFunction qf_mass; 713 CeedOperator op_mass; 714 CeedVector mceed; 715 Vec Mloc; 716 CeedInt ncompq, qdatasize; 717 718 PetscFunctionBeginUser; 719 CeedElemRestrictionGetNumComponents(restrictq, &ncompq); 720 CeedElemRestrictionGetNumComponents(restrictqdi, &qdatasize); 721 // Create the Q-function that defines the action of the mass operator 722 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 723 CeedQFunctionAddInput(qf_mass, "q", ncompq, CEED_EVAL_INTERP); 724 CeedQFunctionAddInput(qf_mass, "qdata", qdatasize, CEED_EVAL_NONE); 725 CeedQFunctionAddOutput(qf_mass, "v", ncompq, CEED_EVAL_INTERP); 726 727 // Create the mass operator 728 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 729 CeedOperatorSetField(op_mass, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 730 CeedOperatorSetField(op_mass, "qdata", restrictqdi, 731 CEED_BASIS_COLLOCATED, qdata); 732 CeedOperatorSetField(op_mass, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 733 734 ierr = DMGetLocalVector(dm, &Mloc); CHKERRQ(ierr); 735 ierr = VecZeroEntries(Mloc); CHKERRQ(ierr); 736 CeedElemRestrictionCreateVector(restrictq, &mceed, NULL); 737 ierr = VectorPlacePetscVec(mceed, Mloc); CHKERRQ(ierr); 738 739 { 740 // Compute a lumped mass matrix 741 CeedVector onesvec; 742 CeedElemRestrictionCreateVector(restrictq, &onesvec, NULL); 743 CeedVectorSetValue(onesvec, 1.0); 744 CeedOperatorApply(op_mass, onesvec, mceed, CEED_REQUEST_IMMEDIATE); 745 CeedVectorDestroy(&onesvec); 746 CeedOperatorDestroy(&op_mass); 747 CeedVectorDestroy(&mceed); 748 } 749 CeedQFunctionDestroy(&qf_mass); 750 751 ierr = VecZeroEntries(M); CHKERRQ(ierr); 752 ierr = DMLocalToGlobal(dm, Mloc, ADD_VALUES, M); CHKERRQ(ierr); 753 ierr = DMRestoreLocalVector(dm, &Mloc); CHKERRQ(ierr); 754 755 // Invert diagonally lumped mass vector for RHS function 756 ierr = VecReciprocal(M); CHKERRQ(ierr); 757 PetscFunctionReturn(0); 758 } 759 760 static PetscErrorCode SetUpDM(DM dm, problemData *problem, PetscInt degree, 761 SimpleBC bc, void *ctxSetup) { 762 PetscErrorCode ierr; 763 764 PetscFunctionBeginUser; 765 { 766 // Configure the finite element space and boundary conditions 767 PetscFE fe; 768 PetscInt ncompq = 5; 769 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, ncompq, 770 PETSC_FALSE, degree, PETSC_DECIDE, 771 &fe); CHKERRQ(ierr); 772 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 773 ierr = DMAddField(dm, NULL,(PetscObject)fe); CHKERRQ(ierr); 774 ierr = DMCreateDS(dm); CHKERRQ(ierr); 775 { 776 PetscInt comps[1] = {1}; 777 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 778 1, comps, (void(*)(void))NULL, bc->nslip[0], 779 bc->slips[0], ctxSetup); CHKERRQ(ierr); 780 comps[0] = 2; 781 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 782 1, comps, (void(*)(void))NULL, bc->nslip[1], 783 bc->slips[1], ctxSetup); CHKERRQ(ierr); 784 comps[0] = 3; 785 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 786 1, comps, (void(*)(void))NULL, bc->nslip[2], 787 bc->slips[2], ctxSetup); CHKERRQ(ierr); 788 } 789 if (bc->userbc == PETSC_TRUE) { 790 for (PetscInt c = 0; c < 3; c++) { 791 for (PetscInt s = 0; s < bc->nslip[c]; s++) { 792 for (PetscInt w = 0; w < bc->nwall; w++) { 793 if (bc->slips[c][s] == bc->walls[w]) 794 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, 795 "Boundary condition already set on face %D!\n", bc->walls[w]); 796 797 } 798 } 799 } 800 } 801 // Wall boundary conditions are zero energy density and zero flux for 802 // velocity in advection/advection2d, and zero velocity and zero flux 803 // for mass density and energy density in density_current 804 { 805 if (problem->bc == Exact_Advection || problem->bc == Exact_Advection2d) { 806 PetscInt comps[1] = {4}; 807 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 808 1, comps, (void(*)(void))problem->bc, 809 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 810 } else if (problem->bc == Exact_DC) { 811 PetscInt comps[3] = {1, 2, 3}; 812 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 813 3, comps, (void(*)(void))problem->bc, 814 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 815 } else 816 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, 817 "Undefined boundary conditions for this problem"); 818 } 819 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 820 CHKERRQ(ierr); 821 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 822 } 823 { 824 // Empty name for conserved field (because there is only one field) 825 PetscSection section; 826 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 827 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 828 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 829 CHKERRQ(ierr); 830 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 831 CHKERRQ(ierr); 832 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 833 CHKERRQ(ierr); 834 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 835 CHKERRQ(ierr); 836 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 837 CHKERRQ(ierr); 838 } 839 PetscFunctionReturn(0); 840 } 841 842 int main(int argc, char **argv) { 843 PetscInt ierr; 844 MPI_Comm comm; 845 DM dm, dmcoord, dmviz; 846 Mat interpviz; 847 TS ts; 848 TSAdapt adapt; 849 User user; 850 Units units; 851 char ceedresource[4096] = "/cpu/self"; 852 PetscInt localNelemVol, lnodes, gnodes, steps; 853 const PetscInt ncompq = 5; 854 PetscMPIInt rank; 855 PetscScalar ftime; 856 Vec Q, Qloc, Xloc; 857 Ceed ceed; 858 CeedInt numP, numQ; 859 CeedVector xcorners, qdata, q0ceed; 860 CeedBasis basisx, basisxc, basisq; 861 CeedElemRestriction restrictx, restrictq, restrictqdi; 862 CeedQFunction qf_setupVol, qf_ics, qf_rhsVol, qf_ifunctionVol; 863 CeedOperator op_setupVol, op_ics; 864 CeedScalar Rd; 865 CeedMemType memtyperequested; 866 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 867 kgpersquaredms, Joulepercubicm, Joule; 868 problemType problemChoice; 869 problemData *problem = NULL; 870 WindType wind_type; 871 StabilizationType stab; 872 testType testChoice; 873 testData *test = NULL; 874 PetscBool implicit; 875 PetscInt viz_refine = 0; 876 struct SimpleBC_ bc = { 877 .nslip = {2, 2, 2}, 878 .slips = {{5, 6}, {3, 4}, {1, 2}} 879 }; 880 double start, cpu_time_used; 881 // Check PETSc CUDA support 882 PetscBool petschavecuda, setmemtyperequest = PETSC_FALSE; 883 // *INDENT-OFF* 884 #ifdef PETSC_HAVE_CUDA 885 petschavecuda = PETSC_TRUE; 886 #else 887 petschavecuda = PETSC_FALSE; 888 #endif 889 // *INDENT-ON* 890 891 // Create the libCEED contexts 892 PetscScalar meter = 1e-2; // 1 meter in scaled length units 893 PetscScalar second = 1e-2; // 1 second in scaled time units 894 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 895 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 896 CeedScalar theta0 = 300.; // K 897 CeedScalar thetaC = -15.; // K 898 CeedScalar P0 = 1.e5; // Pa 899 CeedScalar E_wind = 1.e6; // J 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("-E_wind", "Total energy of inflow wind", 1025 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 1026 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1027 NULL, N, &N, NULL); CHKERRQ(ierr); 1028 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1029 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1030 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1031 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1032 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1033 NULL, g, &g, NULL); CHKERRQ(ierr); 1034 ierr = PetscOptionsScalar("-lambda", 1035 "Stokes hypothesis second viscosity coefficient", 1036 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1037 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1038 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1039 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1040 NULL, k, &k, NULL); CHKERRQ(ierr); 1041 ierr = PetscOptionsScalar("-CtauS", 1042 "Scale coefficient for tau (nondimensional)", 1043 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1044 if (stab == STAB_NONE && CtauS != 0) { 1045 ierr = PetscPrintf(comm, 1046 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1047 CHKERRQ(ierr); 1048 } 1049 ierr = PetscOptionsScalar("-strong_form", 1050 "Strong (1) or weak/integrated by parts (0) advection residual", 1051 NULL, strong_form, &strong_form, NULL); 1052 CHKERRQ(ierr); 1053 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1054 ierr = PetscPrintf(comm, 1055 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1056 CHKERRQ(ierr); 1057 } 1058 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1059 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1060 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1061 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1062 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1063 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1064 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1065 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1066 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1067 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1068 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1069 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1070 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1071 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1072 n = problem->dim; 1073 center[0] = 0.5 * lx; 1074 center[1] = 0.5 * ly; 1075 center[2] = 0.5 * lz; 1076 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1077 NULL, center, &n, NULL); CHKERRQ(ierr); 1078 n = problem->dim; 1079 ierr = PetscOptionsRealArray("-dc_axis", 1080 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1081 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1082 { 1083 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1084 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1085 if (norm > 0) { 1086 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1087 } 1088 } 1089 ierr = PetscOptionsInt("-output_freq", 1090 "Frequency of output, in number of steps", 1091 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1092 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1093 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1094 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1095 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1096 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1097 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1098 PetscBool userQextraSur; 1099 ierr = PetscOptionsInt("-qextra_boundary", "Number of extra quadrature points on in/outflow faces", 1100 NULL, qextraSur, &qextraSur, &userQextraSur); CHKERRQ(ierr); 1101 if ((wind_type == ADVECTION_WIND_ROTATION || problemChoice == NS_DENSITY_CURRENT) && userQextraSur) { 1102 ierr = PetscPrintf(comm, "Warning! Use -qextra_boundary only with -problem_advection_wind translation\n"); 1103 CHKERRQ(ierr); 1104 } 1105 ierr = PetscStrncpy(user->outputfolder, ".", 2); CHKERRQ(ierr); 1106 ierr = PetscOptionsString("-of", "Output folder", 1107 NULL, user->outputfolder, user->outputfolder, 1108 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 1109 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1110 ierr = PetscOptionsEnum("-memtype", 1111 "CEED MemType requested", NULL, 1112 memTypes, (PetscEnum)memtyperequested, 1113 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1114 CHKERRQ(ierr); 1115 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1116 1117 // Define derived units 1118 Pascal = kilogram / (meter * PetscSqr(second)); 1119 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1120 mpersquareds = meter / PetscSqr(second); 1121 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1122 kgpercubicm = kilogram / pow(meter,3); 1123 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1124 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1125 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 1126 1127 // Scale variables to desired units 1128 theta0 *= Kelvin; 1129 thetaC *= Kelvin; 1130 P0 *= Pascal; 1131 E_wind *= Joule; 1132 N *= (1./second); 1133 cv *= JperkgK; 1134 cp *= JperkgK; 1135 Rd = cp - cv; 1136 g *= mpersquareds; 1137 mu *= Pascal * second; 1138 k *= WpermK; 1139 lx = fabs(lx) * meter; 1140 ly = fabs(ly) * meter; 1141 lz = fabs(lz) * meter; 1142 rc = fabs(rc) * meter; 1143 resx = fabs(resx) * meter; 1144 resy = fabs(resy) * meter; 1145 resz = fabs(resz) * meter; 1146 for (int i=0; i<3; i++) center[i] *= meter; 1147 1148 const CeedInt dim = problem->dim, ncompx = problem->dim, 1149 qdatasizeVol = problem->qdatasizeVol; 1150 // Set up the libCEED context 1151 struct SetupContext_ ctxSetup = { 1152 .theta0 = theta0, 1153 .thetaC = thetaC, 1154 .P0 = P0, 1155 .N = N, 1156 .cv = cv, 1157 .cp = cp, 1158 .Rd = Rd, 1159 .g = g, 1160 .rc = rc, 1161 .lx = lx, 1162 .ly = ly, 1163 .lz = lz, 1164 .center[0] = center[0], 1165 .center[1] = center[1], 1166 .center[2] = center[2], 1167 .dc_axis[0] = dc_axis[0], 1168 .dc_axis[1] = dc_axis[1], 1169 .dc_axis[2] = dc_axis[2], 1170 .wind[0] = wind[0], 1171 .wind[1] = wind[1], 1172 .wind[2] = wind[2], 1173 .time = 0, 1174 .wind_type = wind_type, 1175 }; 1176 1177 // Create the mesh 1178 { 1179 const PetscReal scale[3] = {lx, ly, lz}; 1180 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1181 NULL, PETSC_TRUE, &dm); 1182 CHKERRQ(ierr); 1183 } 1184 1185 // Distribute the mesh over processes 1186 { 1187 DM dmDist = NULL; 1188 PetscPartitioner part; 1189 1190 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1191 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1192 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1193 if (dmDist) { 1194 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1195 dm = dmDist; 1196 } 1197 } 1198 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1199 1200 // Setup DM 1201 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1202 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1203 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetup); CHKERRQ(ierr); 1204 1205 // Refine DM for high-order viz 1206 dmviz = NULL; 1207 interpviz = NULL; 1208 if (viz_refine) { 1209 DM dmhierarchy[viz_refine+1]; 1210 1211 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1212 dmhierarchy[0] = dm; 1213 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1214 Mat interp_next; 1215 1216 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1217 CHKERRQ(ierr); 1218 ierr = DMClearDS(dmhierarchy[i+1]);CHKERRQ(ierr); 1219 ierr = DMClearFields(dmhierarchy[i+1]);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_applySur; 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_applySur = NULL; 1447 if (problem->applySur) { 1448 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1449 problem->applySur_loc, &qf_applySur); 1450 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1451 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1452 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1453 CeedQFunctionAddOutput(qf_applySur, "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_applySur, 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_applySur, 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 .E_wind = E_wind, 1475 .strong_form = strong_form, 1476 .implicit = implicit, 1477 }; 1478 switch (problemChoice) { 1479 case NS_DENSITY_CURRENT: 1480 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1481 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1482 sizeof ctxNS); 1483 break; 1484 case NS_ADVECTION: 1485 case NS_ADVECTION2D: 1486 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1487 sizeof ctxAdvection2d); 1488 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1489 sizeof ctxAdvection2d); 1490 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &ctxSurface, sizeof ctxSurface); 1491 } 1492 1493 // Set up PETSc context 1494 // Set up units structure 1495 units->meter = meter; 1496 units->kilogram = kilogram; 1497 units->second = second; 1498 units->Kelvin = Kelvin; 1499 units->Pascal = Pascal; 1500 units->JperkgK = JperkgK; 1501 units->mpersquareds = mpersquareds; 1502 units->WpermK = WpermK; 1503 units->kgpercubicm = kgpercubicm; 1504 units->kgpersquaredms = kgpersquaredms; 1505 units->Joulepercubicm = Joulepercubicm; 1506 units->Joule = Joule; 1507 1508 // Set up user structure 1509 user->comm = comm; 1510 user->outputfreq = outputfreq; 1511 user->contsteps = contsteps; 1512 user->units = units; 1513 user->dm = dm; 1514 user->dmviz = dmviz; 1515 user->interpviz = interpviz; 1516 user->ceed = ceed; 1517 1518 // Calculate qdata and ICs 1519 // Set up state global and local vectors 1520 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1521 1522 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1523 1524 // Apply Setup Ceed Operators 1525 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1526 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1527 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1528 user->M); CHKERRQ(ierr); 1529 1530 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1531 &ctxSetup, 0.0); CHKERRQ(ierr); 1532 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1533 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1534 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1535 // slower execution. 1536 Vec Qbc; 1537 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1538 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1539 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1540 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1541 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1542 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1543 ierr = PetscObjectComposeFunction((PetscObject)dm, 1544 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1545 CHKERRQ(ierr); 1546 } 1547 1548 MPI_Comm_rank(comm, &rank); 1549 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1550 // Gather initial Q values 1551 // In case of continuation of simulation, set up initial values from binary file 1552 if (contsteps) { // continue from existent solution 1553 PetscViewer viewer; 1554 char filepath[PETSC_MAX_PATH_LEN]; 1555 // Read input 1556 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1557 user->outputfolder); 1558 CHKERRQ(ierr); 1559 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1560 CHKERRQ(ierr); 1561 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1562 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1563 } 1564 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1565 1566 // Create and setup TS 1567 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1568 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1569 if (implicit) { 1570 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1571 if (user->op_ifunction) { 1572 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1573 } else { // Implicit integrators can fall back to using an RHSFunction 1574 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1575 } 1576 } else { 1577 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1578 "Problem does not provide RHSFunction"); 1579 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1580 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1581 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1582 } 1583 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1584 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1585 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1586 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1587 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1588 ierr = TSAdaptSetStepLimits(adapt, 1589 1.e-12 * units->second, 1590 1.e2 * units->second); CHKERRQ(ierr); 1591 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1592 if (!contsteps) { // print initial condition 1593 if (testChoice == TEST_NONE) { 1594 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1595 } 1596 } else { // continue from time of last output 1597 PetscReal time; 1598 PetscInt count; 1599 PetscViewer viewer; 1600 char filepath[PETSC_MAX_PATH_LEN]; 1601 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1602 user->outputfolder); CHKERRQ(ierr); 1603 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1604 CHKERRQ(ierr); 1605 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1606 CHKERRQ(ierr); 1607 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1608 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1609 } 1610 if (testChoice == TEST_NONE) { 1611 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1612 } 1613 1614 // Solve 1615 start = MPI_Wtime(); 1616 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1617 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1618 cpu_time_used = MPI_Wtime() - start; 1619 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1620 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1621 comm); CHKERRQ(ierr); 1622 if (testChoice == TEST_NONE) { 1623 ierr = PetscPrintf(PETSC_COMM_WORLD, 1624 "Time taken for solution (sec): %g\n", 1625 (double)cpu_time_used); CHKERRQ(ierr); 1626 } 1627 1628 // Get error 1629 if (problem->non_zero_time && testChoice == TEST_NONE) { 1630 Vec Qexact, Qexactloc; 1631 PetscReal norm; 1632 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1633 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1634 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1635 1636 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1637 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1638 1639 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1640 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1641 CeedVectorDestroy(&q0ceed); 1642 ierr = PetscPrintf(PETSC_COMM_WORLD, 1643 "Max Error: %g\n", 1644 (double)norm); CHKERRQ(ierr); 1645 // Clean up vectors 1646 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1647 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1648 } 1649 1650 // Output Statistics 1651 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1652 if (testChoice == TEST_NONE) { 1653 ierr = PetscPrintf(PETSC_COMM_WORLD, 1654 "Time integrator took %D time steps to reach final time %g\n", 1655 steps, (double)ftime); CHKERRQ(ierr); 1656 } 1657 // Output numerical values from command line 1658 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1659 1660 // compare reference solution values with current run 1661 if (testChoice != TEST_NONE) { 1662 PetscViewer viewer; 1663 // Read reference file 1664 Vec Qref; 1665 PetscReal error, Qrefnorm; 1666 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1667 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1668 CHKERRQ(ierr); 1669 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1670 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1671 1672 // Compute error with respect to reference solution 1673 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1674 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1675 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1676 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1677 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1678 // Check error 1679 if (error > test->testtol) { 1680 ierr = PetscPrintf(PETSC_COMM_WORLD, 1681 "Test failed with error norm %g\n", 1682 (double)error); CHKERRQ(ierr); 1683 } 1684 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1685 } 1686 1687 // Clean up libCEED 1688 CeedVectorDestroy(&qdata); 1689 CeedVectorDestroy(&user->qceed); 1690 CeedVectorDestroy(&user->qdotceed); 1691 CeedVectorDestroy(&user->gceed); 1692 CeedVectorDestroy(&xcorners); 1693 CeedBasisDestroy(&basisq); 1694 CeedBasisDestroy(&basisx); 1695 CeedBasisDestroy(&basisxc); 1696 CeedElemRestrictionDestroy(&restrictq); 1697 CeedElemRestrictionDestroy(&restrictx); 1698 CeedElemRestrictionDestroy(&restrictqdi); 1699 CeedQFunctionDestroy(&qf_setupVol); 1700 CeedQFunctionDestroy(&qf_ics); 1701 CeedQFunctionDestroy(&qf_rhsVol); 1702 CeedQFunctionDestroy(&qf_ifunctionVol); 1703 CeedOperatorDestroy(&op_setupVol); 1704 CeedOperatorDestroy(&op_ics); 1705 CeedOperatorDestroy(&user->op_rhs_vol); 1706 CeedOperatorDestroy(&user->op_ifunction_vol); 1707 CeedDestroy(&ceed); 1708 CeedBasisDestroy(&basisqSur); 1709 CeedBasisDestroy(&basisxSur); 1710 CeedBasisDestroy(&basisxcSur); 1711 CeedQFunctionDestroy(&qf_setupSur); 1712 CeedQFunctionDestroy(&qf_applySur); 1713 CeedOperatorDestroy(&user->op_rhs); 1714 CeedOperatorDestroy(&user->op_ifunction); 1715 1716 // Clean up PETSc 1717 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1718 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1719 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1720 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1721 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1722 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1723 ierr = PetscFree(units); CHKERRQ(ierr); 1724 ierr = PetscFree(user); CHKERRQ(ierr); 1725 return PetscFinalize(); 1726 } 1727 1728