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