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; 867 PetscBool flg; 868 ierr = PetscOptionsIntArray("-bc_outflow", 869 "Use outflow boundary conditions on this list of faces", 870 NULL, bc.outflow, 871 (len = sizeof(bc.outflow) / sizeof(bc.outflow[0]), 872 &len), &flg); CHKERRQ(ierr); 873 if (flg) { 874 bc.noutflow = len; 875 // Using outflow boundaries disables automatic wall/slip boundaries (they must be set explicitly) 876 bc.nwall = 0; 877 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 878 } 879 ierr = PetscOptionsIntArray("-bc_wall", 880 "Use wall boundary conditions on this list of faces", 881 NULL, bc.walls, 882 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 883 &len), &flg); CHKERRQ(ierr); 884 if (flg) { 885 bc.nwall = len; 886 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 887 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 888 } 889 for (PetscInt j=0; j<3; j++) { 890 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 891 ierr = PetscOptionsIntArray(flags[j], 892 "Use slip boundary conditions on this list of faces", 893 NULL, bc.slips[j], 894 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 895 &len), &flg); 896 CHKERRQ(ierr); 897 if (flg) { 898 bc.nslip[j] = len; 899 bc.userbc = PETSC_TRUE; 900 } 901 } 902 } 903 ierr = PetscOptionsInt("-viz_refine", 904 "Regular refinement levels for visualization", 905 NULL, viz_refine, &viz_refine, NULL); 906 CHKERRQ(ierr); 907 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 908 NULL, meter, &meter, NULL); CHKERRQ(ierr); 909 meter = fabs(meter); 910 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 911 NULL, second, &second, NULL); CHKERRQ(ierr); 912 second = fabs(second); 913 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 914 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 915 kilogram = fabs(kilogram); 916 ierr = PetscOptionsScalar("-units_Kelvin", 917 "1 Kelvin in scaled temperature units", 918 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 919 Kelvin = fabs(Kelvin); 920 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 921 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 922 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 923 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 924 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 925 NULL, P0, &P0, NULL); CHKERRQ(ierr); 926 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 927 NULL, N, &N, NULL); CHKERRQ(ierr); 928 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 929 NULL, cv, &cv, NULL); CHKERRQ(ierr); 930 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 931 NULL, cp, &cp, NULL); CHKERRQ(ierr); 932 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 933 NULL, g, &g, NULL); CHKERRQ(ierr); 934 ierr = PetscOptionsScalar("-lambda", 935 "Stokes hypothesis second viscosity coefficient", 936 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 937 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 938 NULL, mu, &mu, NULL); CHKERRQ(ierr); 939 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 940 NULL, k, &k, NULL); CHKERRQ(ierr); 941 ierr = PetscOptionsScalar("-CtauS", 942 "Scale coefficient for tau (nondimensional)", 943 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 944 if (stab == STAB_NONE && CtauS != 0) { 945 ierr = PetscPrintf(comm, 946 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 947 CHKERRQ(ierr); 948 } 949 ierr = PetscOptionsScalar("-strong_form", 950 "Strong (1) or weak/integrated by parts (0) advection residual", 951 NULL, strong_form, &strong_form, NULL); 952 CHKERRQ(ierr); 953 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 954 ierr = PetscPrintf(comm, 955 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 956 CHKERRQ(ierr); 957 } 958 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 959 NULL, lx, &lx, NULL); CHKERRQ(ierr); 960 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 961 NULL, ly, &ly, NULL); CHKERRQ(ierr); 962 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 963 NULL, lz, &lz, NULL); CHKERRQ(ierr); 964 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 965 NULL, rc, &rc, NULL); CHKERRQ(ierr); 966 ierr = PetscOptionsScalar("-resx","Target resolution in x", 967 NULL, resx, &resx, NULL); CHKERRQ(ierr); 968 ierr = PetscOptionsScalar("-resy","Target resolution in y", 969 NULL, resy, &resy, NULL); CHKERRQ(ierr); 970 ierr = PetscOptionsScalar("-resz","Target resolution in z", 971 NULL, resz, &resz, NULL); CHKERRQ(ierr); 972 PetscInt n = problem->dim; 973 center[0] = 0.5 * lx; 974 center[1] = 0.5 * ly; 975 center[2] = 0.5 * lz; 976 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 977 NULL, center, &n, NULL); CHKERRQ(ierr); 978 n = problem->dim; 979 ierr = PetscOptionsRealArray("-dc_axis", 980 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 981 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 982 { 983 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 984 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 985 if (norm > 0) { 986 for (int i=0; i<3; i++) dc_axis[i] /= norm; 987 } 988 } 989 ierr = PetscOptionsInt("-output_freq", 990 "Frequency of output, in number of steps", 991 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 992 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 993 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 994 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 995 NULL, degree, °ree, NULL); CHKERRQ(ierr); 996 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 997 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 998 ierr = PetscStrncpy(user->outputfolder, ".", 2); CHKERRQ(ierr); 999 ierr = PetscOptionsString("-of", "Output folder", 1000 NULL, user->outputfolder, user->outputfolder, 1001 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 1002 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1003 ierr = PetscOptionsEnum("-memtype", 1004 "CEED MemType requested", NULL, 1005 memTypes, (PetscEnum)memtyperequested, 1006 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1007 CHKERRQ(ierr); 1008 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1009 1010 // Define derived units 1011 Pascal = kilogram / (meter * PetscSqr(second)); 1012 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1013 mpersquareds = meter / PetscSqr(second); 1014 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1015 kgpercubicm = kilogram / pow(meter,3); 1016 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1017 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1018 1019 // Scale variables to desired units 1020 theta0 *= Kelvin; 1021 thetaC *= Kelvin; 1022 P0 *= Pascal; 1023 N *= (1./second); 1024 cv *= JperkgK; 1025 cp *= JperkgK; 1026 Rd = cp - cv; 1027 g *= mpersquareds; 1028 mu *= Pascal * second; 1029 k *= WpermK; 1030 lx = fabs(lx) * meter; 1031 ly = fabs(ly) * meter; 1032 lz = fabs(lz) * meter; 1033 rc = fabs(rc) * meter; 1034 resx = fabs(resx) * meter; 1035 resy = fabs(resy) * meter; 1036 resz = fabs(resz) * meter; 1037 for (int i=0; i<3; i++) center[i] *= meter; 1038 1039 const CeedInt dim = problem->dim, ncompx = problem->dim, 1040 qdatasizeVol = problem->qdatasizeVol; 1041 // Set up the libCEED context 1042 struct SetupContext_ ctxSetup = { 1043 .theta0 = theta0, 1044 .thetaC = thetaC, 1045 .P0 = P0, 1046 .N = N, 1047 .cv = cv, 1048 .cp = cp, 1049 .Rd = Rd, 1050 .g = g, 1051 .rc = rc, 1052 .lx = lx, 1053 .ly = ly, 1054 .lz = lz, 1055 .center[0] = center[0], 1056 .center[1] = center[1], 1057 .center[2] = center[2], 1058 .dc_axis[0] = dc_axis[0], 1059 .dc_axis[1] = dc_axis[1], 1060 .dc_axis[2] = dc_axis[2], 1061 .time = 0, 1062 }; 1063 1064 // Create the mesh 1065 { 1066 const PetscReal scale[3] = {lx, ly, lz}; 1067 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1068 NULL, PETSC_TRUE, &dm); 1069 CHKERRQ(ierr); 1070 } 1071 1072 // Distribute the mesh over processes 1073 { 1074 DM dmDist = NULL; 1075 PetscPartitioner part; 1076 1077 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1078 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1079 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1080 if (dmDist) { 1081 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1082 dm = dmDist; 1083 } 1084 } 1085 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1086 1087 // Setup DM 1088 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1089 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1090 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetup); CHKERRQ(ierr); 1091 1092 // Refine DM for high-order viz 1093 dmviz = NULL; 1094 interpviz = NULL; 1095 if (viz_refine) { 1096 DM dmhierarchy[viz_refine+1]; 1097 1098 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1099 dmhierarchy[0] = dm; 1100 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1101 Mat interp_next; 1102 1103 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1104 CHKERRQ(ierr); 1105 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1106 d = (d + 1) / 2; 1107 if (i + 1 == viz_refine) d = 1; 1108 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 1109 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1110 &interp_next, NULL); CHKERRQ(ierr); 1111 if (!i) interpviz = interp_next; 1112 else { 1113 Mat C; 1114 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1115 PETSC_DECIDE, &C); CHKERRQ(ierr); 1116 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1117 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1118 interpviz = C; 1119 } 1120 } 1121 for (PetscInt i=1; i<viz_refine; i++) { 1122 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1123 } 1124 dmviz = dmhierarchy[viz_refine]; 1125 } 1126 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1127 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1128 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1129 lnodes /= ncompq; 1130 1131 // Initialize CEED 1132 CeedInit(ceedresource, &ceed); 1133 // Set memtype 1134 CeedMemType memtypebackend; 1135 CeedGetPreferredMemType(ceed, &memtypebackend); 1136 // Check memtype compatibility 1137 if (!setmemtyperequest) 1138 memtyperequested = memtypebackend; 1139 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1140 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1141 "PETSc was not built with CUDA. " 1142 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1143 1144 // Set number of 1D nodes and quadrature points 1145 numP = degree + 1; 1146 numQ = numP + qextra; 1147 1148 // Print summary 1149 if (testChoice == TEST_NONE) { 1150 CeedInt gdofs, odofs; 1151 int comm_size; 1152 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1153 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1154 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1155 gnodes = gdofs/ncompq; 1156 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1157 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1158 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1159 const char *usedresource; 1160 CeedGetResource(ceed, &usedresource); 1161 1162 ierr = PetscPrintf(comm, 1163 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1164 " rank(s) : %d\n" 1165 " Problem:\n" 1166 " Problem Name : %s\n" 1167 " Stabilization : %s\n" 1168 " PETSc:\n" 1169 " Box Faces : %s\n" 1170 " libCEED:\n" 1171 " libCEED Backend : %s\n" 1172 " libCEED Backend MemType : %s\n" 1173 " libCEED User Requested MemType : %s\n" 1174 " Mesh:\n" 1175 " Number of 1D Basis Nodes (P) : %d\n" 1176 " Number of 1D Quadrature Points (Q) : %d\n" 1177 " Global DoFs : %D\n" 1178 " Owned DoFs : %D\n" 1179 " DoFs per node : %D\n" 1180 " Global nodes : %D\n" 1181 " Owned nodes : %D\n", 1182 comm_size, problemTypes[problemChoice], 1183 StabilizationTypes[stab], box_faces_str, usedresource, 1184 CeedMemTypes[memtypebackend], 1185 (setmemtyperequest) ? 1186 CeedMemTypes[memtyperequested] : "none", 1187 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1188 CHKERRQ(ierr); 1189 } 1190 1191 // Set up global mass vector 1192 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1193 1194 // Set up libCEED 1195 // CEED Bases 1196 CeedInit(ceedresource, &ceed); 1197 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1198 &basisq); 1199 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1200 &basisx); 1201 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1202 CEED_GAUSS_LOBATTO, &basisxc); 1203 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1204 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1205 CHKERRQ(ierr); 1206 1207 // CEED Restrictions 1208 ierr = GetRestrictionForDomain(ceed, dm, ncompx, dim, 0, 0, 0, numP, numQ, 1209 qdatasizeVol, &restrictq, &restrictx, 1210 &restrictqdi); CHKERRQ(ierr); 1211 1212 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1213 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1214 1215 // Create the CEED vectors that will be needed in setup 1216 CeedInt NqptsVol; 1217 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1218 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1219 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1220 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1221 1222 // Create the Q-function that builds the quadrature data for the NS operator 1223 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1224 &qf_setupVol); 1225 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1226 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1227 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1228 1229 // Create the Q-function that sets the ICs of the operator 1230 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1231 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1232 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1233 1234 qf_rhsVol = NULL; 1235 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1236 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1237 problem->applyVol_rhs_loc, &qf_rhsVol); 1238 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1239 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1240 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1241 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1242 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1243 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1244 } 1245 1246 qf_ifunctionVol = NULL; 1247 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1248 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1249 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1250 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1251 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1252 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1253 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1254 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1255 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1256 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1257 } 1258 1259 // Create the operator that builds the quadrature data for the NS operator 1260 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1261 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1262 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1263 basisx, CEED_VECTOR_NONE); 1264 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1265 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1266 1267 // Create the operator that sets the ICs 1268 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1269 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1270 CeedOperatorSetField(op_ics, "q0", restrictq, 1271 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1272 1273 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1274 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1275 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1276 1277 if (qf_rhsVol) { // Create the RHS physics operator 1278 CeedOperator op; 1279 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1280 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1281 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1282 CeedOperatorSetField(op, "qdata", restrictqdi, 1283 CEED_BASIS_COLLOCATED, qdata); 1284 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1285 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1286 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1287 user->op_rhs = op; 1288 } 1289 1290 if (qf_ifunctionVol) { // Create the IFunction operator 1291 CeedOperator op; 1292 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1293 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1294 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1295 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1296 CeedOperatorSetField(op, "qdata", restrictqdi, 1297 CEED_BASIS_COLLOCATED, qdata); 1298 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1299 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1300 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1301 user->op_ifunction = op; 1302 } 1303 1304 //--------------------------------------------------------------------------------------// 1305 // Outflow Boundary Condition 1306 //--------------------------------------------------------------------------------------// 1307 // Set up CEED for the boundaries 1308 CeedInt numP_Sur, numQ_Sur; 1309 CeedInt height = 1; 1310 CeedInt dimSur = dim - height; 1311 numP_Sur = degreeSur + 1; 1312 numQ_Sur = numP_Sur + qextraSur; 1313 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1314 CeedBasis basisxSur, basisxcSur, basisqSur; 1315 CeedElemRestriction restrictxSur[6], restrictqSur[6], restrictqdiSur[6]; 1316 CeedInt NqptsSur; 1317 PetscInt localNelemSur[6]; 1318 CeedVector qdataSur[6], qdataSur_[6]; 1319 CeedQFunction qf_setupSur, qf_rhsSur, qf_ifunctionSur; 1320 CeedOperator op_setupSur[6], op_setupSur_[6]; 1321 PetscInt numOutFlow = bc.noutflow; 1322 DMLabel domainLabel; 1323 1324 // Get Label for the boundaries 1325 ierr = DMGetLabel(dm, "Face Sets", &domainLabel); CHKERRQ(ierr); 1326 1327 // CEED bases for the boundaries 1328 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, CEED_GAUSS, 1329 &basisqSur); 1330 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1331 &basisxSur); 1332 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1333 CEED_GAUSS_LOBATTO, &basisxcSur); 1334 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1335 1336 // Create the Q-function that builds the quadrature data for the Surface operator 1337 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1338 &qf_setupSur); 1339 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1340 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1341 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1342 1343 qf_rhsSur = NULL; 1344 if (problem->applySur_rhs) { // Create the Q-function that defines the action of the RHS operator on the Surface 1345 CeedQFunctionCreateInterior(ceed, 1, problem->applySur_rhs, 1346 problem->applySur_rhs_loc, &qf_rhsSur); 1347 CeedQFunctionAddInput(qf_rhsSur, "q", ncompq, CEED_EVAL_INTERP); 1348 CeedQFunctionAddInput(qf_rhsSur, "dq", ncompq*dim, CEED_EVAL_GRAD); //assumed volumetric elements 1349 CeedQFunctionAddInput(qf_rhsSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1350 CeedQFunctionAddInput(qf_rhsSur, "x", ncompx, CEED_EVAL_INTERP); 1351 CeedQFunctionAddInput(qf_rhsSur, "qdata", qdatasizeVol, CEED_EVAL_NONE); //assumed volumetric elements 1352 CeedQFunctionAddOutput(qf_rhsSur, "v", ncompq, CEED_EVAL_INTERP); 1353 } 1354 qf_ifunctionSur = NULL; 1355 if (problem->applySur_ifunction) { // Create the Q-function that defines the action of the IFunction 1356 CeedQFunctionCreateInterior(ceed, 1, problem->applySur_ifunction, 1357 problem->applySur_ifunction_loc, &qf_ifunctionSur); 1358 CeedQFunctionAddInput(qf_ifunctionSur, "q", ncompq, CEED_EVAL_INTERP); 1359 CeedQFunctionAddInput(qf_ifunctionSur, "dq", ncompq*dim, CEED_EVAL_GRAD); //assumed volumetric elements 1360 CeedQFunctionAddInput(qf_ifunctionSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1361 CeedQFunctionAddInput(qf_ifunctionSur, "x", ncompx, CEED_EVAL_INTERP); 1362 CeedQFunctionAddInput(qf_ifunctionSur, "qdata", qdatasizeVol, CEED_EVAL_NONE); //assumed volumetric elements 1363 CeedQFunctionAddOutput(qf_ifunctionSur, "v", ncompq, CEED_EVAL_INTERP); 1364 } 1365 // Create CEED Operator for each face 1366 for(CeedInt i=0; i<numOutFlow; i++){ 1367 ierr = GetRestrictionForDomain(ceed, dm, ncompx, dimSur, height, domainLabel, bc.outflow[i], numP_Sur, 1368 numQ_Sur, qdatasizeSur, &restrictqSur[i], &restrictxSur[i], 1369 &restrictqdiSur[i]); CHKERRQ(ierr); 1370 // Create the CEED vectors that will be needed in setup 1371 CeedElemRestrictionGetNumElements(restrictqSur[i], &localNelemSur[i]); 1372 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, &qdataSur[i]); 1373 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, &qdataSur_[i]); 1374 1375 // Create the operator that builds the quadrature data for the NS operator 1376 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur[i]); 1377 CeedOperatorSetField(op_setupSur[i], "dx", restrictxSur[i], basisxSur, CEED_VECTOR_ACTIVE); 1378 CeedOperatorSetField(op_setupSur[i], "weight", CEED_ELEMRESTRICTION_NONE, 1379 basisxSur, CEED_VECTOR_NONE); 1380 CeedOperatorSetField(op_setupSur[i], "qdataSur", restrictqdiSur[i], 1381 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1382 1383 // Utility operator that builds the quadrature data for computing viscous terms 1384 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupSur_[i]); 1385 CeedOperatorSetField(op_setupSur_[i], "dx", restrictxSur[i], basisx, CEED_VECTOR_ACTIVE); 1386 CeedOperatorSetField(op_setupSur_[i], "weight", CEED_ELEMRESTRICTION_NONE, 1387 basisx, CEED_VECTOR_NONE); 1388 CeedOperatorSetField(op_setupSur_[i], "qdata", restrictqdiSur[i], 1389 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1390 1391 if (qf_rhsSur) { // Create the RHS physics operator 1392 CeedOperator op; 1393 CeedOperatorCreate(ceed, qf_rhsSur, NULL, NULL, &op); 1394 CeedOperatorSetField(op, "q", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1395 CeedOperatorSetField(op, "dq", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1396 CeedOperatorSetField(op, "qdataSur", restrictqdiSur[i], 1397 CEED_BASIS_COLLOCATED, qdataSur[i]); 1398 CeedOperatorSetField(op, "x", restrictxSur[i], basisxSur, xcorners); 1399 CeedOperatorSetField(op, "qdata", restrictqdiSur[i], 1400 CEED_BASIS_COLLOCATED, qdataSur_[i]); 1401 CeedOperatorSetField(op, "v", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1402 user->op_rhs_sur[i] = op; 1403 } 1404 1405 if (qf_ifunctionSur) { // Create the IFunction operator 1406 CeedOperator op; 1407 CeedOperatorCreate(ceed, qf_ifunctionSur, NULL, NULL, &op); 1408 CeedOperatorSetField(op, "q", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1409 CeedOperatorSetField(op, "dq", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1410 CeedOperatorSetField(op, "qdataSur", restrictqdiSur[i], 1411 CEED_BASIS_COLLOCATED, qdataSur[i]); 1412 CeedOperatorSetField(op, "x", restrictxSur[i], basisxSur, xcorners); 1413 CeedOperatorSetField(op, "qdata", restrictqdiSur[i], 1414 CEED_BASIS_COLLOCATED, qdataSur_[i]); 1415 CeedOperatorSetField(op, "v", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1416 user->op_ifunction_sur[i] = op; 1417 } 1418 } 1419 // Composite Operaters 1420 // IFunction 1421 if (user->op_ifunction_vol) { 1422 if (numOutFlow>0) { 1423 // Composite Operators for the IFunction 1424 CeedCompositeOperatorCreate(ceed, &user->op_ifunction); 1425 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_vol); 1426 for(CeedInt i=0; i<numOutFlow; i++){ 1427 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_sur[i]); 1428 } 1429 } else { 1430 user->op_ifunction = user->op_ifunction_vol; 1431 } 1432 } 1433 // RHS 1434 if (user->op_rhs_vol) { 1435 if (numOutFlow == 1) { 1436 // Composite Operators for the RHS 1437 CeedCompositeOperatorCreate(ceed, &user->op_rhs); 1438 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_vol); 1439 for(CeedInt i=0; i<numOutFlow; i++){ 1440 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_sur[i]); 1441 } 1442 } else { 1443 user->op_rhs = user->op_rhs_vol; 1444 } 1445 } 1446 //--------------------------------------------------------------------------------------// 1447 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1448 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1449 struct Advection2dContext_ ctxAdvection2d = { 1450 .CtauS = CtauS, 1451 .strong_form = strong_form, 1452 .stabilization = stab, 1453 }; 1454 switch (problemChoice) { 1455 case NS_DENSITY_CURRENT: 1456 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1457 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1458 sizeof ctxNS); 1459 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxNS, sizeof ctxNS); 1460 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxNS, 1461 sizeof ctxNS); 1462 break; 1463 case NS_ADVECTION: 1464 case NS_ADVECTION2D: 1465 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1466 sizeof ctxAdvection2d); 1467 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1468 sizeof ctxAdvection2d); 1469 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxAdvection2d, 1470 sizeof ctxAdvection2d); 1471 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxAdvection2d, 1472 sizeof ctxAdvection2d); 1473 } 1474 1475 // Set up PETSc context 1476 // Set up units structure 1477 units->meter = meter; 1478 units->kilogram = kilogram; 1479 units->second = second; 1480 units->Kelvin = Kelvin; 1481 units->Pascal = Pascal; 1482 units->JperkgK = JperkgK; 1483 units->mpersquareds = mpersquareds; 1484 units->WpermK = WpermK; 1485 units->kgpercubicm = kgpercubicm; 1486 units->kgpersquaredms = kgpersquaredms; 1487 units->Joulepercubicm = Joulepercubicm; 1488 1489 // Set up user structure 1490 user->comm = comm; 1491 user->outputfreq = outputfreq; 1492 user->contsteps = contsteps; 1493 user->units = units; 1494 user->dm = dm; 1495 user->dmviz = dmviz; 1496 user->interpviz = interpviz; 1497 user->ceed = ceed; 1498 1499 // Calculate qdata and ICs 1500 // Set up state global and local vectors 1501 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1502 1503 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1504 1505 // Apply Setup Ceed Operators 1506 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1507 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1508 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1509 user->M); CHKERRQ(ierr); 1510 1511 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1512 &ctxSetup, 0.0); CHKERRQ(ierr); 1513 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1514 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1515 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1516 // slower execution. 1517 Vec Qbc; 1518 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1519 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1520 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1521 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1522 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1523 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1524 ierr = PetscObjectComposeFunction((PetscObject)dm, 1525 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1526 CHKERRQ(ierr); 1527 } 1528 1529 MPI_Comm_rank(comm, &rank); 1530 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1531 // Gather initial Q values 1532 // In case of continuation of simulation, set up initial values from binary file 1533 if (contsteps) { // continue from existent solution 1534 PetscViewer viewer; 1535 char filepath[PETSC_MAX_PATH_LEN]; 1536 // Read input 1537 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1538 user->outputfolder); 1539 CHKERRQ(ierr); 1540 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1541 CHKERRQ(ierr); 1542 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1543 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1544 } 1545 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1546 1547 // Create and setup TS 1548 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1549 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1550 if (implicit) { 1551 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1552 if (user->op_ifunction) { 1553 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1554 } else { // Implicit integrators can fall back to using an RHSFunction 1555 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1556 } 1557 } else { 1558 if (!user->op_rhs) SETERRQ(comm,PETSC_ERR_ARG_NULL, 1559 "Problem does not provide RHSFunction"); 1560 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1561 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1562 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1563 } 1564 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1565 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1566 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1567 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1568 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1569 ierr = TSAdaptSetStepLimits(adapt, 1570 1.e-12 * units->second, 1571 1.e2 * units->second); CHKERRQ(ierr); 1572 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1573 if (!contsteps) { // print initial condition 1574 if (testChoice == TEST_NONE) { 1575 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1576 } 1577 } else { // continue from time of last output 1578 PetscReal time; 1579 PetscInt count; 1580 PetscViewer viewer; 1581 char filepath[PETSC_MAX_PATH_LEN]; 1582 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1583 user->outputfolder); CHKERRQ(ierr); 1584 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1585 CHKERRQ(ierr); 1586 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1587 CHKERRQ(ierr); 1588 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1589 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1590 } 1591 if (testChoice == TEST_NONE) { 1592 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1593 } 1594 1595 // Solve 1596 start = MPI_Wtime(); 1597 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1598 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1599 cpu_time_used = MPI_Wtime() - start; 1600 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1601 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1602 comm); CHKERRQ(ierr); 1603 if (testChoice == TEST_NONE) { 1604 ierr = PetscPrintf(PETSC_COMM_WORLD, 1605 "Time taken for solution (sec): %g\n", 1606 (double)cpu_time_used); CHKERRQ(ierr); 1607 } 1608 1609 // Get error 1610 if (problem->non_zero_time && testChoice == TEST_NONE) { 1611 Vec Qexact, Qexactloc; 1612 PetscReal norm; 1613 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1614 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1615 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1616 1617 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1618 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1619 1620 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1621 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1622 CeedVectorDestroy(&q0ceed); 1623 ierr = PetscPrintf(PETSC_COMM_WORLD, 1624 "Max Error: %g\n", 1625 (double)norm); CHKERRQ(ierr); 1626 // Clean up vectors 1627 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1628 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1629 } 1630 1631 // Output Statistics 1632 ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr); 1633 if (testChoice == TEST_NONE) { 1634 ierr = PetscPrintf(PETSC_COMM_WORLD, 1635 "Time integrator took %D time steps to reach final time %g\n", 1636 steps, (double)ftime); CHKERRQ(ierr); 1637 } 1638 // Output numerical values from command line 1639 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1640 1641 // compare reference solution values with current run 1642 if (testChoice != TEST_NONE) { 1643 PetscViewer viewer; 1644 // Read reference file 1645 Vec Qref; 1646 PetscReal error, Qrefnorm; 1647 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1648 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1649 CHKERRQ(ierr); 1650 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1651 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1652 1653 // Compute error with respect to reference solution 1654 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1655 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1656 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1657 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1658 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1659 // Check error 1660 if (error > test->testtol) { 1661 ierr = PetscPrintf(PETSC_COMM_WORLD, 1662 "Test failed with error norm %g\n", 1663 (double)error); CHKERRQ(ierr); 1664 } 1665 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1666 } 1667 1668 // Clean up libCEED 1669 CeedVectorDestroy(&qdata); 1670 CeedVectorDestroy(&user->qceed); 1671 CeedVectorDestroy(&user->qdotceed); 1672 CeedVectorDestroy(&user->gceed); 1673 CeedVectorDestroy(&xcorners); 1674 CeedBasisDestroy(&basisq); 1675 CeedBasisDestroy(&basisx); 1676 CeedBasisDestroy(&basisxc); 1677 CeedElemRestrictionDestroy(&restrictq); 1678 CeedElemRestrictionDestroy(&restrictx); 1679 CeedElemRestrictionDestroy(&restrictqdi); 1680 CeedQFunctionDestroy(&qf_setupVol); 1681 CeedQFunctionDestroy(&qf_ics); 1682 CeedQFunctionDestroy(&qf_rhsVol); 1683 CeedQFunctionDestroy(&qf_ifunctionVol); 1684 CeedOperatorDestroy(&op_setupVol); 1685 CeedOperatorDestroy(&op_ics); 1686 CeedOperatorDestroy(&user->op_rhs_vol); 1687 CeedOperatorDestroy(&user->op_ifunction_vol); 1688 CeedDestroy(&ceed); 1689 CeedBasisDestroy(&basisqSur); 1690 CeedBasisDestroy(&basisxSur); 1691 CeedBasisDestroy(&basisxcSur); 1692 CeedQFunctionDestroy(&qf_setupSur); 1693 CeedQFunctionDestroy(&qf_rhsSur); 1694 CeedQFunctionDestroy(&qf_ifunctionSur); 1695 for(CeedInt i=0; i<numOutFlow; i++){ 1696 CeedVectorDestroy(&qdataSur[i]); 1697 CeedVectorDestroy(&qdataSur_[i]); 1698 CeedElemRestrictionDestroy(&restrictqSur[i]); 1699 CeedElemRestrictionDestroy(&restrictxSur[i]); 1700 CeedElemRestrictionDestroy(&restrictqdiSur[i]); 1701 CeedOperatorDestroy(&op_setupSur[i]); 1702 CeedOperatorDestroy(&op_setupSur_[i]); 1703 CeedOperatorDestroy(&user->op_rhs_sur[i]); 1704 CeedOperatorDestroy(&user->op_ifunction_sur[i]); 1705 } 1706 CeedOperatorDestroy(&user->op_rhs); 1707 CeedOperatorDestroy(&user->op_ifunction); 1708 1709 // Clean up PETSc 1710 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1711 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1712 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1713 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1714 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1715 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1716 ierr = PetscFree(units); CHKERRQ(ierr); 1717 ierr = PetscFree(user); CHKERRQ(ierr); 1718 return PetscFinalize(); 1719 } 1720 1721