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