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