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