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