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