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