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