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