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 = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1219 d = (d + 1) / 2; 1220 if (i + 1 == viz_refine) d = 1; 1221 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 1222 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1223 &interp_next, NULL); CHKERRQ(ierr); 1224 if (!i) interpviz = interp_next; 1225 else { 1226 Mat C; 1227 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1228 PETSC_DECIDE, &C); CHKERRQ(ierr); 1229 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1230 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1231 interpviz = C; 1232 } 1233 } 1234 for (PetscInt i=1; i<viz_refine; i++) { 1235 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1236 } 1237 dmviz = dmhierarchy[viz_refine]; 1238 } 1239 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1240 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1241 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1242 lnodes /= ncompq; 1243 1244 // Initialize CEED 1245 CeedInit(ceedresource, &ceed); 1246 // Set memtype 1247 CeedMemType memtypebackend; 1248 CeedGetPreferredMemType(ceed, &memtypebackend); 1249 // Check memtype compatibility 1250 if (!setmemtyperequest) 1251 memtyperequested = memtypebackend; 1252 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1253 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1254 "PETSc was not built with CUDA. " 1255 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1256 1257 // Set number of 1D nodes and quadrature points 1258 numP = degree + 1; 1259 numQ = numP + qextra; 1260 1261 // Print summary 1262 if (testChoice == TEST_NONE) { 1263 CeedInt gdofs, odofs; 1264 int comm_size; 1265 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1266 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1267 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1268 gnodes = gdofs/ncompq; 1269 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1270 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1271 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1272 const char *usedresource; 1273 CeedGetResource(ceed, &usedresource); 1274 1275 ierr = PetscPrintf(comm, 1276 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1277 " rank(s) : %d\n" 1278 " Problem:\n" 1279 " Problem Name : %s\n" 1280 " Stabilization : %s\n" 1281 " PETSc:\n" 1282 " Box Faces : %s\n" 1283 " libCEED:\n" 1284 " libCEED Backend : %s\n" 1285 " libCEED Backend MemType : %s\n" 1286 " libCEED User Requested MemType : %s\n" 1287 " Mesh:\n" 1288 " Number of 1D Basis Nodes (P) : %d\n" 1289 " Number of 1D Quadrature Points (Q) : %d\n" 1290 " Global DoFs : %D\n" 1291 " Owned DoFs : %D\n" 1292 " DoFs per node : %D\n" 1293 " Global nodes : %D\n" 1294 " Owned nodes : %D\n", 1295 comm_size, problemTypes[problemChoice], 1296 StabilizationTypes[stab], box_faces_str, usedresource, 1297 CeedMemTypes[memtypebackend], 1298 (setmemtyperequest) ? 1299 CeedMemTypes[memtyperequested] : "none", 1300 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1301 CHKERRQ(ierr); 1302 } 1303 1304 // Set up global mass vector 1305 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1306 1307 // Set up libCEED 1308 // CEED Bases 1309 CeedInit(ceedresource, &ceed); 1310 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1311 &basisq); 1312 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1313 &basisx); 1314 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1315 CEED_GAUSS_LOBATTO, &basisxc); 1316 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1317 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1318 CHKERRQ(ierr); 1319 1320 // CEED Restrictions 1321 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1322 qdatasizeVol, &restrictq, &restrictx, 1323 &restrictqdi); CHKERRQ(ierr); 1324 1325 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1326 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1327 1328 // Create the CEED vectors that will be needed in setup 1329 CeedInt NqptsVol; 1330 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1331 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1332 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1333 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1334 1335 // Create the Q-function that builds the quadrature data for the NS operator 1336 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1337 &qf_setupVol); 1338 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1339 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1340 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1341 1342 // Create the Q-function that sets the ICs of the operator 1343 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1344 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1345 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1346 1347 qf_rhsVol = NULL; 1348 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1349 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1350 problem->applyVol_rhs_loc, &qf_rhsVol); 1351 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1352 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1353 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1354 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1355 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1356 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1357 } 1358 1359 qf_ifunctionVol = NULL; 1360 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1361 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1362 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1363 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1364 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1365 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1366 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1367 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1368 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1369 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1370 } 1371 1372 // Create the operator that builds the quadrature data for the NS operator 1373 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1374 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1375 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1376 basisx, CEED_VECTOR_NONE); 1377 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1378 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1379 1380 // Create the operator that sets the ICs 1381 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1382 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1383 CeedOperatorSetField(op_ics, "q0", restrictq, 1384 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1385 1386 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1387 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1388 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1389 1390 if (qf_rhsVol) { // Create the RHS physics operator 1391 CeedOperator op; 1392 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1393 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1394 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1395 CeedOperatorSetField(op, "qdata", restrictqdi, 1396 CEED_BASIS_COLLOCATED, qdata); 1397 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1398 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1399 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1400 user->op_rhs_vol = op; 1401 } 1402 1403 if (qf_ifunctionVol) { // Create the IFunction operator 1404 CeedOperator op; 1405 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1406 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1407 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1408 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1409 CeedOperatorSetField(op, "qdata", restrictqdi, 1410 CEED_BASIS_COLLOCATED, qdata); 1411 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1412 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1413 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1414 user->op_ifunction_vol = op; 1415 } 1416 1417 // Set up CEED for the boundaries 1418 CeedInt height = 1; 1419 CeedInt dimSur = dim - height; 1420 CeedInt numP_Sur = degree + 1; 1421 CeedInt numQ_Sur = numP_Sur + qextraSur; 1422 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1423 CeedBasis basisxSur, basisxcSur, basisqSur; 1424 CeedInt NqptsSur; 1425 CeedQFunction qf_setupSur, qf_applySur; 1426 1427 // CEED bases for the boundaries 1428 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, CEED_GAUSS, 1429 &basisqSur); 1430 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1431 &basisxSur); 1432 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1433 CEED_GAUSS_LOBATTO, &basisxcSur); 1434 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1435 1436 // Create the Q-function that builds the quadrature data for the Surface operator 1437 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1438 &qf_setupSur); 1439 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1440 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1441 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1442 1443 // Creat Q-Function for Boundaries 1444 qf_applySur = NULL; 1445 if (problem->applySur) { 1446 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1447 problem->applySur_loc, &qf_applySur); 1448 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1449 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1450 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1451 CeedQFunctionAddOutput(qf_applySur, "v", ncompq, CEED_EVAL_INTERP); 1452 } 1453 1454 // Create CEED Operator for the whole domain 1455 if (!implicit) 1456 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, user->op_rhs_vol, qf_applySur, qf_setupSur, 1457 height, numP_Sur, numQ_Sur, qdatasizeSur, NqptsSur, basisxSur, 1458 basisqSur, &user->op_rhs); CHKERRQ(ierr); 1459 if (implicit) 1460 ierr = CreateOperatorForDomain(ceed, dm, &bc, wind_type, user->op_ifunction_vol, qf_applySur, qf_setupSur, 1461 height, numP_Sur, numQ_Sur, qdatasizeSur, NqptsSur, basisxSur, 1462 basisqSur, &user->op_ifunction); CHKERRQ(ierr); 1463 // Set up contex for QFunctions 1464 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1465 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1466 struct Advection2dContext_ ctxAdvection2d = { 1467 .CtauS = CtauS, 1468 .strong_form = strong_form, 1469 .stabilization = stab, 1470 }; 1471 struct SurfaceContext_ ctxSurface = { 1472 .E_wind = E_wind, 1473 .strong_form = strong_form, 1474 .implicit = implicit, 1475 }; 1476 switch (problemChoice) { 1477 case NS_DENSITY_CURRENT: 1478 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1479 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1480 sizeof ctxNS); 1481 break; 1482 case NS_ADVECTION: 1483 case NS_ADVECTION2D: 1484 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1485 sizeof ctxAdvection2d); 1486 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1487 sizeof ctxAdvection2d); 1488 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &ctxSurface, sizeof ctxSurface); 1489 } 1490 1491 // Set up PETSc context 1492 // Set up units structure 1493 units->meter = meter; 1494 units->kilogram = kilogram; 1495 units->second = second; 1496 units->Kelvin = Kelvin; 1497 units->Pascal = Pascal; 1498 units->JperkgK = JperkgK; 1499 units->mpersquareds = mpersquareds; 1500 units->WpermK = WpermK; 1501 units->kgpercubicm = kgpercubicm; 1502 units->kgpersquaredms = kgpersquaredms; 1503 units->Joulepercubicm = Joulepercubicm; 1504 units->Joule = Joule; 1505 1506 // Set up user structure 1507 user->comm = comm; 1508 user->outputfreq = outputfreq; 1509 user->contsteps = contsteps; 1510 user->units = units; 1511 user->dm = dm; 1512 user->dmviz = dmviz; 1513 user->interpviz = interpviz; 1514 user->ceed = ceed; 1515 1516 // Calculate qdata and ICs 1517 // Set up state global and local vectors 1518 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1519 1520 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1521 1522 // Apply Setup Ceed Operators 1523 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1524 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1525 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1526 user->M); CHKERRQ(ierr); 1527 1528 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1529 &ctxSetup, 0.0); CHKERRQ(ierr); 1530 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1531 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1532 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1533 // slower execution. 1534 Vec Qbc; 1535 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1536 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1537 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1538 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1539 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1540 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1541 ierr = PetscObjectComposeFunction((PetscObject)dm, 1542 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1543 CHKERRQ(ierr); 1544 } 1545 1546 MPI_Comm_rank(comm, &rank); 1547 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1548 // Gather initial Q values 1549 // In case of continuation of simulation, set up initial values from binary file 1550 if (contsteps) { // continue from existent solution 1551 PetscViewer viewer; 1552 char filepath[PETSC_MAX_PATH_LEN]; 1553 // Read input 1554 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1555 user->outputfolder); 1556 CHKERRQ(ierr); 1557 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1558 CHKERRQ(ierr); 1559 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1560 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1561 } 1562 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1563 1564 // Create and setup TS 1565 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1566 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1567 if (implicit) { 1568 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1569 if (user->op_ifunction) { 1570 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1571 } else { // Implicit integrators can fall back to using an RHSFunction 1572 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1573 } 1574 } else { 1575 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1576 "Problem does not provide RHSFunction"); 1577 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1578 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1579 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1580 } 1581 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1582 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1583 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1584 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1585 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1586 ierr = TSAdaptSetStepLimits(adapt, 1587 1.e-12 * units->second, 1588 1.e2 * units->second); CHKERRQ(ierr); 1589 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1590 if (!contsteps) { // print initial condition 1591 if (testChoice == TEST_NONE) { 1592 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1593 } 1594 } else { // continue from time of last output 1595 PetscReal time; 1596 PetscInt count; 1597 PetscViewer viewer; 1598 char filepath[PETSC_MAX_PATH_LEN]; 1599 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1600 user->outputfolder); CHKERRQ(ierr); 1601 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1602 CHKERRQ(ierr); 1603 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1604 CHKERRQ(ierr); 1605 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1606 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1607 } 1608 if (testChoice == TEST_NONE) { 1609 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1610 } 1611 1612 // Solve 1613 start = MPI_Wtime(); 1614 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1615 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1616 cpu_time_used = MPI_Wtime() - start; 1617 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1618 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1619 comm); CHKERRQ(ierr); 1620 if (testChoice == TEST_NONE) { 1621 ierr = PetscPrintf(PETSC_COMM_WORLD, 1622 "Time taken for solution (sec): %g\n", 1623 (double)cpu_time_used); CHKERRQ(ierr); 1624 } 1625 1626 // Get error 1627 if (problem->non_zero_time && testChoice == TEST_NONE) { 1628 Vec Qexact, Qexactloc; 1629 PetscReal norm; 1630 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1631 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1632 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1633 1634 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1635 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1636 1637 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1638 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1639 CeedVectorDestroy(&q0ceed); 1640 ierr = PetscPrintf(PETSC_COMM_WORLD, 1641 "Max Error: %g\n", 1642 (double)norm); CHKERRQ(ierr); 1643 // Clean up vectors 1644 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1645 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1646 } 1647 1648 // Output Statistics 1649 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1650 if (testChoice == TEST_NONE) { 1651 ierr = PetscPrintf(PETSC_COMM_WORLD, 1652 "Time integrator took %D time steps to reach final time %g\n", 1653 steps, (double)ftime); CHKERRQ(ierr); 1654 } 1655 // Output numerical values from command line 1656 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1657 1658 // compare reference solution values with current run 1659 if (testChoice != TEST_NONE) { 1660 PetscViewer viewer; 1661 // Read reference file 1662 Vec Qref; 1663 PetscReal error, Qrefnorm; 1664 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1665 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1666 CHKERRQ(ierr); 1667 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1668 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1669 1670 // Compute error with respect to reference solution 1671 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1672 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1673 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1674 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1675 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1676 // Check error 1677 if (error > test->testtol) { 1678 ierr = PetscPrintf(PETSC_COMM_WORLD, 1679 "Test failed with error norm %g\n", 1680 (double)error); CHKERRQ(ierr); 1681 } 1682 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1683 } 1684 1685 // Clean up libCEED 1686 CeedVectorDestroy(&qdata); 1687 CeedVectorDestroy(&user->qceed); 1688 CeedVectorDestroy(&user->qdotceed); 1689 CeedVectorDestroy(&user->gceed); 1690 CeedVectorDestroy(&xcorners); 1691 CeedBasisDestroy(&basisq); 1692 CeedBasisDestroy(&basisx); 1693 CeedBasisDestroy(&basisxc); 1694 CeedElemRestrictionDestroy(&restrictq); 1695 CeedElemRestrictionDestroy(&restrictx); 1696 CeedElemRestrictionDestroy(&restrictqdi); 1697 CeedQFunctionDestroy(&qf_setupVol); 1698 CeedQFunctionDestroy(&qf_ics); 1699 CeedQFunctionDestroy(&qf_rhsVol); 1700 CeedQFunctionDestroy(&qf_ifunctionVol); 1701 CeedOperatorDestroy(&op_setupVol); 1702 CeedOperatorDestroy(&op_ics); 1703 CeedOperatorDestroy(&user->op_rhs_vol); 1704 CeedOperatorDestroy(&user->op_ifunction_vol); 1705 CeedDestroy(&ceed); 1706 CeedBasisDestroy(&basisqSur); 1707 CeedBasisDestroy(&basisxSur); 1708 CeedBasisDestroy(&basisxcSur); 1709 CeedQFunctionDestroy(&qf_setupSur); 1710 CeedQFunctionDestroy(&qf_applySur); 1711 CeedOperatorDestroy(&user->op_rhs); 1712 CeedOperatorDestroy(&user->op_ifunction); 1713 1714 // Clean up PETSc 1715 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1716 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1717 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1718 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1719 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1720 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1721 ierr = PetscFree(units); CHKERRQ(ierr); 1722 ierr = PetscFree(user); CHKERRQ(ierr); 1723 return PetscFinalize(); 1724 } 1725 1726