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, op_rhs, 169 op_ifunction_vol, op_ifunction_sur, 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]; 194 PetscInt walls[10], slips[3][10]; 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, CeedInt 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, CeedInt 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, localNelemSur, 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, numP_Sur, numQ_Vol, numQ_Sur; 687 CeedVector xcorners, qdata, qdataSur, q0ceedVol, q0ceedSur; 688 CeedBasis basisxVol, basisxcVol, basisqVol, basisxSur, basisxcSur, basisqSur; 689 CeedElemRestriction restrictxVol, restrictqVol, restrictqdiVol, 690 restrictxSur, restrictqSur, restrictqdiSur; 691 CeedQFunction qf_setupVol, qf_setupSur, qf_ics, qf_rhsVol, qf_rhsSur, 692 qf_ifunctionVol, qf_ifunctionSur; 693 CeedOperator op_setupVol, op_setupSur, op_ics; 694 CeedScalar Rd; 695 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 696 kgpersquaredms, Joulepercubicm; 697 problemType problemChoice; 698 problemData *problem = NULL; 699 StabilizationType stab; 700 PetscBool test, implicit; 701 PetscInt viz_refine = 0; 702 struct SimpleBC_ bc = { 703 .nwall = 6, 704 .walls = {1,2,3,4,5,6}, 705 }; 706 double start, cpu_time_used; 707 708 // Create the libCEED contexts 709 PetscScalar meter = 1e-2; // 1 meter in scaled length units 710 PetscScalar second = 1e-2; // 1 second in scaled time units 711 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 712 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 713 CeedScalar theta0 = 300.; // K 714 CeedScalar thetaC = -15.; // K 715 CeedScalar P0 = 1.e5; // Pa 716 CeedScalar N = 0.01; // 1/s 717 CeedScalar cv = 717.; // J/(kg K) 718 CeedScalar cp = 1004.; // J/(kg K) 719 CeedScalar g = 9.81; // m/s^2 720 CeedScalar lambda = -2./3.; // - 721 CeedScalar mu = 75.; // Pa s, dynamic viscosity 722 // mu = 75 is not physical for air, but is good for numerical stability 723 CeedScalar k = 0.02638; // W/(m K) 724 CeedScalar CtauS = 0.; // dimensionless 725 CeedScalar strong_form = 0.; // [0,1] 726 PetscScalar lx = 8000.; // m 727 PetscScalar ly = 8000.; // m 728 PetscScalar lz = 4000.; // m 729 CeedScalar rc = 1000.; // m (Radius of bubble) 730 PetscScalar resx = 1000.; // m (resolution in x) 731 PetscScalar resy = 1000.; // m (resolution in y) 732 PetscScalar resz = 1000.; // m (resolution in z) 733 PetscInt outputfreq = 10; // - 734 PetscInt contsteps = 0; // - 735 PetscInt degreeVol = 1; // - 736 PetscInt degreeSur = 1; // - 737 PetscInt qextraVol = 2; // - 738 PetscInt qextraSur = 2; // - 739 DMBoundaryType periodicity[] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 740 DM_BOUNDARY_NONE 741 }; 742 PetscReal center[3], dc_axis[3] = {0, 0, 0}; 743 744 ierr = PetscInitialize(&argc, &argv, NULL, help); 745 if (ierr) return ierr; 746 747 // Allocate PETSc context 748 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 749 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 750 751 // Parse command line options 752 comm = PETSC_COMM_WORLD; 753 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 754 NULL); CHKERRQ(ierr); 755 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 756 NULL, ceedresource, ceedresource, 757 sizeof(ceedresource), NULL); CHKERRQ(ierr); 758 ierr = PetscOptionsBool("-test", "Run in test mode", 759 NULL, test=PETSC_FALSE, &test, NULL); CHKERRQ(ierr); 760 problemChoice = NS_DENSITY_CURRENT; 761 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 762 problemTypes, (PetscEnum)problemChoice, 763 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 764 problem = &problemOptions[problemChoice]; 765 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 766 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 767 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 768 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 769 NULL, implicit=PETSC_FALSE, &implicit, NULL); 770 CHKERRQ(ierr); 771 { 772 PetscInt len; 773 PetscBool flg; 774 ierr = PetscOptionsIntArray("-bc_wall", 775 "Use wall boundary conditions on this list of faces", 776 NULL, bc.walls, 777 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 778 &len), &flg); CHKERRQ(ierr); 779 if (flg) bc.nwall = len; 780 for (PetscInt j=0; j<3; j++) { 781 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 782 ierr = PetscOptionsIntArray(flags[j], 783 "Use slip boundary conditions on this list of faces", 784 NULL, bc.slips[j], 785 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 786 &len), &flg); 787 CHKERRQ(ierr); 788 if (flg) bc.nslip[j] = len; 789 } 790 } 791 ierr = PetscOptionsInt("-viz_refine", 792 "Regular refinement levels for visualization", 793 NULL, viz_refine, &viz_refine, NULL); 794 CHKERRQ(ierr); 795 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 796 NULL, meter, &meter, NULL); CHKERRQ(ierr); 797 meter = fabs(meter); 798 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 799 NULL, second, &second, NULL); CHKERRQ(ierr); 800 second = fabs(second); 801 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 802 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 803 kilogram = fabs(kilogram); 804 ierr = PetscOptionsScalar("-units_Kelvin", 805 "1 Kelvin in scaled temperature units", 806 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 807 Kelvin = fabs(Kelvin); 808 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 809 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 810 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 811 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 812 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 813 NULL, P0, &P0, NULL); CHKERRQ(ierr); 814 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 815 NULL, N, &N, NULL); CHKERRQ(ierr); 816 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 817 NULL, cv, &cv, NULL); CHKERRQ(ierr); 818 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 819 NULL, cp, &cp, NULL); CHKERRQ(ierr); 820 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 821 NULL, g, &g, NULL); CHKERRQ(ierr); 822 ierr = PetscOptionsScalar("-lambda", 823 "Stokes hypothesis second viscosity coefficient", 824 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 825 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 826 NULL, mu, &mu, NULL); CHKERRQ(ierr); 827 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 828 NULL, k, &k, NULL); CHKERRQ(ierr); 829 ierr = PetscOptionsScalar("-CtauS", 830 "Scale coefficient for tau (nondimensional)", 831 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 832 ierr = PetscOptionsScalar("-strong_form", 833 "Strong (1) or weak/integrated by parts (0) advection residual", 834 NULL, strong_form, &strong_form, NULL); 835 CHKERRQ(ierr); 836 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 837 NULL, lx, &lx, NULL); CHKERRQ(ierr); 838 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 839 NULL, ly, &ly, NULL); CHKERRQ(ierr); 840 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 841 NULL, lz, &lz, NULL); CHKERRQ(ierr); 842 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 843 NULL, rc, &rc, NULL); CHKERRQ(ierr); 844 ierr = PetscOptionsScalar("-resx","Target resolution in x", 845 NULL, resx, &resx, NULL); CHKERRQ(ierr); 846 ierr = PetscOptionsScalar("-resy","Target resolution in y", 847 NULL, resy, &resy, NULL); CHKERRQ(ierr); 848 ierr = PetscOptionsScalar("-resz","Target resolution in z", 849 NULL, resz, &resz, NULL); CHKERRQ(ierr); 850 PetscInt n = problem->dim; 851 ierr = PetscOptionsEnumArray("-periodicity", "Periodicity per direction", 852 NULL, DMBoundaryTypes, (PetscEnum *)periodicity, 853 &n, NULL); CHKERRQ(ierr); 854 n = problem->dim; 855 center[0] = 0.5 * lx; 856 center[1] = 0.5 * ly; 857 center[2] = 0.5 * lz; 858 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 859 NULL, center, &n, NULL); CHKERRQ(ierr); 860 n = problem->dim; 861 ierr = PetscOptionsRealArray("-dc_axis", 862 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 863 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 864 { 865 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 866 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 867 if (norm > 0) { 868 for (int i=0; i<3; i++) dc_axis[i] /= norm; 869 } 870 } 871 ierr = PetscOptionsInt("-output_freq", 872 "Frequency of output, in number of steps", 873 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 874 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 875 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 876 ierr = PetscOptionsInt("-degree_volume", "Polynomial degree of finite elements for the Volume", 877 NULL, degreeVol, °reeVol, NULL); CHKERRQ(ierr); 878 ierr = PetscOptionsInt("-qextra_volume", "Number of extra quadrature points for the Volume", 879 NULL, qextraVol, &qextraVol, NULL); CHKERRQ(ierr); 880 PetscStrncpy(user->outputfolder, ".", 2); 881 ierr = PetscOptionsString("-of", "Output folder", 882 NULL, user->outputfolder, user->outputfolder, 883 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 884 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 885 886 // Define derived units 887 Pascal = kilogram / (meter * PetscSqr(second)); 888 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 889 mpersquareds = meter / PetscSqr(second); 890 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 891 kgpercubicm = kilogram / pow(meter,3); 892 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 893 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 894 895 // Scale variables to desired units 896 theta0 *= Kelvin; 897 thetaC *= Kelvin; 898 P0 *= Pascal; 899 N *= (1./second); 900 cv *= JperkgK; 901 cp *= JperkgK; 902 Rd = cp - cv; 903 g *= mpersquareds; 904 mu *= Pascal * second; 905 k *= WpermK; 906 lx = fabs(lx) * meter; 907 ly = fabs(ly) * meter; 908 lz = fabs(lz) * meter; 909 rc = fabs(rc) * meter; 910 resx = fabs(resx) * meter; 911 resy = fabs(resy) * meter; 912 resz = fabs(resz) * meter; 913 for (int i=0; i<3; i++) center[i] *= meter; 914 915 const CeedInt dim = problem->dim, ncompx = problem->dim, 916 qdatasizeVol = problem->qdatasizeVol, 917 qdatasizeSur = problem->qdatasizeSur; 918 // Set up the libCEED context 919 struct SetupContext_ ctxSetup = { 920 .theta0 = theta0, 921 .thetaC = thetaC, 922 .P0 = P0, 923 .N = N, 924 .cv = cv, 925 .cp = cp, 926 .Rd = Rd, 927 .g = g, 928 .rc = rc, 929 .lx = lx, 930 .ly = ly, 931 .lz = lz, 932 .periodicity0 = periodicity[0], 933 .periodicity1 = periodicity[1], 934 .periodicity2 = periodicity[2], 935 .center[0] = center[0], 936 .center[1] = center[1], 937 .center[2] = center[2], 938 .dc_axis[0] = dc_axis[0], 939 .dc_axis[1] = dc_axis[1], 940 .dc_axis[2] = dc_axis[2], 941 .time = 0, 942 }; 943 944 { 945 const PetscReal scale[3] = {lx, ly, lz}; 946 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 947 periodicity, PETSC_TRUE, &dm); 948 CHKERRQ(ierr); 949 } 950 if (1) { 951 DM dmDist = NULL; 952 PetscPartitioner part; 953 954 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 955 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 956 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 957 if (dmDist) { 958 ierr = DMDestroy(&dm); CHKERRQ(ierr); 959 dm = dmDist; 960 } 961 } 962 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 963 964 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 965 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 966 ierr = SetUpDM(dm, problem, degreeVol, &bc, &ctxSetup); CHKERRQ(ierr); 967 if (!test) { 968 ierr = PetscPrintf(PETSC_COMM_WORLD, 969 "Degree of Volumetric FEM Space: %D\n", 970 degreeVol); CHKERRQ(ierr); 971 } 972 dmviz = NULL; 973 interpviz = NULL; 974 if (viz_refine) { 975 DM dmhierarchy[viz_refine+1]; 976 977 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 978 dmhierarchy[0] = dm; 979 for (PetscInt i = 0, d = degreeVol; i < viz_refine; i++) { 980 Mat interp_next; 981 982 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 983 CHKERRQ(ierr); 984 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 985 d = (d + 1) / 2; 986 if (i + 1 == viz_refine) d = 1; 987 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 988 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 989 &interp_next, NULL); CHKERRQ(ierr); 990 if (!i) interpviz = interp_next; 991 else { 992 Mat C; 993 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 994 PETSC_DECIDE, &C); CHKERRQ(ierr); 995 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 996 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 997 interpviz = C; 998 } 999 } 1000 for (PetscInt i=1; i<viz_refine; i++) { 1001 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1002 } 1003 dmviz = dmhierarchy[viz_refine]; 1004 } 1005 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1006 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1007 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1008 lnodes /= ncompq; 1009 1010 { 1011 // Print grid information 1012 CeedInt gdofs, odofs; 1013 int comm_size; 1014 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1015 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1016 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1017 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1018 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1019 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1020 if (!test) { 1021 ierr = PetscPrintf(comm, "Global FEM dofs: %D (%D owned) on %d ranks\n", 1022 gdofs, odofs, comm_size); CHKERRQ(ierr); 1023 ierr = PetscPrintf(comm, "Local FEM nodes: %D\n", lnodes); CHKERRQ(ierr); 1024 ierr = PetscPrintf(comm, "dm_plex_box_faces: %s\n", box_faces_str); 1025 CHKERRQ(ierr); 1026 } 1027 1028 } 1029 1030 // Set up global mass vector 1031 ierr = VecDuplicate(Q,&user->M); CHKERRQ(ierr); 1032 1033 // Set up CEED 1034 // CEED Bases 1035 CeedInit(ceedresource, &ceed); 1036 numP_Vol = degreeVol + 1; 1037 numQ_Vol = numP_Vol + qextraVol; 1038 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP_Vol, numQ_Vol, CEED_GAUSS, 1039 &basisqVol); 1040 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ_Vol, CEED_GAUSS, 1041 &basisxVol); 1042 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP_Vol, 1043 CEED_GAUSS_LOBATTO, &basisxcVol); 1044 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1045 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1046 CHKERRQ(ierr); 1047 1048 // CEED Restrictions 1049 // Restrictions on the Volume 1050 ierr = GetRestrictionForDomain(ceed, dm, ncompx, dim, 0, 0, 0, numP_Vol, numQ_Vol, 1051 qdatasizeVol, &restrictqVol, &restrictxVol, 1052 &restrictqdiVol); CHKERRQ(ierr); 1053 1054 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1055 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1056 1057 // Create the CEED vectors that will be needed in setup 1058 CeedInt NqptsVol; 1059 CeedBasisGetNumQuadraturePoints(basisqVol, &NqptsVol); 1060 CeedElemRestrictionGetNumElements(restrictqVol, &localNelemVol); 1061 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1062 CeedElemRestrictionCreateVector(restrictqVol, &q0ceedVol, NULL); 1063 1064 // Create the Q-function that builds the quadrature data for the NS operator 1065 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1066 &qf_setupVol); 1067 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1068 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1069 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1070 1071 // Create the Q-function that sets the ICs of the operator 1072 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1073 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1074 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1075 1076 qf_rhsVol = NULL; 1077 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1078 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1079 problem->applyVol_rhs_loc, &qf_rhsVol); 1080 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1081 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1082 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1083 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1084 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1085 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1086 } 1087 1088 qf_ifunctionVol = NULL; 1089 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1090 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1091 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1092 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1093 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1094 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1095 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1096 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1097 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1098 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1099 } 1100 1101 // Create the operator that builds the quadrature data for the NS operator 1102 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1103 CeedOperatorSetField(op_setupVol, "dx", restrictxVol, basisxVol, CEED_VECTOR_ACTIVE); 1104 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1105 basisxVol, CEED_VECTOR_NONE); 1106 CeedOperatorSetField(op_setupVol, "qdata", restrictqdiVol, 1107 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1108 1109 // Create the operator that sets the ICs 1110 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1111 CeedOperatorSetField(op_ics, "x", restrictxVol, basisxcVol, CEED_VECTOR_ACTIVE); 1112 CeedOperatorSetField(op_ics, "q0", restrictqVol, 1113 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1114 1115 CeedElemRestrictionCreateVector(restrictqVol, &user->qceed, NULL); 1116 CeedElemRestrictionCreateVector(restrictqVol, &user->qdotceed, NULL); 1117 CeedElemRestrictionCreateVector(restrictqVol, &user->gceed, NULL); 1118 1119 if (qf_rhsVol) { // Create the RHS physics operator 1120 CeedOperator op; 1121 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1122 CeedOperatorSetField(op, "q", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1123 CeedOperatorSetField(op, "dq", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1124 CeedOperatorSetField(op, "qdata", restrictqdiVol, 1125 CEED_BASIS_COLLOCATED, qdata); 1126 CeedOperatorSetField(op, "x", restrictxVol, basisxVol, xcorners); 1127 CeedOperatorSetField(op, "v", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1128 CeedOperatorSetField(op, "dv", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1129 user->op_rhs = op; 1130 } 1131 1132 if (qf_ifunctionVol) { // Create the IFunction operator 1133 CeedOperator op; 1134 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1135 CeedOperatorSetField(op, "q", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1136 CeedOperatorSetField(op, "dq", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1137 CeedOperatorSetField(op, "qdot", restrictqVol, basisqVol, user->qdotceed); 1138 CeedOperatorSetField(op, "qdata", restrictqdiVol, 1139 CEED_BASIS_COLLOCATED, qdata); 1140 CeedOperatorSetField(op, "x", restrictxVol, basisxVol, xcorners); 1141 CeedOperatorSetField(op, "v", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1142 CeedOperatorSetField(op, "dv", restrictqVol, basisqVol, CEED_VECTOR_ACTIVE); 1143 user->op_ifunction = op; 1144 } 1145 1146 //**************************************************************************************// 1147 // Add boundary Integral (TODO: create a function for all faces) 1148 //--------------------------------------------------------------------------------------// 1149 // Set up CEED for the boundaries 1150 // CEED bases 1151 CeedInt height = 1; 1152 CeedInt dimSur = dim - height; 1153 numP_Sur = degreeSur + 1; 1154 numQ_Sur = numP_Sur + qextraSur; 1155 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, CEED_GAUSS, 1156 &basisqSur); 1157 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1158 &basisxSur); 1159 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1160 CEED_GAUSS_LOBATTO, &basisxcSur); 1161 // CEED Restrictions 1162 // Restriction on one face 1163 DMLabel domainLabel; 1164 ierr = DMGetLabel(dm, "Face Sets", &domainLabel); CHKERRQ(ierr); 1165 ierr = GetRestrictionForDomain(ceed, dm, ncompx, dimSur, height, domainLabel, 2, numP_Sur, 1166 numQ_Sur, qdatasizeSur, &restrictqSur, &restrictxSur, 1167 &restrictqdiSur); CHKERRQ(ierr); 1168 1169 // Create the CEED vectors that will be needed in setup 1170 CeedInt NqptsSur; 1171 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1172 CeedElemRestrictionGetNumElements(restrictqSur, &localNelemSur); 1173 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur*NqptsSur, &qdataSur); 1174 1175 // Create the Q-function that builds the quadrature data for the NS operator 1176 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1177 &qf_setupSur); 1178 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1179 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1180 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1181 1182 1183 qf_rhsSur = NULL; 1184 if (problem->applySur_rhs) { // Create the Q-function that defines the action of the RHS operator 1185 CeedQFunctionCreateInterior(ceed, 1, problem->applySur_rhs, 1186 problem->applySur_rhs_loc, &qf_rhsSur); 1187 CeedQFunctionAddInput(qf_rhsSur, "q", ncompq, CEED_EVAL_INTERP); 1188 CeedQFunctionAddInput(qf_rhsSur, "dq", ncompq*dim, CEED_EVAL_GRAD); //assumed volumetric elements 1189 CeedQFunctionAddInput(qf_rhsSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1190 CeedQFunctionAddInput(qf_rhsSur, "x", ncompx, CEED_EVAL_INTERP); 1191 CeedQFunctionAddInput(qf_rhsSur, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1192 CeedQFunctionAddOutput(qf_rhsSur, "v", ncompq, CEED_EVAL_INTERP); 1193 } 1194 1195 qf_ifunctionSur = NULL; 1196 if (problem->applySur_ifunction) { // Create the Q-function that defines the action of the IFunction 1197 CeedQFunctionCreateInterior(ceed, 1, problem->applySur_ifunction, 1198 problem->applySur_ifunction_loc, &qf_ifunctionSur); 1199 CeedQFunctionAddInput(qf_ifunctionSur, "q", ncompq, CEED_EVAL_INTERP); 1200 CeedQFunctionAddInput(qf_ifunctionSur, "dq", ncompq*dim, CEED_EVAL_GRAD); //assumed volumetric elements 1201 CeedQFunctionAddInput(qf_ifunctionSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1202 CeedQFunctionAddInput(qf_ifunctionSur, "x", ncompx, CEED_EVAL_INTERP); 1203 CeedQFunctionAddInput(qf_rhsSur, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1204 CeedQFunctionAddOutput(qf_ifunctionSur, "v", ncompq, CEED_EVAL_INTERP); 1205 } 1206 1207 // Create the operator that builds the quadrature data for the NS operator 1208 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur); 1209 CeedOperatorSetField(op_setupSur, "dx", restrictxSur, basisxSur, CEED_VECTOR_ACTIVE); 1210 CeedOperatorSetField(op_setupSur, "weight", CEED_ELEMRESTRICTION_NONE, 1211 basisxSur, CEED_VECTOR_NONE); 1212 CeedOperatorSetField(op_setupSur, "qdataSur", restrictqdiSur, 1213 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1214 1215 1216 if (qf_rhsSur) { // Create the RHS physics operator 1217 CeedOperator op; 1218 CeedOperatorCreate(ceed, qf_rhsSur, NULL, NULL, &op); 1219 CeedOperatorSetField(op, "q", restrictqSur, basisqSur, CEED_VECTOR_ACTIVE); 1220 CeedOperatorSetField(op, "dq", restrictqSur, basisqSur, CEED_VECTOR_ACTIVE); 1221 CeedOperatorSetField(op, "qdataSur", restrictqdiSur, 1222 CEED_BASIS_COLLOCATED, qdataSur); 1223 CeedOperatorSetField(op, "x", restrictxSur, basisxSur, xcorners); 1224 CeedOperatorSetField(op, "v", restrictqSur, basisqSur, CEED_VECTOR_ACTIVE); 1225 user->op_rhs_sur = op; 1226 } 1227 1228 if (qf_ifunctionSur) { // Create the IFunction operator 1229 CeedOperator op; 1230 CeedOperatorCreate(ceed, qf_ifunctionSur, NULL, NULL, &op); 1231 CeedOperatorSetField(op, "q", restrictqSur, basisqSur, CEED_VECTOR_ACTIVE); 1232 CeedOperatorSetField(op, "dq", restrictqSur, basisqSur, CEED_VECTOR_ACTIVE); 1233 CeedOperatorSetField(op, "qdataSur", restrictqdiSur, 1234 CEED_BASIS_COLLOCATED, qdataSur); 1235 CeedOperatorSetField(op, "x", restrictxSur, basisxSur, xcorners); 1236 CeedOperatorSetField(op, "v", restrictqSur, basisqSur, CEED_VECTOR_ACTIVE); 1237 user->op_ifunction_sur = op; 1238 } 1239 // Composite Operaters 1240 if (user->op_ifunction_vol) { 1241 if (user->op_ifunction_sur) { 1242 // Composite Operators for the IFunction 1243 CeedCompositeOperatorCreate(ceed, &user->op_ifunction); 1244 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_vol); 1245 CeedCompositeOperatorAddSub(user->op_ifunction, user->op_ifunction_sur); 1246 } else { 1247 user->op_ifunction = user->op_ifunction_vol; 1248 } 1249 } 1250 if (user->op_rhs_vol) { 1251 if (user->op_rhs_sur) { 1252 // Composite Operators for the RHS 1253 CeedCompositeOperatorCreate(ceed, &user->op_rhs); 1254 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_vol); 1255 CeedCompositeOperatorAddSub(user->op_rhs, user->op_rhs_sur); 1256 } else { 1257 user->op_rhs = user->op_rhs_vol; 1258 } 1259 } 1260 //**************************************************************************************// 1261 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1262 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1263 struct Advection2dContext_ ctxAdvection2d = { 1264 .CtauS = CtauS, 1265 .strong_form = strong_form, 1266 .stabilization = stab, 1267 }; 1268 switch (problemChoice) { 1269 case NS_DENSITY_CURRENT: 1270 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1271 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1272 sizeof ctxNS); 1273 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxNS, sizeof ctxNS); 1274 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxNS, 1275 sizeof ctxNS); 1276 break; 1277 case NS_ADVECTION: 1278 case NS_ADVECTION2D: 1279 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1280 sizeof ctxAdvection2d); 1281 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1282 sizeof ctxAdvection2d); 1283 if (qf_rhsSur) CeedQFunctionSetContext(qf_rhsSur, &ctxAdvection2d, 1284 sizeof ctxAdvection2d); 1285 if (qf_ifunctionSur) CeedQFunctionSetContext(qf_ifunctionSur, &ctxAdvection2d, 1286 sizeof ctxAdvection2d); 1287 } 1288 1289 // Set up PETSc context 1290 // Set up units structure 1291 units->meter = meter; 1292 units->kilogram = kilogram; 1293 units->second = second; 1294 units->Kelvin = Kelvin; 1295 units->Pascal = Pascal; 1296 units->JperkgK = JperkgK; 1297 units->mpersquareds = mpersquareds; 1298 units->WpermK = WpermK; 1299 units->kgpercubicm = kgpercubicm; 1300 units->kgpersquaredms = kgpersquaredms; 1301 units->Joulepercubicm = Joulepercubicm; 1302 1303 // Set up user structure 1304 user->comm = comm; 1305 user->outputfreq = outputfreq; 1306 user->contsteps = contsteps; 1307 user->units = units; 1308 user->dm = dm; 1309 user->dmviz = dmviz; 1310 user->interpviz = interpviz; 1311 user->ceed = ceed; 1312 1313 // Calculate qdata and ICs 1314 // Set up state global and local vectors 1315 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1316 1317 ierr = VectorPlacePetscVec(q0ceedVol, Qloc); CHKERRQ(ierr); 1318 1319 // Apply Setup Ceed Operators 1320 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1321 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1322 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictqVol, basisqVol, restrictqdiVol, qdata, 1323 user->M); CHKERRQ(ierr); 1324 1325 ierr = ICs_PetscMultiplicity(op_ics, xcorners, q0ceedVol, dm, Qloc, Q, restrictqVol, 1326 &ctxSetup, 0.0); 1327 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1328 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1329 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1330 // slower execution. 1331 Vec Qbc; 1332 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1333 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1334 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1335 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1336 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1337 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1338 ierr = PetscObjectComposeFunction((PetscObject)dm, 1339 "DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_NS); CHKERRQ(ierr); 1340 } 1341 1342 MPI_Comm_rank(comm, &rank); 1343 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1344 // Gather initial Q values 1345 // In case of continuation of simulation, set up initial values from binary file 1346 if (contsteps) { // continue from existent solution 1347 PetscViewer viewer; 1348 char filepath[PETSC_MAX_PATH_LEN]; 1349 // Read input 1350 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1351 user->outputfolder); 1352 CHKERRQ(ierr); 1353 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1354 CHKERRQ(ierr); 1355 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1356 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1357 } else { 1358 //ierr = DMLocalToGlobal(dm, Qloc, INSERT_VALUES, Q);CHKERRQ(ierr); 1359 } 1360 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1361 1362 // Create and setup TS 1363 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1364 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1365 if (implicit) { 1366 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1367 if (user->op_ifunction) { 1368 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1369 } else { // Implicit integrators can fall back to using an RHSFunction 1370 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1371 } 1372 } else { 1373 if (!user->op_rhs) SETERRQ(comm,PETSC_ERR_ARG_NULL, 1374 "Problem does not provide RHSFunction"); 1375 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1376 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1377 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1378 } 1379 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1380 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1381 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1382 if (test) {ierr = TSSetMaxSteps(ts, 1); CHKERRQ(ierr);} 1383 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1384 ierr = TSAdaptSetStepLimits(adapt, 1385 1.e-12 * units->second, 1386 1.e2 * units->second); CHKERRQ(ierr); 1387 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1388 if (!contsteps) { // print initial condition 1389 if (!test) { 1390 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1391 } 1392 } else { // continue from time of last output 1393 PetscReal time; 1394 PetscInt count; 1395 PetscViewer viewer; 1396 char filepath[PETSC_MAX_PATH_LEN]; 1397 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1398 user->outputfolder); CHKERRQ(ierr); 1399 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1400 CHKERRQ(ierr); 1401 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1402 CHKERRQ(ierr); 1403 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1404 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1405 } 1406 if (!test) { 1407 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1408 } 1409 1410 // Solve 1411 start = MPI_Wtime(); 1412 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1413 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1414 cpu_time_used = MPI_Wtime() - start; 1415 ierr = TSGetSolveTime(ts,&ftime); CHKERRQ(ierr); 1416 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1417 comm); CHKERRQ(ierr); 1418 if (!test) { 1419 ierr = PetscPrintf(PETSC_COMM_WORLD, 1420 "Time taken for solution: %g\n", 1421 (double)cpu_time_used); CHKERRQ(ierr); 1422 } 1423 1424 // Get error 1425 if (problem->non_zero_time && !test) { 1426 Vec Qexact, Qexactloc; 1427 PetscReal norm; 1428 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1429 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1430 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1431 1432 ierr = ICs_PetscMultiplicity(op_ics, xcorners, q0ceedVol, dm, Qexactloc, Qexact, 1433 restrictqVol, &ctxSetup, ftime); CHKERRQ(ierr); 1434 1435 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1436 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1437 CeedVectorDestroy(&q0ceedVol); 1438 ierr = PetscPrintf(PETSC_COMM_WORLD, 1439 "Max Error: %g\n", 1440 (double)norm); CHKERRQ(ierr); 1441 } 1442 1443 // Output Statistics 1444 ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr); 1445 if (!test) { 1446 ierr = PetscPrintf(PETSC_COMM_WORLD, 1447 "Time integrator took %D time steps to reach final time %g\n", 1448 steps,(double)ftime); CHKERRQ(ierr); 1449 } 1450 1451 // Clean up libCEED 1452 CeedVectorDestroy(&qdata); 1453 CeedVectorDestroy(&user->qceed); 1454 CeedVectorDestroy(&user->qdotceed); 1455 CeedVectorDestroy(&user->gceed); 1456 CeedVectorDestroy(&xcorners); 1457 CeedBasisDestroy(&basisqVol); 1458 CeedBasisDestroy(&basisxVol); 1459 CeedBasisDestroy(&basisxcVol); 1460 CeedElemRestrictionDestroy(&restrictqVol); 1461 CeedElemRestrictionDestroy(&restrictxVol); 1462 CeedElemRestrictionDestroy(&restrictqdiVol); 1463 CeedQFunctionDestroy(&qf_setupVol); 1464 CeedQFunctionDestroy(&qf_ics); 1465 CeedQFunctionDestroy(&qf_rhsVol); 1466 CeedQFunctionDestroy(&qf_ifunctionVol); 1467 CeedOperatorDestroy(&op_setupVol); 1468 CeedOperatorDestroy(&op_ics); 1469 CeedOperatorDestroy(&user->op_rhs_vol); 1470 CeedOperatorDestroy(&user->op_ifunction_vol); 1471 CeedDestroy(&ceed); 1472 1473 CeedVectorDestroy(&qdataSur); 1474 CeedBasisDestroy(&basisqSur); 1475 CeedBasisDestroy(&basisxSur); 1476 CeedBasisDestroy(&basisxcSur); 1477 CeedElemRestrictionDestroy(&restrictqSur); 1478 CeedElemRestrictionDestroy(&restrictxSur); 1479 CeedElemRestrictionDestroy(&restrictqdiSur); 1480 CeedQFunctionDestroy(&qf_setupSur); 1481 CeedQFunctionDestroy(&qf_rhsSur); 1482 CeedQFunctionDestroy(&qf_ifunctionSur); 1483 CeedOperatorDestroy(&op_setupSur); 1484 CeedOperatorDestroy(&user->op_rhs_sur); 1485 CeedOperatorDestroy(&user->op_ifunction_sur); 1486 1487 CeedOperatorDestroy(&user->op_rhs); 1488 CeedOperatorDestroy(&user->op_ifunction); 1489 1490 // Clean up PETSc 1491 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1492 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1493 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1494 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1495 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1496 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1497 ierr = PetscFree(units); CHKERRQ(ierr); 1498 ierr = PetscFree(user); CHKERRQ(ierr); 1499 return PetscFinalize(); 1500 } 1501 1502