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