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