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