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