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