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