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