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