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