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]; 1319 CeedQFunction qf_setupSur, qf_rhsSur, qf_ifunctionSur; 1320 CeedOperator 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, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1349 CeedQFunctionAddInput(qf_rhsSur, "x", ncompx, CEED_EVAL_INTERP); 1350 CeedQFunctionAddOutput(qf_rhsSur, "v", ncompq, CEED_EVAL_INTERP); 1351 } 1352 qf_ifunctionSur = NULL; 1353 if (problem->applySur_ifunction) { // Create the Q-function that defines the action of the IFunction 1354 CeedQFunctionCreateInterior(ceed, 1, problem->applySur_ifunction, 1355 problem->applySur_ifunction_loc, &qf_ifunctionSur); 1356 CeedQFunctionAddInput(qf_ifunctionSur, "q", ncompq, CEED_EVAL_INTERP); 1357 CeedQFunctionAddInput(qf_ifunctionSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1358 CeedQFunctionAddInput(qf_ifunctionSur, "x", ncompx, CEED_EVAL_INTERP); 1359 CeedQFunctionAddOutput(qf_ifunctionSur, "v", ncompq, CEED_EVAL_INTERP); 1360 } 1361 // Create CEED Operator for each face 1362 for(CeedInt i=0; i<numOutFlow; i++){ 1363 ierr = GetRestrictionForDomain(ceed, dm, ncompx, dimSur, height, domainLabel, bc.outflow[i], numP_Sur, 1364 numQ_Sur, qdatasizeSur, &restrictqSur[i], &restrictxSur[i], 1365 &restrictqdiSur[i]); CHKERRQ(ierr); 1366 // Create the CEED vectors that will be needed in setup 1367 CeedElemRestrictionGetNumElements(restrictqSur[i], &localNelemSur[i]); 1368 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, &qdataSur[i]); 1369 1370 // Create the operator that builds the quadrature data for the NS operator 1371 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur[i]); 1372 CeedOperatorSetField(op_setupSur[i], "dx", restrictxSur[i], basisxSur, CEED_VECTOR_ACTIVE); 1373 CeedOperatorSetField(op_setupSur[i], "weight", CEED_ELEMRESTRICTION_NONE, 1374 basisxSur, CEED_VECTOR_NONE); 1375 CeedOperatorSetField(op_setupSur[i], "qdataSur", restrictqdiSur[i], 1376 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1377 1378 if (qf_rhsSur) { // Create the RHS physics operator 1379 CeedOperator op; 1380 CeedOperatorCreate(ceed, qf_rhsSur, NULL, NULL, &op); 1381 CeedOperatorSetField(op, "q", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1382 CeedOperatorSetField(op, "qdataSur", restrictqdiSur[i], 1383 CEED_BASIS_COLLOCATED, qdataSur[i]); 1384 CeedOperatorSetField(op, "x", restrictxSur[i], basisxSur, xcorners); 1385 CeedOperatorSetField(op, "v", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1386 user->op_rhs_sur[i] = op; 1387 } 1388 1389 if (qf_ifunctionSur) { // Create the IFunction operator 1390 CeedOperator op; 1391 CeedOperatorCreate(ceed, qf_ifunctionSur, NULL, NULL, &op); 1392 CeedOperatorSetField(op, "q", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1393 CeedOperatorSetField(op, "qdataSur", restrictqdiSur[i], 1394 CEED_BASIS_COLLOCATED, qdataSur[i]); 1395 CeedOperatorSetField(op, "x", restrictxSur[i], basisxSur, xcorners); 1396 CeedOperatorSetField(op, "v", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1397 user->op_ifunction_sur[i] = op; 1398 } 1399 } 1400 // Composite Operaters 1401 // IFunction 1402 if (user->op_ifunction_vol) { 1403 if (numOutFlow>0) { 1404 // Composite Operators for the IFunction 1405 CeedCompositeOperatorCreate(ceed, &user->op_ifunction); 1406 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_vol); 1407 for(CeedInt i=0; i<numOutFlow; i++){ 1408 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_sur[i]); 1409 } 1410 } else { 1411 user->op_ifunction = user->op_ifunction_vol; 1412 } 1413 } 1414 // RHS 1415 if (user->op_rhs_vol) { 1416 if (numOutFlow == 1) { 1417 // Composite Operators for the RHS 1418 CeedCompositeOperatorCreate(ceed, &user->op_rhs); 1419 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_vol); 1420 for(CeedInt i=0; i<numOutFlow; i++){ 1421 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_sur[i]); 1422 } 1423 } else { 1424 user->op_rhs = user->op_rhs_vol; 1425 } 1426 } 1427 //--------------------------------------------------------------------------------------// 1428 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1429 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1430 struct Advection2dContext_ ctxAdvection2d = { 1431 .CtauS = CtauS, 1432 .strong_form = strong_form, 1433 .stabilization = stab, 1434 }; 1435 switch (problemChoice) { 1436 case NS_DENSITY_CURRENT: 1437 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1438 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1439 sizeof ctxNS); 1440 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxNS, sizeof ctxNS); 1441 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxNS, 1442 sizeof ctxNS); 1443 break; 1444 case NS_ADVECTION: 1445 case NS_ADVECTION2D: 1446 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1447 sizeof ctxAdvection2d); 1448 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1449 sizeof ctxAdvection2d); 1450 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxAdvection2d, 1451 sizeof ctxAdvection2d); 1452 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxAdvection2d, 1453 sizeof ctxAdvection2d); 1454 } 1455 1456 // Set up PETSc context 1457 // Set up units structure 1458 units->meter = meter; 1459 units->kilogram = kilogram; 1460 units->second = second; 1461 units->Kelvin = Kelvin; 1462 units->Pascal = Pascal; 1463 units->JperkgK = JperkgK; 1464 units->mpersquareds = mpersquareds; 1465 units->WpermK = WpermK; 1466 units->kgpercubicm = kgpercubicm; 1467 units->kgpersquaredms = kgpersquaredms; 1468 units->Joulepercubicm = Joulepercubicm; 1469 1470 // Set up user structure 1471 user->comm = comm; 1472 user->outputfreq = outputfreq; 1473 user->contsteps = contsteps; 1474 user->units = units; 1475 user->dm = dm; 1476 user->dmviz = dmviz; 1477 user->interpviz = interpviz; 1478 user->ceed = ceed; 1479 1480 // Calculate qdata and ICs 1481 // Set up state global and local vectors 1482 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1483 1484 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1485 1486 // Apply Setup Ceed Operators 1487 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1488 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1489 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1490 user->M); CHKERRQ(ierr); 1491 1492 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1493 &ctxSetup, 0.0); CHKERRQ(ierr); 1494 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1495 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1496 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1497 // slower execution. 1498 Vec Qbc; 1499 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1500 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1501 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1502 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1503 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1504 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1505 ierr = PetscObjectComposeFunction((PetscObject)dm, 1506 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1507 CHKERRQ(ierr); 1508 } 1509 1510 MPI_Comm_rank(comm, &rank); 1511 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1512 // Gather initial Q values 1513 // In case of continuation of simulation, set up initial values from binary file 1514 if (contsteps) { // continue from existent solution 1515 PetscViewer viewer; 1516 char filepath[PETSC_MAX_PATH_LEN]; 1517 // Read input 1518 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1519 user->outputfolder); 1520 CHKERRQ(ierr); 1521 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1522 CHKERRQ(ierr); 1523 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1524 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1525 } 1526 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1527 1528 // Create and setup TS 1529 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1530 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1531 if (implicit) { 1532 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1533 if (user->op_ifunction) { 1534 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1535 } else { // Implicit integrators can fall back to using an RHSFunction 1536 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1537 } 1538 } else { 1539 if (!user->op_rhs) SETERRQ(comm,PETSC_ERR_ARG_NULL, 1540 "Problem does not provide RHSFunction"); 1541 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1542 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1543 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1544 } 1545 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1546 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1547 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1548 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1549 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1550 ierr = TSAdaptSetStepLimits(adapt, 1551 1.e-12 * units->second, 1552 1.e2 * units->second); CHKERRQ(ierr); 1553 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1554 if (!contsteps) { // print initial condition 1555 if (testChoice == TEST_NONE) { 1556 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1557 } 1558 } else { // continue from time of last output 1559 PetscReal time; 1560 PetscInt count; 1561 PetscViewer viewer; 1562 char filepath[PETSC_MAX_PATH_LEN]; 1563 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1564 user->outputfolder); CHKERRQ(ierr); 1565 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1566 CHKERRQ(ierr); 1567 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1568 CHKERRQ(ierr); 1569 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1570 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1571 } 1572 if (testChoice == TEST_NONE) { 1573 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1574 } 1575 1576 // Solve 1577 start = MPI_Wtime(); 1578 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1579 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1580 cpu_time_used = MPI_Wtime() - start; 1581 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1582 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1583 comm); CHKERRQ(ierr); 1584 if (testChoice == TEST_NONE) { 1585 ierr = PetscPrintf(PETSC_COMM_WORLD, 1586 "Time taken for solution (sec): %g\n", 1587 (double)cpu_time_used); CHKERRQ(ierr); 1588 } 1589 1590 // Get error 1591 if (problem->non_zero_time && testChoice == TEST_NONE) { 1592 Vec Qexact, Qexactloc; 1593 PetscReal norm; 1594 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1595 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1596 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1597 1598 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1599 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1600 1601 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1602 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1603 CeedVectorDestroy(&q0ceed); 1604 ierr = PetscPrintf(PETSC_COMM_WORLD, 1605 "Max Error: %g\n", 1606 (double)norm); CHKERRQ(ierr); 1607 // Clean up vectors 1608 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1609 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1610 } 1611 1612 // Output Statistics 1613 ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr); 1614 if (testChoice == TEST_NONE) { 1615 ierr = PetscPrintf(PETSC_COMM_WORLD, 1616 "Time integrator took %D time steps to reach final time %g\n", 1617 steps, (double)ftime); CHKERRQ(ierr); 1618 } 1619 // Output numerical values from command line 1620 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1621 1622 // compare reference solution values with current run 1623 if (testChoice != TEST_NONE) { 1624 PetscViewer viewer; 1625 // Read reference file 1626 Vec Qref; 1627 PetscReal error, Qrefnorm; 1628 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1629 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1630 CHKERRQ(ierr); 1631 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1632 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1633 1634 // Compute error with respect to reference solution 1635 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1636 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1637 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1638 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1639 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1640 // Check error 1641 if (error > test->testtol) { 1642 ierr = PetscPrintf(PETSC_COMM_WORLD, 1643 "Test failed with error norm %g\n", 1644 (double)error); CHKERRQ(ierr); 1645 } 1646 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1647 } 1648 1649 // Clean up libCEED 1650 CeedVectorDestroy(&qdata); 1651 CeedVectorDestroy(&user->qceed); 1652 CeedVectorDestroy(&user->qdotceed); 1653 CeedVectorDestroy(&user->gceed); 1654 CeedVectorDestroy(&xcorners); 1655 CeedBasisDestroy(&basisq); 1656 CeedBasisDestroy(&basisx); 1657 CeedBasisDestroy(&basisxc); 1658 CeedElemRestrictionDestroy(&restrictq); 1659 CeedElemRestrictionDestroy(&restrictx); 1660 CeedElemRestrictionDestroy(&restrictqdi); 1661 CeedQFunctionDestroy(&qf_setupVol); 1662 CeedQFunctionDestroy(&qf_ics); 1663 CeedQFunctionDestroy(&qf_rhsVol); 1664 CeedQFunctionDestroy(&qf_ifunctionVol); 1665 CeedOperatorDestroy(&op_setupVol); 1666 CeedOperatorDestroy(&op_ics); 1667 CeedOperatorDestroy(&user->op_rhs_vol); 1668 CeedOperatorDestroy(&user->op_ifunction_vol); 1669 CeedDestroy(&ceed); 1670 CeedBasisDestroy(&basisqSur); 1671 CeedBasisDestroy(&basisxSur); 1672 CeedBasisDestroy(&basisxcSur); 1673 CeedQFunctionDestroy(&qf_setupSur); 1674 CeedQFunctionDestroy(&qf_rhsSur); 1675 CeedQFunctionDestroy(&qf_ifunctionSur); 1676 for(CeedInt i=0; i<numOutFlow; i++){ 1677 CeedVectorDestroy(&qdataSur[i]); 1678 CeedElemRestrictionDestroy(&restrictqSur[i]); 1679 CeedElemRestrictionDestroy(&restrictxSur[i]); 1680 CeedElemRestrictionDestroy(&restrictqdiSur[i]); 1681 CeedOperatorDestroy(&op_setupSur[i]); 1682 CeedOperatorDestroy(&user->op_rhs_sur[i]); 1683 CeedOperatorDestroy(&user->op_ifunction_sur[i]); 1684 } 1685 CeedOperatorDestroy(&user->op_rhs); 1686 CeedOperatorDestroy(&user->op_ifunction); 1687 1688 // Clean up PETSc 1689 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1690 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1691 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1692 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1693 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1694 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1695 ierr = PetscFree(units); CHKERRQ(ierr); 1696 ierr = PetscFree(user); CHKERRQ(ierr); 1697 return PetscFinalize(); 1698 } 1699 1700