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_volume 1 32 // ./navierstokes -ceed /gpu/occa -problem advection -degree_volume 1 33 // 34 //TESTARGS -ceed {ceed_resource} -test -degree_volume 1 35 36 /// @file 37 /// Navier-Stokes example using PETSc 38 39 const char help[] = "Solve Navier-Stokes using PETSc and libCEED\n"; 40 41 #include <petscts.h> 42 #include <petscdmplex.h> 43 #include <ceed.h> 44 #include <stdbool.h> 45 #include <petscsys.h> 46 #include "common.h" 47 #include "setup-boundary.h" 48 #include "advection.h" 49 #include "advection2d.h" 50 #include "densitycurrent.h" 51 52 // Problem Options 53 typedef enum { 54 NS_DENSITY_CURRENT = 0, 55 NS_ADVECTION = 1, 56 NS_ADVECTION2D = 2, 57 } problemType; 58 static const char *const problemTypes[] = { 59 "density_current", 60 "advection", 61 "advection2d", 62 "problemType","NS_",0 63 }; 64 65 typedef enum { 66 STAB_NONE = 0, 67 STAB_SU = 1, // Streamline Upwind 68 STAB_SUPG = 2, // Streamline Upwind Petrov-Galerkin 69 } StabilizationType; 70 static const char *const StabilizationTypes[] = { 71 "NONE", 72 "SU", 73 "SUPG", 74 "StabilizationType", "STAB_", NULL 75 }; 76 77 // Problem specific data 78 typedef struct { 79 CeedInt dim, qdatasizeVol, qdatasizeSur; 80 CeedQFunctionUser setupVol, setupSur, ics, applyVol_rhs, applySur_rhs, 81 applyVol_ifunction, applySur_ifunction; 82 PetscErrorCode (*bc)(PetscInt, PetscReal, const PetscReal[], PetscInt, 83 PetscScalar[], void *); 84 const char *setupVol_loc, *setupSur_loc, *ics_loc, *applyVol_rhs_loc, 85 *applySur_rhs_loc, *applyVol_ifunction_loc, *applySur_ifunction_loc; 86 const bool non_zero_time; 87 } problemData; 88 89 problemData problemOptions[] = { 90 [NS_DENSITY_CURRENT] = { 91 .dim = 3, 92 .qdatasizeVol = 10, 93 .qdatasizeSur = 4, 94 .setupVol = Setup, 95 .setupVol_loc = Setup_loc, 96 .setupSur = SetupBoundary, 97 .setupSur_loc = SetupBoundary_loc, 98 .ics = ICsDC, 99 .ics_loc = ICsDC_loc, 100 .applyVol_rhs = DC, 101 .applyVol_rhs_loc = DC_loc, 102 .applySur_rhs = DC_Sur, 103 .applySur_rhs_loc = DC_Sur_loc, 104 .applyVol_ifunction = IFunction_DC, 105 .applyVol_ifunction_loc = IFunction_DC_loc, 106 //.applySur_ifunction = IFunction_DC_Sur, 107 //.applySur_ifunction_loc = IFunction_DC_Sur_loc, 108 .bc = Exact_DC, 109 .non_zero_time = false, 110 }, 111 [NS_ADVECTION] = { 112 .dim = 3, 113 .qdatasizeVol = 10, 114 .qdatasizeSur = 4, 115 .setupVol = Setup, 116 .setupVol_loc = Setup_loc, 117 .setupSur = SetupBoundary, 118 .setupSur_loc = SetupBoundary_loc, 119 .ics = ICsAdvection, 120 .ics_loc = ICsAdvection_loc, 121 .applyVol_rhs = Advection, 122 .applyVol_rhs_loc = Advection_loc, 123 .applySur_rhs = Advection_Sur, 124 .applySur_rhs_loc = Advection_Sur_loc, 125 .applyVol_ifunction = IFunction_Advection, 126 .applyVol_ifunction_loc = IFunction_Advection_loc, 127 //.applySur_ifunction = IFunction_Advection_Sur, 128 //.applySur_ifunction_loc = IFunction_Advection_Sur_loc, 129 .bc = Exact_Advection, 130 .non_zero_time = true, 131 }, 132 [NS_ADVECTION2D] = { 133 .dim = 2, 134 .qdatasizeVol = 5, 135 .qdatasizeSur = 3, 136 .setupVol = Setup2d, 137 .setupVol_loc = Setup2d_loc, 138 .setupSur = SetupBoundary2d, 139 .setupSur_loc = SetupBoundary2d_loc, 140 .ics = ICsAdvection2d, 141 .ics_loc = ICsAdvection2d_loc, 142 .applyVol_rhs = Advection2d, 143 .applyVol_rhs_loc = Advection2d_loc, 144 .applySur_rhs = Advection2d_Sur, 145 .applySur_rhs_loc = Advection2d_Sur_loc, 146 .applyVol_ifunction = IFunction_Advection2d, 147 .applyVol_ifunction_loc = IFunction_Advection2d_loc, 148 //.applySur_ifunction = IFunction_Advection2d_Sur, 149 //.applySur_ifunction_loc = IFunction_Advection2d_Sur_loc, 150 .bc = Exact_Advection2d, 151 .non_zero_time = true, 152 }, 153 }; 154 155 // PETSc user data 156 typedef struct User_ *User; 157 typedef struct Units_ *Units; 158 159 struct User_ { 160 MPI_Comm comm; 161 PetscInt outputfreq; 162 DM dm; 163 DM dmviz; 164 Mat interpviz; 165 Ceed ceed; 166 Units units; 167 CeedVector qceed, qdotceed, gceed; 168 CeedOperator op_rhs_vol, op_rhs_sur[6], op_rhs, 169 op_ifunction_vol, op_ifunction_sur[6], op_ifunction; 170 Vec M; 171 char outputfolder[PETSC_MAX_PATH_LEN]; 172 PetscInt contsteps; 173 }; 174 175 struct Units_ { 176 // fundamental units 177 PetscScalar meter; 178 PetscScalar kilogram; 179 PetscScalar second; 180 PetscScalar Kelvin; 181 // derived units 182 PetscScalar Pascal; 183 PetscScalar JperkgK; 184 PetscScalar mpersquareds; 185 PetscScalar WpermK; 186 PetscScalar kgpercubicm; 187 PetscScalar kgpersquaredms; 188 PetscScalar Joulepercubicm; 189 }; 190 191 typedef struct SimpleBC_ *SimpleBC; 192 struct SimpleBC_ { 193 PetscInt nwall, nslip[3], noutflow; 194 PetscInt walls[10], slips[3][10], outflow[6]; 195 }; 196 197 // Essential BC dofs are encoded in closure indices as -(i+1). 198 static PetscInt Involute(PetscInt i) { 199 return i >= 0 ? i : -(i+1); 200 } 201 202 // Utility function to create local CEED restriction 203 static PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 204 CeedInt height, DMLabel domainLabel, PetscInt value, 205 CeedElemRestriction *Erestrict) { 206 207 PetscSection section; 208 PetscInt p, Nelem, Ndof, *erestrict, eoffset, nfields, dim, 209 depth; 210 DMLabel depthLabel; 211 IS depthIS, iterIS; 212 const PetscInt *iterIndices; 213 PetscErrorCode ierr; 214 Vec Uloc; 215 216 PetscFunctionBeginUser; 217 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 218 dim -= height; 219 ierr = DMGetLocalSection(dm,§ion); CHKERRQ(ierr); 220 ierr = PetscSectionGetNumFields(section, &nfields); CHKERRQ(ierr); 221 PetscInt ncomp[nfields], fieldoff[nfields+1]; 222 fieldoff[0] = 0; 223 for (PetscInt f=0; f<nfields; f++) { 224 ierr = PetscSectionGetFieldComponents(section, f, &ncomp[f]); CHKERRQ(ierr); 225 fieldoff[f+1] = fieldoff[f] + ncomp[f]; 226 } 227 228 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 229 ierr = DMPlexGetDepthLabel(dm, &depthLabel); CHKERRQ(ierr); 230 ierr = DMLabelGetStratumIS(depthLabel, depth - height, &depthIS); CHKERRQ(ierr); 231 if (domainLabel) { 232 IS domainIS; 233 ierr = DMLabelGetStratumIS(domainLabel, value, &domainIS); CHKERRQ(ierr); 234 ierr = ISIntersect(depthIS, domainIS, &iterIS); CHKERRQ(ierr); 235 ierr = ISDestroy(&domainIS); CHKERRQ(ierr); 236 ierr = ISDestroy(&depthIS); CHKERRQ(ierr); 237 } else { 238 iterIS = depthIS; 239 } 240 ierr = ISGetLocalSize(iterIS, &Nelem); CHKERRQ(ierr); 241 ierr = ISGetIndices(iterIS, &iterIndices); CHKERRQ(ierr); 242 ierr = PetscMalloc1(Nelem*PetscPowInt(P, dim), &erestrict); CHKERRQ(ierr); 243 for (p=0,eoffset=0; p<Nelem; p++) { 244 PetscInt c = iterIndices[p]; 245 PetscInt numindices, *indices, nnodes; 246 ierr = DMPlexGetClosureIndices(dm, section, section, c, &numindices, 247 &indices, NULL); CHKERRQ(ierr); 248 if (numindices % fieldoff[nfields]) 249 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 250 "Number of closure indices not compatible with Cell %D", c); 251 nnodes = numindices / fieldoff[nfields]; 252 for (PetscInt i=0; i<nnodes; i++) { 253 // Check that indices are blocked by node and thus can be coalesced as a single field with 254 // fieldoff[nfields] = sum(ncomp) components. 255 for (PetscInt f=0; f<nfields; f++) { 256 for (PetscInt j=0; j<ncomp[f]; j++) { 257 if (Involute(indices[fieldoff[f]*nnodes + i*ncomp[f] + j]) 258 != Involute(indices[i*ncomp[0]]) + fieldoff[f] + j) 259 SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, 260 "Cell %D closure indices not interlaced for node %D field %D component %D", 261 c, i, f, j); 262 } 263 } 264 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 265 PetscInt loc = Involute(indices[i*ncomp[0]]); 266 erestrict[eoffset++] = loc; 267 } 268 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, &numindices, 269 &indices, NULL); CHKERRQ(ierr); 270 } 271 if (eoffset != Nelem*PetscPowInt(P, dim)) 272 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 273 "ElemRestriction of size (%D,%D) initialized %D nodes", Nelem, 274 PetscPowInt(P, dim),eoffset); 275 ierr = ISRestoreIndices(iterIS, &iterIndices); CHKERRQ(ierr); 276 ierr = ISDestroy(&iterIS); CHKERRQ(ierr); 277 278 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 279 ierr = VecGetLocalSize(Uloc, &Ndof); CHKERRQ(ierr); 280 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 281 CeedElemRestrictionCreate(ceed, Nelem, PetscPowInt(P, dim), fieldoff[nfields], 282 1, Ndof, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 283 Erestrict); 284 ierr = PetscFree(erestrict); CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 // Utility function to get Ceed Restriction for each domain 289 static PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt ncompx, CeedInt dim, 290 CeedInt height, DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, 291 CeedInt qdatasize, CeedElemRestriction *restrictq, 292 CeedElemRestriction *restrictx, CeedElemRestriction *restrictqdi) { 293 294 DM dmcoord; 295 CeedInt localNelem; 296 CeedInt Qdim = CeedIntPow(Q, dim); 297 PetscErrorCode ierr; 298 299 PetscFunctionBeginUser; 300 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 301 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 302 CHKERRQ(ierr); 303 ierr = CreateRestrictionFromPlex(ceed, dm, P, height, domainLabel, value, restrictq); 304 CHKERRQ(ierr); 305 ierr = CreateRestrictionFromPlex(ceed, dmcoord, 2, height, domainLabel, value, restrictx); 306 CHKERRQ(ierr); 307 CeedElemRestrictionGetNumElements(*restrictq, &localNelem); 308 CeedElemRestrictionCreateStrided(ceed, localNelem, Qdim, 309 qdatasize, qdatasize*localNelem*Qdim, 310 CEED_STRIDES_BACKEND, restrictqdi); 311 PetscFunctionReturn(0); 312 } 313 314 static int CreateVectorFromPetscVec(Ceed ceed, Vec p, CeedVector *v) { 315 PetscErrorCode ierr; 316 PetscInt m; 317 318 PetscFunctionBeginUser; 319 ierr = VecGetLocalSize(p, &m); CHKERRQ(ierr); 320 ierr = CeedVectorCreate(ceed, m, v); CHKERRQ(ierr); 321 PetscFunctionReturn(0); 322 } 323 324 static int VectorPlacePetscVec(CeedVector c, Vec p) { 325 PetscErrorCode ierr; 326 PetscInt mceed,mpetsc; 327 PetscScalar *a; 328 329 PetscFunctionBeginUser; 330 ierr = CeedVectorGetLength(c, &mceed); CHKERRQ(ierr); 331 ierr = VecGetLocalSize(p, &mpetsc); CHKERRQ(ierr); 332 if (mceed != mpetsc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, 333 "Cannot place PETSc Vec of length %D in CeedVector of length %D",mpetsc,mceed); 334 ierr = VecGetArray(p, &a); CHKERRQ(ierr); 335 CeedVectorSetArray(c, CEED_MEM_HOST, CEED_USE_POINTER, a); 336 PetscFunctionReturn(0); 337 } 338 339 static PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 340 PetscBool insertEssential, Vec Qloc, PetscReal time, Vec faceGeomFVM, 341 Vec cellGeomFVM, Vec gradFVM) { 342 PetscErrorCode ierr; 343 Vec Qbc; 344 345 PetscFunctionBegin; 346 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 347 ierr = VecAXPY(Qloc, 1., Qbc); CHKERRQ(ierr); 348 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 349 PetscFunctionReturn(0); 350 } 351 352 // This is the RHS of the ODE, given as u_t = G(t,u) 353 // This function takes in a state vector Q and writes into G 354 static PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *userData) { 355 PetscErrorCode ierr; 356 User user = *(User *)userData; 357 PetscScalar *q, *g; 358 Vec Qloc, Gloc; 359 360 // Global-to-local 361 PetscFunctionBeginUser; 362 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 363 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 364 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 365 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 366 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 367 NULL, NULL, NULL); CHKERRQ(ierr); 368 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 369 370 // Ceed Vectors 371 ierr = VecGetArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 372 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 373 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); 374 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 375 376 // Apply CEED operator 377 CeedOperatorApply(user->op_rhs, user->qceed, user->gceed, 378 CEED_REQUEST_IMMEDIATE); 379 380 // Restore vectors 381 ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 382 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 383 384 ierr = VecZeroEntries(G); CHKERRQ(ierr); 385 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 386 387 // Inverse of the lumped mass matrix 388 ierr = VecPointwiseMult(G, G, user->M); // M is Minv 389 CHKERRQ(ierr); 390 391 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 392 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395 396 static PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Qdot, Vec G, 397 void *userData) { 398 PetscErrorCode ierr; 399 User user = *(User *)userData; 400 const PetscScalar *q, *qdot; 401 PetscScalar *g; 402 Vec Qloc, Qdotloc, Gloc; 403 404 // Global-to-local 405 PetscFunctionBeginUser; 406 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 407 ierr = DMGetLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 408 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 409 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 410 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 411 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 412 NULL, NULL, NULL); CHKERRQ(ierr); 413 ierr = VecZeroEntries(Qdotloc); CHKERRQ(ierr); 414 ierr = DMGlobalToLocal(user->dm, Qdot, INSERT_VALUES, Qdotloc); CHKERRQ(ierr); 415 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 416 417 // Ceed Vectors 418 ierr = VecGetArrayRead(Qloc, &q); CHKERRQ(ierr); 419 ierr = VecGetArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 420 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 421 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, 422 (PetscScalar *)q); 423 CeedVectorSetArray(user->qdotceed, CEED_MEM_HOST, CEED_USE_POINTER, 424 (PetscScalar *)qdot); 425 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 426 427 // Apply CEED operator 428 CeedOperatorApply(user->op_ifunction, user->qceed, user->gceed, 429 CEED_REQUEST_IMMEDIATE); 430 431 // Restore vectors 432 ierr = VecRestoreArrayRead(Qloc, &q); CHKERRQ(ierr); 433 ierr = VecRestoreArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 434 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 435 436 ierr = VecZeroEntries(G); CHKERRQ(ierr); 437 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 438 439 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 440 ierr = DMRestoreLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 441 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 442 PetscFunctionReturn(0); 443 } 444 445 // User provided TS Monitor 446 static PetscErrorCode TSMonitor_NS(TS ts, PetscInt stepno, PetscReal time, 447 Vec Q, void *ctx) { 448 User user = ctx; 449 Vec Qloc; 450 char filepath[PETSC_MAX_PATH_LEN]; 451 PetscViewer viewer; 452 PetscErrorCode ierr; 453 454 // Set up output 455 PetscFunctionBeginUser; 456 // Print every 'outputfreq' steps 457 if (stepno % user->outputfreq != 0) 458 PetscFunctionReturn(0); 459 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 460 ierr = PetscObjectSetName((PetscObject)Qloc, "StateVec"); CHKERRQ(ierr); 461 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 462 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 463 464 // Output 465 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-%03D.vtu", 466 user->outputfolder, stepno + user->contsteps); 467 CHKERRQ(ierr); 468 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), filepath, 469 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 470 ierr = VecView(Qloc, viewer); CHKERRQ(ierr); 471 if (user->dmviz) { 472 Vec Qrefined, Qrefined_loc; 473 char filepath_refined[PETSC_MAX_PATH_LEN]; 474 PetscViewer viewer_refined; 475 476 ierr = DMGetGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 477 ierr = DMGetLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 478 ierr = PetscObjectSetName((PetscObject)Qrefined_loc, "Refined"); 479 CHKERRQ(ierr); 480 ierr = MatInterpolate(user->interpviz, Q, Qrefined); CHKERRQ(ierr); 481 ierr = VecZeroEntries(Qrefined_loc); CHKERRQ(ierr); 482 ierr = DMGlobalToLocal(user->dmviz, Qrefined, INSERT_VALUES, Qrefined_loc); 483 CHKERRQ(ierr); 484 ierr = PetscSNPrintf(filepath_refined, sizeof filepath_refined, 485 "%s/nsrefined-%03D.vtu", 486 user->outputfolder, stepno + user->contsteps); 487 CHKERRQ(ierr); 488 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Qrefined), 489 filepath_refined, 490 FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 491 ierr = VecView(Qrefined_loc, viewer_refined); CHKERRQ(ierr); 492 ierr = DMRestoreLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 493 ierr = DMRestoreGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 494 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 495 } 496 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 497 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 498 499 // Save data in a binary file for continuation of simulations 500 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 501 user->outputfolder); CHKERRQ(ierr); 502 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 503 CHKERRQ(ierr); 504 ierr = VecView(Q, viewer); CHKERRQ(ierr); 505 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 506 507 // Save time stamp 508 // Dimensionalize time back 509 time /= user->units->second; 510 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 511 user->outputfolder); CHKERRQ(ierr); 512 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 513 CHKERRQ(ierr); 514 #if PETSC_VERSION_GE(3,13,0) 515 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 516 #else 517 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 518 #endif 519 CHKERRQ(ierr); 520 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 521 522 PetscFunctionReturn(0); 523 } 524 525 static PetscErrorCode ICs_PetscMultiplicity(CeedOperator op_ics, 526 CeedVector xcorners, CeedVector q0ceed, DM dm, Vec Qloc, Vec Q, 527 CeedElemRestriction restrictq, SetupContext ctxSetup, CeedScalar time) { 528 PetscErrorCode ierr; 529 CeedVector multlvec; 530 Vec Multiplicity, MultiplicityLoc; 531 532 ctxSetup->time = time; 533 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 534 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 535 CeedOperatorApply(op_ics, xcorners, q0ceed, CEED_REQUEST_IMMEDIATE); 536 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 537 ierr = DMLocalToGlobal(dm, Qloc, ADD_VALUES, Q); CHKERRQ(ierr); 538 539 // Fix multiplicity for output of ICs 540 ierr = DMGetLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 541 CeedElemRestrictionCreateVector(restrictq, &multlvec, NULL); 542 ierr = VectorPlacePetscVec(multlvec, MultiplicityLoc); CHKERRQ(ierr); 543 CeedElemRestrictionGetMultiplicity(restrictq, multlvec); 544 CeedVectorDestroy(&multlvec); 545 ierr = DMGetGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 546 ierr = VecZeroEntries(Multiplicity); CHKERRQ(ierr); 547 ierr = DMLocalToGlobal(dm, MultiplicityLoc, ADD_VALUES, Multiplicity); 548 CHKERRQ(ierr); 549 ierr = VecPointwiseDivide(Q, Q, Multiplicity); CHKERRQ(ierr); 550 ierr = VecPointwiseDivide(Qloc, Qloc, MultiplicityLoc); CHKERRQ(ierr); 551 ierr = DMRestoreLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 552 ierr = DMRestoreGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 553 554 PetscFunctionReturn(0); 555 } 556 557 static PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, 558 CeedElemRestriction restrictq, CeedBasis basisq, 559 CeedElemRestriction restrictqdi, CeedVector qdata, Vec M) { 560 PetscErrorCode ierr; 561 CeedQFunction qf_mass; 562 CeedOperator op_mass; 563 CeedVector mceed; 564 Vec Mloc; 565 CeedInt ncompq, qdatasize; 566 567 PetscFunctionBeginUser; 568 CeedElemRestrictionGetNumComponents(restrictq, &ncompq); 569 CeedElemRestrictionGetNumComponents(restrictqdi, &qdatasize); 570 // Create the Q-function that defines the action of the mass operator 571 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 572 CeedQFunctionAddInput(qf_mass, "q", ncompq, CEED_EVAL_INTERP); 573 CeedQFunctionAddInput(qf_mass, "qdata", qdatasize, CEED_EVAL_NONE); 574 CeedQFunctionAddOutput(qf_mass, "v", ncompq, CEED_EVAL_INTERP); 575 576 // Create the mass operator 577 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 578 CeedOperatorSetField(op_mass, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 579 CeedOperatorSetField(op_mass, "qdata", restrictqdi, 580 CEED_BASIS_COLLOCATED, qdata); 581 CeedOperatorSetField(op_mass, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 582 583 ierr = DMGetLocalVector(dm, &Mloc); CHKERRQ(ierr); 584 ierr = VecZeroEntries(Mloc); CHKERRQ(ierr); 585 CeedElemRestrictionCreateVector(restrictq, &mceed, NULL); 586 ierr = VectorPlacePetscVec(mceed, Mloc); CHKERRQ(ierr); 587 588 { 589 // Compute a lumped mass matrix 590 CeedVector onesvec; 591 CeedElemRestrictionCreateVector(restrictq, &onesvec, NULL); 592 CeedVectorSetValue(onesvec, 1.0); 593 CeedOperatorApply(op_mass, onesvec, mceed, CEED_REQUEST_IMMEDIATE); 594 CeedVectorDestroy(&onesvec); 595 CeedOperatorDestroy(&op_mass); 596 CeedVectorDestroy(&mceed); 597 } 598 CeedQFunctionDestroy(&qf_mass); 599 600 ierr = VecZeroEntries(M); CHKERRQ(ierr); 601 ierr = DMLocalToGlobal(dm, Mloc, ADD_VALUES, M); CHKERRQ(ierr); 602 ierr = DMRestoreLocalVector(dm, &Mloc); CHKERRQ(ierr); 603 604 // Invert diagonally lumped mass vector for RHS function 605 ierr = VecReciprocal(M); CHKERRQ(ierr); 606 PetscFunctionReturn(0); 607 } 608 609 PetscErrorCode SetUpDM(DM dm, problemData *problem, PetscInt degree, 610 SimpleBC bc, void *ctxSetup) { 611 PetscErrorCode ierr; 612 613 PetscFunctionBeginUser; 614 { 615 // Configure the finite element space and boundary conditions 616 PetscFE fe; 617 PetscSpace fespace; 618 PetscInt ncompq = 5; 619 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, ncompq, 620 PETSC_FALSE, degree, PETSC_DECIDE, 621 &fe); CHKERRQ(ierr); 622 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 623 ierr = DMAddField(dm,NULL,(PetscObject)fe); CHKERRQ(ierr); 624 ierr = DMCreateDS(dm); CHKERRQ(ierr); 625 /* Wall boundary conditions are zero velocity and zero flux for density and energy */ 626 { 627 PetscInt comps[3] = {1, 2, 3}; 628 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 629 3, comps, (void(*)(void))problem->bc, 630 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 631 } 632 { 633 PetscInt comps[1] = {1}; 634 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 635 1, comps, (void(*)(void))NULL, bc->nslip[0], 636 bc->slips[0], ctxSetup); CHKERRQ(ierr); 637 comps[0] = 2; 638 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 639 1, comps, (void(*)(void))NULL, bc->nslip[1], 640 bc->slips[1], ctxSetup); CHKERRQ(ierr); 641 comps[0] = 3; 642 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 643 1, comps, (void(*)(void))NULL, bc->nslip[2], 644 bc->slips[2], ctxSetup); CHKERRQ(ierr); 645 } 646 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE,NULL); 647 CHKERRQ(ierr); 648 ierr = PetscFEGetBasisSpace(fe, &fespace); CHKERRQ(ierr); 649 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 650 } 651 { 652 // Empty name for conserved field (because there is only one field) 653 PetscSection section; 654 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 655 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 656 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 657 CHKERRQ(ierr); 658 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 659 CHKERRQ(ierr); 660 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 661 CHKERRQ(ierr); 662 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 663 CHKERRQ(ierr); 664 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 665 CHKERRQ(ierr); 666 } 667 PetscFunctionReturn(0); 668 } 669 670 int main(int argc, char **argv) { 671 PetscInt ierr; 672 MPI_Comm comm; 673 DM dm, dmcoord, dmviz, dmvizfine; 674 Mat interpviz; 675 TS ts; 676 TSAdapt adapt; 677 User user; 678 Units units; 679 char ceedresource[4096] = "/cpu/self"; 680 PetscInt localNelemVol, lnodes, steps; 681 const PetscInt ncompq = 5; 682 PetscMPIInt rank; 683 PetscScalar ftime; 684 Vec Q, Qloc, Xloc; 685 Ceed ceed; 686 CeedInt numP_Vol, numQ_Vol; 687 CeedVector xcorners, qdata, q0ceed; 688 CeedBasis basisxVol, basisxcVol, basisqVol; 689 CeedElemRestriction restrictxVol, restrictqVol, restrictqdiVol; 690 CeedQFunction qf_setupVol, qf_ics, qf_rhsVol, qf_ifunctionVol; 691 CeedOperator op_setupVol, op_ics; 692 CeedScalar Rd; 693 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 694 kgpersquaredms, Joulepercubicm; 695 problemType problemChoice; 696 problemData *problem = NULL; 697 StabilizationType stab; 698 PetscBool test, implicit; 699 PetscInt viz_refine = 0; 700 struct SimpleBC_ bc = { 701 .nwall = 6, 702 .walls = {1,2,3,4,5,6}, 703 }; 704 double start, cpu_time_used; 705 706 // Create the libCEED contexts 707 PetscScalar meter = 1e-2; // 1 meter in scaled length units 708 PetscScalar second = 1e-2; // 1 second in scaled time units 709 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 710 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 711 CeedScalar theta0 = 300.; // K 712 CeedScalar thetaC = -15.; // K 713 CeedScalar P0 = 1.e5; // Pa 714 CeedScalar N = 0.01; // 1/s 715 CeedScalar cv = 717.; // J/(kg K) 716 CeedScalar cp = 1004.; // J/(kg K) 717 CeedScalar g = 9.81; // m/s^2 718 CeedScalar lambda = -2./3.; // - 719 CeedScalar mu = 75.; // Pa s, dynamic viscosity 720 // mu = 75 is not physical for air, but is good for numerical stability 721 CeedScalar k = 0.02638; // W/(m K) 722 CeedScalar CtauS = 0.; // dimensionless 723 CeedScalar strong_form = 0.; // [0,1] 724 PetscScalar lx = 8000.; // m 725 PetscScalar ly = 8000.; // m 726 PetscScalar lz = 4000.; // m 727 CeedScalar rc = 1000.; // m (Radius of bubble) 728 PetscScalar resx = 1000.; // m (resolution in x) 729 PetscScalar resy = 1000.; // m (resolution in y) 730 PetscScalar resz = 1000.; // m (resolution in z) 731 PetscInt outputfreq = 10; // - 732 PetscInt contsteps = 0; // - 733 PetscInt degreeVol = 1; // - 734 PetscInt degreeSur = 1; // - 735 PetscInt qextraVol = 2; // - 736 PetscInt qextraSur = 2; // - 737 DMBoundaryType periodicity[] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 738 DM_BOUNDARY_NONE 739 }; 740 PetscReal center[3], dc_axis[3] = {0, 0, 0}; 741 742 ierr = PetscInitialize(&argc, &argv, NULL, help); 743 if (ierr) return ierr; 744 745 // Allocate PETSc context 746 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 747 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 748 749 // Parse command line options 750 comm = PETSC_COMM_WORLD; 751 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 752 NULL); CHKERRQ(ierr); 753 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 754 NULL, ceedresource, ceedresource, 755 sizeof(ceedresource), NULL); CHKERRQ(ierr); 756 ierr = PetscOptionsBool("-test", "Run in test mode", 757 NULL, test=PETSC_FALSE, &test, NULL); CHKERRQ(ierr); 758 problemChoice = NS_DENSITY_CURRENT; 759 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 760 problemTypes, (PetscEnum)problemChoice, 761 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 762 problem = &problemOptions[problemChoice]; 763 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 764 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 765 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 766 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 767 NULL, implicit=PETSC_FALSE, &implicit, NULL); 768 CHKERRQ(ierr); 769 { 770 PetscInt len, len1; 771 PetscBool flg, flg1; 772 ierr = PetscOptionsIntArray("-bc_wall", 773 "Use wall boundary conditions on this list of faces", 774 NULL, bc.walls, 775 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 776 &len), &flg); CHKERRQ(ierr); 777 if (flg) bc.nwall = len; 778 for (PetscInt j=0; j<3; j++) { 779 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 780 ierr = PetscOptionsIntArray(flags[j], 781 "Use slip boundary conditions on this list of faces", 782 NULL, bc.slips[j], 783 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 784 &len), &flg); 785 CHKERRQ(ierr); 786 if (flg) bc.nslip[j] = len; 787 } 788 ierr = PetscOptionsIntArray("-bc_outflow", 789 "Use outflow boundary conditions on this list of faces", 790 NULL, bc.outflow, 791 (len1 = sizeof(bc.outflow) / sizeof(bc.outflow[0]), 792 &len1), &flg1); CHKERRQ(ierr); 793 if (flg1) bc.noutflow = len1; 794 } 795 ierr = PetscOptionsInt("-viz_refine", 796 "Regular refinement levels for visualization", 797 NULL, viz_refine, &viz_refine, NULL); 798 CHKERRQ(ierr); 799 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 800 NULL, meter, &meter, NULL); CHKERRQ(ierr); 801 meter = fabs(meter); 802 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 803 NULL, second, &second, NULL); CHKERRQ(ierr); 804 second = fabs(second); 805 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 806 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 807 kilogram = fabs(kilogram); 808 ierr = PetscOptionsScalar("-units_Kelvin", 809 "1 Kelvin in scaled temperature units", 810 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 811 Kelvin = fabs(Kelvin); 812 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 813 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 814 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 815 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 816 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 817 NULL, P0, &P0, NULL); CHKERRQ(ierr); 818 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 819 NULL, N, &N, NULL); CHKERRQ(ierr); 820 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 821 NULL, cv, &cv, NULL); CHKERRQ(ierr); 822 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 823 NULL, cp, &cp, NULL); CHKERRQ(ierr); 824 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 825 NULL, g, &g, NULL); CHKERRQ(ierr); 826 ierr = PetscOptionsScalar("-lambda", 827 "Stokes hypothesis second viscosity coefficient", 828 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 829 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 830 NULL, mu, &mu, NULL); CHKERRQ(ierr); 831 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 832 NULL, k, &k, NULL); CHKERRQ(ierr); 833 ierr = PetscOptionsScalar("-CtauS", 834 "Scale coefficient for tau (nondimensional)", 835 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 836 ierr = PetscOptionsScalar("-strong_form", 837 "Strong (1) or weak/integrated by parts (0) advection residual", 838 NULL, strong_form, &strong_form, NULL); 839 CHKERRQ(ierr); 840 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 841 NULL, lx, &lx, NULL); CHKERRQ(ierr); 842 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 843 NULL, ly, &ly, NULL); CHKERRQ(ierr); 844 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 845 NULL, lz, &lz, NULL); CHKERRQ(ierr); 846 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 847 NULL, rc, &rc, NULL); CHKERRQ(ierr); 848 ierr = PetscOptionsScalar("-resx","Target resolution in x", 849 NULL, resx, &resx, NULL); CHKERRQ(ierr); 850 ierr = PetscOptionsScalar("-resy","Target resolution in y", 851 NULL, resy, &resy, NULL); CHKERRQ(ierr); 852 ierr = PetscOptionsScalar("-resz","Target resolution in z", 853 NULL, resz, &resz, NULL); CHKERRQ(ierr); 854 PetscInt n = problem->dim; 855 ierr = PetscOptionsEnumArray("-periodicity", "Periodicity per direction", 856 NULL, DMBoundaryTypes, (PetscEnum *)periodicity, 857 &n, NULL); CHKERRQ(ierr); 858 n = problem->dim; 859 center[0] = 0.5 * lx; 860 center[1] = 0.5 * ly; 861 center[2] = 0.5 * lz; 862 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 863 NULL, center, &n, NULL); CHKERRQ(ierr); 864 n = problem->dim; 865 ierr = PetscOptionsRealArray("-dc_axis", 866 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 867 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 868 { 869 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 870 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 871 if (norm > 0) { 872 for (int i=0; i<3; i++) dc_axis[i] /= norm; 873 } 874 } 875 ierr = PetscOptionsInt("-output_freq", 876 "Frequency of output, in number of steps", 877 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 878 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 879 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 880 ierr = PetscOptionsInt("-degree_volume", "Polynomial degree of finite elements for the Volume", 881 NULL, degreeVol, °reeVol, NULL); CHKERRQ(ierr); 882 ierr = PetscOptionsInt("-qextra_volume", "Number of extra quadrature points for the Volume", 883 NULL, qextraVol, &qextraVol, NULL); CHKERRQ(ierr); 884 PetscStrncpy(user->outputfolder, ".", 2); 885 ierr = PetscOptionsString("-of", "Output folder", 886 NULL, user->outputfolder, user->outputfolder, 887 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 888 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 889 890 // Define derived units 891 Pascal = kilogram / (meter * PetscSqr(second)); 892 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 893 mpersquareds = meter / PetscSqr(second); 894 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 895 kgpercubicm = kilogram / pow(meter,3); 896 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 897 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 898 899 // Scale variables to desired units 900 theta0 *= Kelvin; 901 thetaC *= Kelvin; 902 P0 *= Pascal; 903 N *= (1./second); 904 cv *= JperkgK; 905 cp *= JperkgK; 906 Rd = cp - cv; 907 g *= mpersquareds; 908 mu *= Pascal * second; 909 k *= WpermK; 910 lx = fabs(lx) * meter; 911 ly = fabs(ly) * meter; 912 lz = fabs(lz) * meter; 913 rc = fabs(rc) * meter; 914 resx = fabs(resx) * meter; 915 resy = fabs(resy) * meter; 916 resz = fabs(resz) * meter; 917 for (int i=0; i<3; i++) center[i] *= meter; 918 919 const CeedInt dim = problem->dim, ncompx = problem->dim, 920 qdatasizeVol = problem->qdatasizeVol; 921 // Set up the libCEED context 922 struct SetupContext_ ctxSetup = { 923 .theta0 = theta0, 924 .thetaC = thetaC, 925 .P0 = P0, 926 .N = N, 927 .cv = cv, 928 .cp = cp, 929 .Rd = Rd, 930 .g = g, 931 .rc = rc, 932 .lx = lx, 933 .ly = ly, 934 .lz = lz, 935 .periodicity0 = periodicity[0], 936 .periodicity1 = periodicity[1], 937 .periodicity2 = periodicity[2], 938 .center[0] = center[0], 939 .center[1] = center[1], 940 .center[2] = center[2], 941 .dc_axis[0] = dc_axis[0], 942 .dc_axis[1] = dc_axis[1], 943 .dc_axis[2] = dc_axis[2], 944 .time = 0, 945 }; 946 947 { 948 const PetscReal scale[3] = {lx, ly, lz}; 949 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 950 periodicity, PETSC_TRUE, &dm); 951 CHKERRQ(ierr); 952 } 953 if (1) { 954 DM dmDist = NULL; 955 PetscPartitioner part; 956 957 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 958 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 959 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 960 if (dmDist) { 961 ierr = DMDestroy(&dm); CHKERRQ(ierr); 962 dm = dmDist; 963 } 964 } 965 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 966 967 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 968 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 969 ierr = SetUpDM(dm, problem, degreeVol, &bc, &ctxSetup); CHKERRQ(ierr); 970 if (!test) { 971 ierr = PetscPrintf(PETSC_COMM_WORLD, 972 "Degree of Volumetric FEM Space: %D\n", 973 degreeVol); CHKERRQ(ierr); 974 } 975 dmviz = NULL; 976 interpviz = NULL; 977 if (viz_refine) { 978 DM dmhierarchy[viz_refine+1]; 979 980 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 981 dmhierarchy[0] = dm; 982 for (PetscInt i = 0, d = degreeVol; i < viz_refine; i++) { 983 Mat interp_next; 984 985 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 986 CHKERRQ(ierr); 987 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 988 d = (d + 1) / 2; 989 if (i + 1 == viz_refine) d = 1; 990 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 991 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 992 &interp_next, NULL); CHKERRQ(ierr); 993 if (!i) interpviz = interp_next; 994 else { 995 Mat C; 996 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 997 PETSC_DECIDE, &C); CHKERRQ(ierr); 998 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 999 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1000 interpviz = C; 1001 } 1002 } 1003 for (PetscInt i=1; i<viz_refine; i++) { 1004 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1005 } 1006 dmviz = dmhierarchy[viz_refine]; 1007 } 1008 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1009 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1010 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1011 lnodes /= ncompq; 1012 1013 { 1014 // Print grid information 1015 CeedInt gdofs, odofs; 1016 int comm_size; 1017 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1018 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1019 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1020 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1021 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1022 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1023 if (!test) { 1024 ierr = PetscPrintf(comm, "Global FEM dofs: %D (%D owned) on %d ranks\n", 1025 gdofs, odofs, comm_size); CHKERRQ(ierr); 1026 ierr = PetscPrintf(comm, "Local FEM nodes: %D\n", lnodes); CHKERRQ(ierr); 1027 ierr = PetscPrintf(comm, "dm_plex_box_faces: %s\n", box_faces_str); 1028 CHKERRQ(ierr); 1029 } 1030 1031 } 1032 1033 // Set up global mass vector 1034 ierr = VecDuplicate(Q,&user->M); CHKERRQ(ierr); 1035 1036 // Set up CEED 1037 // CEED Bases 1038 CeedInit(ceedresource, &ceed); 1039 numP_Vol = degreeVol + 1; 1040 numQ_Vol = numP_Vol + qextraVol; 1041 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP_Vol, numQ_Vol, CEED_GAUSS, 1042 &basisqVol); 1043 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ_Vol, CEED_GAUSS, 1044 &basisxVol); 1045 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP_Vol, 1046 CEED_GAUSS_LOBATTO, &basisxcVol); 1047 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1048 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1049 CHKERRQ(ierr); 1050 1051 // CEED Restrictions 1052 // Restrictions on the Volume 1053 ierr = GetRestrictionForDomain(ceed, dm, ncompx, dim, 0, 0, 0, numP_Vol, numQ_Vol, 1054 qdatasizeVol, &restrictqVol, &restrictxVol, 1055 &restrictqdiVol); CHKERRQ(ierr); 1056 1057 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1058 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1059 1060 // Create the CEED vectors that will be needed in setup 1061 CeedInt NqptsVol; 1062 CeedBasisGetNumQuadraturePoints(basisqVol, &NqptsVol); 1063 CeedElemRestrictionGetNumElements(restrictqVol, &localNelemVol); 1064 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1065 CeedElemRestrictionCreateVector(restrictqVol, &q0ceed, NULL); 1066 1067 // Create the Q-function that builds the quadrature data for the NS operator 1068 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1069 &qf_setupVol); 1070 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1071 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1072 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1073 1074 // Create the Q-function that sets the ICs of the operator 1075 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1076 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1077 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1078 1079 qf_rhsVol = NULL; 1080 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1081 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1082 problem->applyVol_rhs_loc, &qf_rhsVol); 1083 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1084 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1085 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1086 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1087 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1088 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1089 } 1090 1091 qf_ifunctionVol = NULL; 1092 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1093 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1094 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1095 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1096 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1097 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1098 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1099 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1100 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1101 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1102 } 1103 1104 // Create the operator that builds the quadrature data for the NS operator 1105 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1106 CeedOperatorSetField(op_setupVol, "dx", restrictxVol, basisxVol, CEED_VECTOR_ACTIVE); 1107 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1108 basisxVol, CEED_VECTOR_NONE); 1109 CeedOperatorSetField(op_setupVol, "qdata", restrictqdiVol, 1110 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1111 1112 // Create the operator that sets the ICs 1113 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1114 CeedOperatorSetField(op_ics, "x", restrictxVol, basisxcVol, CEED_VECTOR_ACTIVE); 1115 CeedOperatorSetField(op_ics, "q0", restrictqVol, 1116 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1117 1118 CeedElemRestrictionCreateVector(restrictqVol, &user->qceed, NULL); 1119 CeedElemRestrictionCreateVector(restrictqVol, &user->qdotceed, NULL); 1120 CeedElemRestrictionCreateVector(restrictqVol, &user->gceed, NULL); 1121 1122 if (qf_rhsVol) { // Create the RHS physics operator 1123 CeedOperator op; 1124 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1125 CeedOperatorSetField(op, "q", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1126 CeedOperatorSetField(op, "dq", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1127 CeedOperatorSetField(op, "qdata", restrictqdiVol, 1128 CEED_BASIS_COLLOCATED, qdata); 1129 CeedOperatorSetField(op, "x", restrictxVol, basisxVol, xcorners); 1130 CeedOperatorSetField(op, "v", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1131 CeedOperatorSetField(op, "dv", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1132 user->op_rhs = op; 1133 } 1134 1135 if (qf_ifunctionVol) { // Create the IFunction operator 1136 CeedOperator op; 1137 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1138 CeedOperatorSetField(op, "q", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1139 CeedOperatorSetField(op, "dq", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1140 CeedOperatorSetField(op, "qdot", restrictqVol, basisqVol, user->qdotceed); 1141 CeedOperatorSetField(op, "qdata", restrictqdiVol, 1142 CEED_BASIS_COLLOCATED, qdata); 1143 CeedOperatorSetField(op, "x", restrictxVol, basisxVol, xcorners); 1144 CeedOperatorSetField(op, "v", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1145 CeedOperatorSetField(op, "dv", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1146 user->op_ifunction = op; 1147 } 1148 1149 //--------------------------------------------------------------------------------------// 1150 // Outflow Boundary Condition 1151 //--------------------------------------------------------------------------------------// 1152 // Set up CEED for the boundaries 1153 CeedInt numP_Sur, numQ_Sur; 1154 CeedInt height = 1; 1155 CeedInt dimSur = dim - height; 1156 numP_Sur = degreeSur + 1; 1157 numQ_Sur = numP_Sur + qextraSur; 1158 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1159 CeedBasis basisxSur, basisxcSur, basisqSur; 1160 CeedElemRestriction restrictxSur[6], restrictqSur[6], restrictqdiSur[6]; 1161 CeedInt NqptsSur; 1162 PetscInt localNelemSur[6]; 1163 CeedVector qdataSur[6], qdataSur_[6]; 1164 CeedQFunction qf_setupSur, qf_rhsSur, qf_ifunctionSur; 1165 CeedOperator op_setupSur[6], op_setupSur_[6]; 1166 PetscInt numOutFlow = bc.noutflow; 1167 DMLabel domainLabel; 1168 1169 // Get Label for the boundaries 1170 ierr = DMGetLabel(dm, "Face Sets", &domainLabel); CHKERRQ(ierr); 1171 1172 // CEED bases for the boundaries 1173 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, CEED_GAUSS, 1174 &basisqSur); 1175 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1176 &basisxSur); 1177 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1178 CEED_GAUSS_LOBATTO, &basisxcSur); 1179 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1180 1181 // Create the Q-function that builds the quadrature data for the Surface operator 1182 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1183 &qf_setupSur); 1184 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1185 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1186 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1187 1188 qf_rhsSur = NULL; 1189 if (problem->applySur_rhs) { // Create the Q-function that defines the action of the RHS operator on the Surface 1190 CeedQFunctionCreateInterior(ceed, 1, problem->applySur_rhs, 1191 problem->applySur_rhs_loc, &qf_rhsSur); 1192 CeedQFunctionAddInput(qf_rhsSur, "q", ncompq, CEED_EVAL_INTERP); 1193 CeedQFunctionAddInput(qf_rhsSur, "dq", ncompq*dim, CEED_EVAL_GRAD); //assumed volumetric elements 1194 CeedQFunctionAddInput(qf_rhsSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1195 CeedQFunctionAddInput(qf_rhsSur, "x", ncompx, CEED_EVAL_INTERP); 1196 CeedQFunctionAddInput(qf_rhsSur, "qdata", qdatasizeVol, CEED_EVAL_NONE); //assumed volumetric elements 1197 CeedQFunctionAddOutput(qf_rhsSur, "v", ncompq, CEED_EVAL_INTERP); 1198 } 1199 qf_ifunctionSur = NULL; 1200 if (problem->applySur_ifunction) { // Create the Q-function that defines the action of the IFunction 1201 CeedQFunctionCreateInterior(ceed, 1, problem->applySur_ifunction, 1202 problem->applySur_ifunction_loc, &qf_ifunctionSur); 1203 CeedQFunctionAddInput(qf_ifunctionSur, "q", ncompq, CEED_EVAL_INTERP); 1204 CeedQFunctionAddInput(qf_ifunctionSur, "dq", ncompq*dim, CEED_EVAL_GRAD); //assumed volumetric elements 1205 CeedQFunctionAddInput(qf_ifunctionSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1206 CeedQFunctionAddInput(qf_ifunctionSur, "x", ncompx, CEED_EVAL_INTERP); 1207 CeedQFunctionAddInput(qf_ifunctionSur, "qdata", qdatasizeVol, CEED_EVAL_NONE); //assumed volumetric elements 1208 CeedQFunctionAddOutput(qf_ifunctionSur, "v", ncompq, CEED_EVAL_INTERP); 1209 } 1210 // Create CEED Operator for each face 1211 for(CeedInt i=0; i<numOutFlow; i++){ 1212 1213 ierr = GetRestrictionForDomain(ceed, dm, ncompx, dimSur, height, domainLabel, bc.outflow[i], numP_Sur, 1214 numQ_Sur, qdatasizeSur, &restrictqSur[i], &restrictxSur[i], 1215 &restrictqdiSur[i]); CHKERRQ(ierr); 1216 // Create the CEED vectors that will be needed in setup 1217 CeedElemRestrictionGetNumElements(restrictqSur[i], &localNelemSur[i]); 1218 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, &qdataSur[i]); 1219 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, &qdataSur_[i]); 1220 1221 // Create the operator that builds the quadrature data for the NS operator 1222 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur[i]); 1223 CeedOperatorSetField(op_setupSur[i], "dx", restrictxSur[i], basisxSur, CEED_VECTOR_ACTIVE); 1224 CeedOperatorSetField(op_setupSur[i], "weight", CEED_ELEMRESTRICTION_NONE, 1225 basisxSur, CEED_VECTOR_NONE); 1226 CeedOperatorSetField(op_setupSur[i], "qdataSur", restrictqdiSur[i], 1227 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1228 1229 // Utility operator that builds the quadrature data for computing viscous terms 1230 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupSur_[i]); 1231 CeedOperatorSetField(op_setupSur_[i], "dx", restrictxSur[i], basisxVol, CEED_VECTOR_ACTIVE); 1232 CeedOperatorSetField(op_setupSur_[i], "weight", CEED_ELEMRESTRICTION_NONE, 1233 basisxVol, CEED_VECTOR_NONE); 1234 CeedOperatorSetField(op_setupSur_[i], "qdata", restrictqdiSur[i], 1235 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1236 1237 if (qf_rhsSur) { // Create the RHS physics operator 1238 CeedOperator op; 1239 CeedOperatorCreate(ceed, qf_rhsSur, NULL, NULL, &op); 1240 CeedOperatorSetField(op, "q", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1241 CeedOperatorSetField(op, "dq", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1242 CeedOperatorSetField(op, "qdataSur", restrictqdiSur[i], 1243 CEED_BASIS_COLLOCATED, qdataSur[i]); 1244 CeedOperatorSetField(op, "x", restrictxSur[i], basisxSur, xcorners); 1245 CeedOperatorSetField(op, "qdata", restrictqdiSur[i], 1246 CEED_BASIS_COLLOCATED, qdataSur_[i]); 1247 CeedOperatorSetField(op, "v", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1248 user->op_rhs_sur[i] = op; 1249 } 1250 1251 if (qf_ifunctionSur) { // Create the IFunction operator 1252 CeedOperator op; 1253 CeedOperatorCreate(ceed, qf_ifunctionSur, NULL, NULL, &op); 1254 CeedOperatorSetField(op, "q", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1255 CeedOperatorSetField(op, "dq", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1256 CeedOperatorSetField(op, "qdataSur", restrictqdiSur[i], 1257 CEED_BASIS_COLLOCATED, qdataSur[i]); 1258 CeedOperatorSetField(op, "x", restrictxSur[i], basisxSur, xcorners); 1259 CeedOperatorSetField(op, "qdata", restrictqdiSur[i], 1260 CEED_BASIS_COLLOCATED, qdataSur_[i]); 1261 CeedOperatorSetField(op, "v", restrictqSur[i], basisqSur, CEED_VECTOR_ACTIVE); 1262 user->op_ifunction_sur[i] = op; 1263 } 1264 } 1265 // Composite Operaters 1266 if (user->op_ifunction_vol) { 1267 if (numOutFlow>0) { 1268 // Composite Operators for the IFunction 1269 CeedCompositeOperatorCreate(ceed, &user->op_ifunction); 1270 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_vol); 1271 for(CeedInt i=0; i<numOutFlow; i++){ 1272 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_sur[i]); 1273 } 1274 } else { 1275 user->op_ifunction = user->op_ifunction_vol; 1276 } 1277 } 1278 if (user->op_rhs_vol) { 1279 if (numOutFlow == 1) { 1280 // Composite Operators for the RHS 1281 CeedCompositeOperatorCreate(ceed, &user->op_rhs); 1282 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_vol); 1283 for(CeedInt i=0; i<numOutFlow; i++){ 1284 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_sur[i]); 1285 } 1286 } else { 1287 user->op_rhs = user->op_rhs_vol; 1288 } 1289 } 1290 //--------------------------------------------------------------------------------------// 1291 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1292 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1293 struct Advection2dContext_ ctxAdvection2d = { 1294 .CtauS = CtauS, 1295 .strong_form = strong_form, 1296 .stabilization = stab, 1297 }; 1298 switch (problemChoice) { 1299 case NS_DENSITY_CURRENT: 1300 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1301 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1302 sizeof ctxNS); 1303 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxNS, sizeof ctxNS); 1304 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxNS, 1305 sizeof ctxNS); 1306 break; 1307 case NS_ADVECTION: 1308 case NS_ADVECTION2D: 1309 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1310 sizeof ctxAdvection2d); 1311 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1312 sizeof ctxAdvection2d); 1313 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxAdvection2d, 1314 sizeof ctxAdvection2d); 1315 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxAdvection2d, 1316 sizeof ctxAdvection2d); 1317 } 1318 1319 // Set up PETSc context 1320 // Set up units structure 1321 units->meter = meter; 1322 units->kilogram = kilogram; 1323 units->second = second; 1324 units->Kelvin = Kelvin; 1325 units->Pascal = Pascal; 1326 units->JperkgK = JperkgK; 1327 units->mpersquareds = mpersquareds; 1328 units->WpermK = WpermK; 1329 units->kgpercubicm = kgpercubicm; 1330 units->kgpersquaredms = kgpersquaredms; 1331 units->Joulepercubicm = Joulepercubicm; 1332 1333 // Set up user structure 1334 user->comm = comm; 1335 user->outputfreq = outputfreq; 1336 user->contsteps = contsteps; 1337 user->units = units; 1338 user->dm = dm; 1339 user->dmviz = dmviz; 1340 user->interpviz = interpviz; 1341 user->ceed = ceed; 1342 1343 // Calculate qdata and ICs 1344 // Set up state global and local vectors 1345 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1346 1347 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1348 1349 // Apply Setup Ceed Operators 1350 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1351 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1352 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictqVol, basisqVol, restrictqdiVol, qdata, 1353 user->M); CHKERRQ(ierr); 1354 1355 ierr = ICs_PetscMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictqVol, 1356 &ctxSetup, 0.0); 1357 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1358 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1359 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1360 // slower execution. 1361 Vec Qbc; 1362 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1363 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1364 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1365 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1366 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1367 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1368 ierr = PetscObjectComposeFunction((PetscObject)dm, 1369 "DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_NS); CHKERRQ(ierr); 1370 } 1371 1372 MPI_Comm_rank(comm, &rank); 1373 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1374 // Gather initial Q values 1375 // In case of continuation of simulation, set up initial values from binary file 1376 if (contsteps) { // continue from existent solution 1377 PetscViewer viewer; 1378 char filepath[PETSC_MAX_PATH_LEN]; 1379 // Read input 1380 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1381 user->outputfolder); 1382 CHKERRQ(ierr); 1383 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1384 CHKERRQ(ierr); 1385 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1386 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1387 } else { 1388 //ierr = DMLocalToGlobal(dm, Qloc, INSERT_VALUES, Q);CHKERRQ(ierr); 1389 } 1390 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1391 1392 // Create and setup TS 1393 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1394 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1395 if (implicit) { 1396 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1397 if (user->op_ifunction) { 1398 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1399 } else { // Implicit integrators can fall back to using an RHSFunction 1400 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1401 } 1402 } else { 1403 if (!user->op_rhs) SETERRQ(comm,PETSC_ERR_ARG_NULL, 1404 "Problem does not provide RHSFunction"); 1405 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1406 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1407 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1408 } 1409 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1410 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1411 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1412 if (test) {ierr = TSSetMaxSteps(ts, 1); CHKERRQ(ierr);} 1413 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1414 ierr = TSAdaptSetStepLimits(adapt, 1415 1.e-12 * units->second, 1416 1.e2 * units->second); CHKERRQ(ierr); 1417 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1418 if (!contsteps) { // print initial condition 1419 if (!test) { 1420 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1421 } 1422 } else { // continue from time of last output 1423 PetscReal time; 1424 PetscInt count; 1425 PetscViewer viewer; 1426 char filepath[PETSC_MAX_PATH_LEN]; 1427 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1428 user->outputfolder); CHKERRQ(ierr); 1429 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1430 CHKERRQ(ierr); 1431 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1432 CHKERRQ(ierr); 1433 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1434 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1435 } 1436 if (!test) { 1437 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1438 } 1439 1440 // Solve 1441 start = MPI_Wtime(); 1442 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1443 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1444 cpu_time_used = MPI_Wtime() - start; 1445 ierr = TSGetSolveTime(ts,&ftime); CHKERRQ(ierr); 1446 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1447 comm); CHKERRQ(ierr); 1448 if (!test) { 1449 ierr = PetscPrintf(PETSC_COMM_WORLD, 1450 "Time taken for solution: %g\n", 1451 (double)cpu_time_used); CHKERRQ(ierr); 1452 } 1453 1454 // Get error 1455 if (problem->non_zero_time && !test) { 1456 Vec Qexact, Qexactloc; 1457 PetscReal norm; 1458 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1459 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1460 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1461 1462 ierr = ICs_PetscMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1463 restrictqVol, &ctxSetup, ftime); CHKERRQ(ierr); 1464 1465 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1466 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1467 CeedVectorDestroy(&q0ceed); 1468 ierr = PetscPrintf(PETSC_COMM_WORLD, 1469 "Max Error: %g\n", 1470 (double)norm); CHKERRQ(ierr); 1471 } 1472 1473 // Output Statistics 1474 ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr); 1475 if (!test) { 1476 ierr = PetscPrintf(PETSC_COMM_WORLD, 1477 "Time integrator took %D time steps to reach final time %g\n", 1478 steps,(double)ftime); CHKERRQ(ierr); 1479 } 1480 1481 // Clean up libCEED 1482 CeedVectorDestroy(&qdata); 1483 CeedVectorDestroy(&user->qceed); 1484 CeedVectorDestroy(&user->qdotceed); 1485 CeedVectorDestroy(&user->gceed); 1486 CeedVectorDestroy(&xcorners); 1487 CeedBasisDestroy(&basisqVol); 1488 CeedBasisDestroy(&basisxVol); 1489 CeedBasisDestroy(&basisxcVol); 1490 CeedElemRestrictionDestroy(&restrictqVol); 1491 CeedElemRestrictionDestroy(&restrictxVol); 1492 CeedElemRestrictionDestroy(&restrictqdiVol); 1493 CeedQFunctionDestroy(&qf_setupVol); 1494 CeedQFunctionDestroy(&qf_ics); 1495 CeedQFunctionDestroy(&qf_rhsVol); 1496 CeedQFunctionDestroy(&qf_ifunctionVol); 1497 CeedOperatorDestroy(&op_setupVol); 1498 CeedOperatorDestroy(&op_ics); 1499 CeedOperatorDestroy(&user->op_rhs_vol); 1500 CeedOperatorDestroy(&user->op_ifunction_vol); 1501 CeedDestroy(&ceed); 1502 CeedBasisDestroy(&basisqSur); 1503 CeedBasisDestroy(&basisxSur); 1504 CeedBasisDestroy(&basisxcSur); 1505 CeedQFunctionDestroy(&qf_setupSur); 1506 CeedQFunctionDestroy(&qf_rhsSur); 1507 CeedQFunctionDestroy(&qf_ifunctionSur); 1508 for(CeedInt i=0; i<numOutFlow; i++){ 1509 CeedVectorDestroy(&qdataSur[i]); 1510 CeedVectorDestroy(&qdataSur_[i]); 1511 CeedElemRestrictionDestroy(&restrictqSur[i]); 1512 CeedElemRestrictionDestroy(&restrictxSur[i]); 1513 CeedElemRestrictionDestroy(&restrictqdiSur[i]); 1514 CeedOperatorDestroy(&op_setupSur[i]); 1515 CeedOperatorDestroy(&op_setupSur_[i]); 1516 CeedOperatorDestroy(&user->op_rhs_sur[i]); 1517 CeedOperatorDestroy(&user->op_ifunction_sur[i]); 1518 } 1519 CeedOperatorDestroy(&user->op_rhs); 1520 CeedOperatorDestroy(&user->op_ifunction); 1521 1522 // Clean up PETSc 1523 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1524 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1525 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1526 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1527 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1528 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1529 ierr = PetscFree(units); CHKERRQ(ierr); 1530 ierr = PetscFree(user); CHKERRQ(ierr); 1531 return PetscFinalize(); 1532 } 1533 1534