1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 // libCEED + PETSc Example: Navier-Stokes 18 // 19 // This example demonstrates a simple usage of libCEED with PETSc to solve a 20 // Navier-Stokes problem. 21 // 22 // The code is intentionally "raw", using only low-level communication 23 // primitives. 24 // 25 // Build with: 26 // 27 // make [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] navierstokes 28 // 29 // Sample runs: 30 // 31 // ./navierstokes -ceed /cpu/self -problem density_current -degree 1 32 // ./navierstokes -ceed /gpu/occa -problem advection -degree 1 33 // 34 //TESTARGS -ceed {ceed_resource} -test explicit -degree 3 -dm_plex_box_faces 1,1,2 -units_kilogram 1e-9 -lx 125 -ly 125 -lz 250 -center 62.5,62.5,187.5 -rc 100. -thetaC -35. -ksp_atol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 35 //TESTARGS -ceed {ceed_resource} -test implicit_stab_none -degree 3 -dm_plex_box_faces 1,1,2 -units_kilogram 1e-9 -lx 125 -ly 125 -lz 250 -center 62.5,62.5,187.5 -rc 100. -thetaC -35. -ksp_atol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 -implicit -ts_type alpha 36 //TESTARGS -ceed {ceed_resource} -test implicit_stab_supg -degree 3 -dm_plex_box_faces 1,1,2 -units_kilogram 1e-9 -lx 125 -ly 125 -lz 250 -center 62.5,62.5,187.5 -rc 100. -thetaC -35. -ksp_atol 1e-4 -ksp_rtol 1e-3 -ksp_type bcgs -snes_atol 1e-3 -snes_lag_jacobian 100 -snes_lag_jacobian_persists -snes_mf_operator -ts_dt 1e-3 -implicit -ts_type alpha -stab supg 37 38 /// @file 39 /// Navier-Stokes example using PETSc 40 41 const char help[] = "Solve Navier-Stokes using PETSc and libCEED\n"; 42 43 #include <petscts.h> 44 #include <petscdmplex.h> 45 #include <ceed.h> 46 #include <stdbool.h> 47 #include <petscsys.h> 48 #include "common.h" 49 #include "setup-boundary.h" 50 #include "advection.h" 51 #include "advection2d.h" 52 #include "densitycurrent.h" 53 54 #if PETSC_VERSION_LT(3,14,0) 55 # define DMPlexGetClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexGetClosureIndices(a,b,c,d,f,g,i) 56 # define DMPlexRestoreClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexRestoreClosureIndices(a,b,c,d,f,g,i) 57 #endif 58 59 // MemType Options 60 static const char *const memTypes[] = { 61 "host", 62 "device", 63 "memType", "CEED_MEM_", NULL 64 }; 65 66 // Problem Options 67 typedef enum { 68 NS_DENSITY_CURRENT = 0, 69 NS_ADVECTION = 1, 70 NS_ADVECTION2D = 2, 71 } problemType; 72 static const char *const problemTypes[] = { 73 "density_current", 74 "advection", 75 "advection2d", 76 "problemType", "NS_", NULL 77 }; 78 79 typedef enum { 80 STAB_NONE = 0, 81 STAB_SU = 1, // Streamline Upwind 82 STAB_SUPG = 2, // Streamline Upwind Petrov-Galerkin 83 } StabilizationType; 84 static const char *const StabilizationTypes[] = { 85 "none", 86 "SU", 87 "SUPG", 88 "StabilizationType", "STAB_", NULL 89 }; 90 91 // Test Options 92 typedef enum { 93 TEST_NONE = 0, // Non test mode 94 TEST_EXPLICIT = 1, // Explicit test 95 TEST_IMPLICIT_STAB_NONE = 2, // Implicit test no stab 96 TEST_IMPLICIT_STAB_SUPG = 3, // Implicit test supg stab 97 } testType; 98 static const char *const testTypes[] = { 99 "none", 100 "explicit", 101 "implicit_stab_none", 102 "implicit_stab_supg", 103 "testType", "TEST_", NULL 104 }; 105 106 // Tests specific data 107 typedef struct { 108 PetscScalar testtol; 109 const char *filepath; 110 } testData; 111 112 testData testOptions[] = { 113 [TEST_NONE] = { 114 .testtol = 0., 115 .filepath = NULL 116 }, 117 [TEST_EXPLICIT] = { 118 .testtol = 1E-5, 119 .filepath = "examples/fluids/tests-output/fluids-navierstokes-explicit.bin" 120 }, 121 [TEST_IMPLICIT_STAB_NONE] = { 122 .testtol = 5E-4, 123 .filepath = "examples/fluids/tests-output/fluids-navierstokes-implicit-stab-none.bin" 124 }, 125 [TEST_IMPLICIT_STAB_SUPG] = { 126 .testtol = 5E-4, 127 .filepath = "examples/fluids/tests-output/fluids-navierstokes-implicit-stab-supg.bin" 128 } 129 }; 130 131 // Problem specific data 132 typedef struct { 133 CeedInt dim, qdatasizeVol, qdatasizeSur; 134 CeedQFunctionUser setupVol, setupSur, ics, applyVol_rhs, applyVol_ifunction, 135 applyOut_rhs, applyOut_ifunction, applyIn_rhs, applyIn_ifunction; 136 PetscErrorCode (*bc)(PetscInt, PetscReal, const PetscReal[], PetscInt, 137 PetscScalar[], void *); 138 const char *setupVol_loc, *setupSur_loc, *ics_loc, *applyVol_rhs_loc, *applyVol_ifunction_loc, 139 *applyOut_rhs_loc, *applyOut_ifunction_loc, *applyIn_rhs_loc, *applyIn_ifunction_loc; 140 const bool non_zero_time; 141 } problemData; 142 143 problemData problemOptions[] = { 144 [NS_DENSITY_CURRENT] = { 145 .dim = 3, 146 .qdatasizeVol = 10, 147 .qdatasizeSur = 4, 148 .setupVol = Setup, 149 .setupVol_loc = Setup_loc, 150 .setupSur = SetupBoundary, 151 .setupSur_loc = SetupBoundary_loc, 152 .ics = ICsDC, 153 .ics_loc = ICsDC_loc, 154 .applyVol_rhs = DC, 155 .applyVol_rhs_loc = DC_loc, 156 //.applyOut_rhs = DC_Out, 157 //.applyOut_rhs_loc = DC_Out_loc, 158 //.applyIn_rhs = DC_In, 159 //.applyIn_rhs_loc = DC_In_loc, 160 .applyVol_ifunction = IFunction_DC, 161 .applyVol_ifunction_loc = IFunction_DC_loc, 162 //.applyOut_ifunction = IFunction_DC_Out, 163 //.applyOut_ifunction_loc = IFunction_DC_Out_loc, 164 //.applyIn_ifunction = IFunction_DC_In, 165 //.applyIn_ifunction_loc = IFunction_DC_In_loc, 166 .bc = Exact_DC, 167 .non_zero_time = PETSC_FALSE, 168 }, 169 [NS_ADVECTION] = { 170 .dim = 3, 171 .qdatasizeVol = 10, 172 .qdatasizeSur = 4, 173 .setupVol = Setup, 174 .setupVol_loc = Setup_loc, 175 .setupSur = SetupBoundary, 176 .setupSur_loc = SetupBoundary_loc, 177 .ics = ICsAdvection, 178 .ics_loc = ICsAdvection_loc, 179 .applyVol_rhs = Advection, 180 .applyVol_rhs_loc = Advection_loc, 181 .applyOut_rhs = Advection_Out, 182 .applyOut_rhs_loc = Advection_Out_loc, 183 .applyIn_rhs = Advection_In, 184 .applyIn_rhs_loc = Advection_In_loc, 185 .applyVol_ifunction = IFunction_Advection, 186 .applyVol_ifunction_loc = IFunction_Advection_loc, 187 .applyOut_ifunction = IFunction_Advection_Out, 188 .applyOut_ifunction_loc = IFunction_Advection_Out_loc, 189 .applyIn_ifunction = IFunction_Advection_In, 190 .applyIn_ifunction_loc = IFunction_Advection_In_loc, 191 .bc = Exact_Advection, 192 .non_zero_time = PETSC_FALSE, 193 }, 194 [NS_ADVECTION2D] = { 195 .dim = 2, 196 .qdatasizeVol = 5, 197 .qdatasizeSur = 3, 198 .setupVol = Setup2d, 199 .setupVol_loc = Setup2d_loc, 200 .setupSur = SetupBoundary2d, 201 .setupSur_loc = SetupBoundary2d_loc, 202 .ics = ICsAdvection2d, 203 .ics_loc = ICsAdvection2d_loc, 204 .applyVol_rhs = Advection2d, 205 .applyVol_rhs_loc = Advection2d_loc, 206 .applyOut_rhs = Advection2d_Out, 207 .applyOut_rhs_loc = Advection2d_Out_loc, 208 .applyIn_rhs = Advection2d_In, 209 .applyIn_rhs_loc = Advection2d_In_loc, 210 .applyVol_ifunction = IFunction_Advection2d, 211 .applyVol_ifunction_loc = IFunction_Advection2d_loc, 212 .applyOut_ifunction = IFunction_Advection2d_Out, 213 .applyOut_ifunction_loc = IFunction_Advection2d_Out_loc, 214 .applyIn_ifunction = IFunction_Advection2d_In, 215 .applyIn_ifunction_loc = IFunction_Advection2d_In_loc, 216 .bc = Exact_Advection2d, 217 .non_zero_time = PETSC_TRUE, 218 }, 219 }; 220 221 // PETSc user data 222 typedef struct User_ *User; 223 typedef struct Units_ *Units; 224 225 struct User_ { 226 MPI_Comm comm; 227 PetscInt outputfreq; 228 DM dm; 229 DM dmviz; 230 Mat interpviz; 231 Ceed ceed; 232 Units units; 233 CeedVector qceed, qdotceed, gceed; 234 CeedOperator op_rhs_vol, op_rhs, op_ifunction_vol, op_ifunction; 235 Vec M; 236 char outputfolder[PETSC_MAX_PATH_LEN]; 237 PetscInt contsteps; 238 }; 239 240 struct Units_ { 241 // fundamental units 242 PetscScalar meter; 243 PetscScalar kilogram; 244 PetscScalar second; 245 PetscScalar Kelvin; 246 // derived units 247 PetscScalar Pascal; 248 PetscScalar JperkgK; 249 PetscScalar mpersquareds; 250 PetscScalar WpermK; 251 PetscScalar kgpercubicm; 252 PetscScalar kgpersquaredms; 253 PetscScalar Joulepercubicm; 254 }; 255 256 typedef struct SimpleBC_ *SimpleBC; 257 struct SimpleBC_ { 258 PetscInt nwall, nslip[3], noutflow, ninflow; 259 PetscInt walls[6], slips[3][6], outflow[6], inflow[6]; 260 PetscBool userbc; 261 }; 262 263 // Essential BC dofs are encoded in closure indices as -(i+1). 264 static PetscInt Involute(PetscInt i) { 265 return i >= 0 ? i : -(i+1); 266 } 267 268 // Utility function to create local CEED restriction 269 static PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 270 CeedInt height, DMLabel domainLabel, CeedInt value, 271 CeedElemRestriction *Erestrict) { 272 273 PetscSection section; 274 PetscInt p, Nelem, Ndof, *erestrict, eoffset, nfields, dim, 275 depth; 276 DMLabel depthLabel; 277 IS depthIS, iterIS; 278 Vec Uloc; 279 const PetscInt *iterIndices; 280 PetscErrorCode ierr; 281 282 PetscFunctionBeginUser; 283 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 284 dim -= height; 285 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 286 ierr = PetscSectionGetNumFields(section, &nfields); CHKERRQ(ierr); 287 PetscInt ncomp[nfields], fieldoff[nfields+1]; 288 fieldoff[0] = 0; 289 for (PetscInt f=0; f<nfields; f++) { 290 ierr = PetscSectionGetFieldComponents(section, f, &ncomp[f]); CHKERRQ(ierr); 291 fieldoff[f+1] = fieldoff[f] + ncomp[f]; 292 } 293 294 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 295 ierr = DMPlexGetDepthLabel(dm, &depthLabel); CHKERRQ(ierr); 296 ierr = DMLabelGetStratumIS(depthLabel, depth - height, &depthIS); CHKERRQ(ierr); 297 if (domainLabel) { 298 IS domainIS; 299 ierr = DMLabelGetStratumIS(domainLabel, value, &domainIS); CHKERRQ(ierr); 300 ierr = ISIntersect(depthIS, domainIS, &iterIS); CHKERRQ(ierr); 301 ierr = ISDestroy(&domainIS); CHKERRQ(ierr); 302 ierr = ISDestroy(&depthIS); CHKERRQ(ierr); 303 } else { 304 iterIS = depthIS; 305 } 306 ierr = ISGetLocalSize(iterIS, &Nelem); CHKERRQ(ierr); 307 ierr = ISGetIndices(iterIS, &iterIndices); CHKERRQ(ierr); 308 ierr = PetscMalloc1(Nelem*PetscPowInt(P, dim), &erestrict); CHKERRQ(ierr); 309 for (p=0,eoffset=0; p<Nelem; p++) { 310 PetscInt c = iterIndices[p]; 311 PetscInt numindices, *indices, nnodes; 312 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, 313 &numindices, &indices, NULL, NULL); 314 CHKERRQ(ierr); 315 bool flip = false; 316 if (height > 0) { 317 PetscInt numCells, numFaces, start = -1; 318 const PetscInt *orients, *faces, *cells; 319 ierr = DMPlexGetSupport(dm, c, &cells); CHKERRQ(ierr); 320 ierr = DMPlexGetSupportSize(dm, c, &numCells); CHKERRQ(ierr); 321 if (numCells != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Expected one cell in support of exterior face, but got %D cells", numCells); 322 ierr = DMPlexGetCone(dm, cells[0], &faces); CHKERRQ(ierr); 323 ierr = DMPlexGetConeSize(dm, cells[0], &numFaces); CHKERRQ(ierr); 324 for (PetscInt i=0; i<numFaces; i++) {if (faces[i] == c) start = i;} 325 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Could not find face %D in cone of its support", c); 326 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients); CHKERRQ(ierr); 327 if (orients[start] < 0) flip = true; 328 } 329 if (numindices % fieldoff[nfields]) SETERRQ1(PETSC_COMM_SELF, 330 PETSC_ERR_ARG_INCOMP, "Number of closure indices not compatible with Cell %D", 331 c); 332 nnodes = numindices / fieldoff[nfields]; 333 for (PetscInt i=0; i<nnodes; i++) { 334 PetscInt ii = i; 335 if (flip) { 336 if (P == nnodes) ii = nnodes - 1 - i; 337 else if (P*P == nnodes) { 338 PetscInt row = i / P, col = i % P; 339 ii = row + col * P; 340 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for flipping point with %D nodes != P (%D) or P^2", nnodes, P); 341 } 342 // Check that indices are blocked by node and thus can be coalesced as a single field with 343 // fieldoff[nfields] = sum(ncomp) components. 344 for (PetscInt f=0; f<nfields; f++) { 345 for (PetscInt j=0; j<ncomp[f]; j++) { 346 if (Involute(indices[fieldoff[f]*nnodes + ii*ncomp[f] + j]) 347 != Involute(indices[ii*ncomp[0]]) + fieldoff[f] + j) 348 SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, 349 "Cell %D closure indices not interlaced for node %D field %D component %D", 350 c, ii, f, j); 351 } 352 } 353 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 354 PetscInt loc = Involute(indices[ii*ncomp[0]]); 355 erestrict[eoffset++] = loc; 356 } 357 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, 358 &numindices, &indices, NULL, NULL); 359 CHKERRQ(ierr); 360 } 361 if (eoffset != Nelem*PetscPowInt(P, dim)) 362 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 363 "ElemRestriction of size (%D,%D) initialized %D nodes", Nelem, 364 PetscPowInt(P, dim),eoffset); 365 ierr = ISRestoreIndices(iterIS, &iterIndices); CHKERRQ(ierr); 366 ierr = ISDestroy(&iterIS); CHKERRQ(ierr); 367 368 ierr = DMGetLocalVector(dm, &Uloc); CHKERRQ(ierr); 369 ierr = VecGetLocalSize(Uloc, &Ndof); CHKERRQ(ierr); 370 ierr = DMRestoreLocalVector(dm, &Uloc); CHKERRQ(ierr); 371 CeedElemRestrictionCreate(ceed, Nelem, PetscPowInt(P, dim), fieldoff[nfields], 372 1, Ndof, CEED_MEM_HOST, CEED_COPY_VALUES, erestrict, 373 Erestrict); 374 ierr = PetscFree(erestrict); CHKERRQ(ierr); 375 PetscFunctionReturn(0); 376 } 377 378 // Utility function to get Ceed Restriction for each domain 379 static PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 380 DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, CeedInt qdatasize, 381 CeedElemRestriction *restrictq, CeedElemRestriction *restrictx, 382 CeedElemRestriction *restrictqdi) { 383 384 DM dmcoord; 385 CeedInt dim, localNelem; 386 CeedInt Qdim; 387 PetscErrorCode ierr; 388 389 PetscFunctionBeginUser; 390 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 391 dim -= height; 392 Qdim = CeedIntPow(Q, dim); 393 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 394 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 395 CHKERRQ(ierr); 396 ierr = CreateRestrictionFromPlex(ceed, dm, P, height, domainLabel, value, restrictq); 397 CHKERRQ(ierr); 398 ierr = CreateRestrictionFromPlex(ceed, dmcoord, 2, height, domainLabel, value, restrictx); 399 CHKERRQ(ierr); 400 CeedElemRestrictionGetNumElements(*restrictq, &localNelem); 401 CeedElemRestrictionCreateStrided(ceed, localNelem, Qdim, 402 qdatasize, qdatasize*localNelem*Qdim, 403 CEED_STRIDES_BACKEND, restrictqdi); 404 PetscFunctionReturn(0); 405 } 406 407 // Utility function to create CEED Composite Operator for the entire domain 408 static PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, CeedOperator op_applyVol, 409 CeedQFunction qf_applyOut, CeedQFunction qf_applyIn, CeedQFunction qf_setupSur, CeedInt height, 410 PetscInt nOut, PetscInt valueOut[6], PetscInt nIn, PetscInt valueIn[6], CeedInt numP_Sur, 411 CeedInt numQ_Sur, CeedInt qdatasizeSur, CeedInt NqptsSur, 412 CeedBasis basisxSur, CeedBasis basisqSur, CeedOperator *op_apply) { 413 414 CeedElemRestriction restrictxOut[6], restrictqOut[6], restrictqdiOut[6], 415 restrictxIn[6], restrictqIn[6], restrictqdiIn[6]; 416 PetscInt localNelemOut[6], localNelemIn[6]; 417 Vec Xloc; 418 CeedVector xcorners, qdataOut[6], qdataIn[6]; 419 CeedOperator op_setupOut[6], op_applyOut[6], op_setupIn[6], op_applyIn[6]; 420 DMLabel domainLabel; 421 PetscScalar *x; 422 PetscInt lsize; 423 PetscErrorCode ierr; 424 425 PetscFunctionBeginUser; 426 // Composite Operaters 427 CeedCompositeOperatorCreate(ceed, op_apply); 428 CeedCompositeOperatorAddSub(*op_apply, op_applyVol); // Apply a Sub-Operator for the volume 429 430 if (nOut || nIn) { 431 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 432 ierr = VecGetLocalSize(Xloc, &lsize); CHKERRQ(ierr); 433 ierr = CeedVectorCreate(ceed, lsize, &xcorners); CHKERRQ(ierr); 434 ierr = VecGetArray(Xloc, &x); CHKERRQ(ierr); 435 CeedVectorSetArray(xcorners, CEED_MEM_HOST, CEED_USE_POINTER, x); 436 ierr = DMGetLabel(dm, "Face Sets", &domainLabel); CHKERRQ(ierr); 437 438 // Create CEED Operator for each OutFlow faces 439 if (nOut) { 440 for (CeedInt i=0; i<nOut; i++) { 441 ierr = GetRestrictionForDomain(ceed, dm, height, domainLabel, valueOut[i], numP_Sur, 442 numQ_Sur, qdatasizeSur, &restrictqOut[i], &restrictxOut[i], 443 &restrictqdiOut[i]); CHKERRQ(ierr); 444 // Create the CEED vectors that will be needed in boundary setup 445 CeedElemRestrictionGetNumElements(restrictqOut[i], &localNelemOut[i]); 446 CeedVectorCreate(ceed, qdatasizeSur*localNelemOut[i]*NqptsSur, &qdataOut[i]); 447 // Create the operator that builds the quadrature data for the OutFlow operator 448 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupOut[i]); 449 CeedOperatorSetField(op_setupOut[i], "dx", restrictxOut[i], basisxSur, CEED_VECTOR_ACTIVE); 450 CeedOperatorSetField(op_setupOut[i], "weight", CEED_ELEMRESTRICTION_NONE, 451 basisxSur, CEED_VECTOR_NONE); 452 CeedOperatorSetField(op_setupOut[i], "qdataSur", restrictqdiOut[i], 453 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 454 // Create OutFlow operator 455 CeedOperatorCreate(ceed, qf_applyOut, NULL, NULL, &op_applyOut[i]); 456 CeedOperatorSetField(op_applyOut[i], "q", restrictqOut[i], basisqSur, CEED_VECTOR_ACTIVE); 457 CeedOperatorSetField(op_applyOut[i], "qdataSur", restrictqdiOut[i], 458 CEED_BASIS_COLLOCATED, qdataOut[i]); 459 CeedOperatorSetField(op_applyOut[i], "x", restrictxOut[i], basisxSur, xcorners); 460 CeedOperatorSetField(op_applyOut[i], "v", restrictqOut[i], basisqSur, CEED_VECTOR_ACTIVE); 461 // Apply CEED operator for OutFlow setup 462 CeedOperatorApply(op_setupOut[i], xcorners, qdataOut[i], CEED_REQUEST_IMMEDIATE); 463 // Apply Sub-Operator for OutFlow BCs 464 CeedCompositeOperatorAddSub(*op_apply, op_applyOut[i]); 465 } 466 } 467 // Create CEED Operator for each InFlow faces 468 if (nIn) { 469 for (CeedInt i=0; i<nIn; i++) { 470 ierr = GetRestrictionForDomain(ceed, dm, height, domainLabel, valueIn[i], numP_Sur, 471 numQ_Sur, qdatasizeSur, &restrictqIn[i], &restrictxIn[i], 472 &restrictqdiIn[i]); CHKERRQ(ierr); 473 // Create the CEED vectors that will be needed in boundary setup 474 CeedElemRestrictionGetNumElements(restrictqIn[i], &localNelemIn[i]); 475 CeedVectorCreate(ceed, qdatasizeSur*localNelemIn[i]*NqptsSur, &qdataIn[i]); 476 // Create the operator that builds the quadrature data for the InFlow operator 477 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupIn[i]); 478 CeedOperatorSetField(op_setupIn[i], "dx", restrictxIn[i], basisxSur, CEED_VECTOR_ACTIVE); 479 CeedOperatorSetField(op_setupIn[i], "weight", CEED_ELEMRESTRICTION_NONE, 480 basisxSur, CEED_VECTOR_NONE); 481 CeedOperatorSetField(op_setupIn[i], "qdataSur", restrictqdiIn[i], 482 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 483 // Create InFlow operator 484 CeedOperatorCreate(ceed, qf_applyIn, NULL, NULL, &op_applyIn[i]); 485 CeedOperatorSetField(op_applyIn[i], "q", restrictqIn[i], basisqSur, CEED_VECTOR_ACTIVE); 486 CeedOperatorSetField(op_applyIn[i], "qdataSur", restrictqdiIn[i], 487 CEED_BASIS_COLLOCATED, qdataIn[i]); 488 CeedOperatorSetField(op_applyIn[i], "x", restrictxIn[i], basisxSur, xcorners); 489 CeedOperatorSetField(op_applyIn[i], "v", restrictqIn[i], basisqSur, CEED_VECTOR_ACTIVE); 490 // Apply CEED operator for InFlow setup 491 CeedOperatorApply(op_setupIn[i], xcorners, qdataIn[i], CEED_REQUEST_IMMEDIATE); 492 // Apply Sub-Operator for InFlow BCs 493 CeedCompositeOperatorAddSub(*op_apply, op_applyIn[i]); 494 } 495 } 496 CeedVectorDestroy(&xcorners); 497 } 498 PetscFunctionReturn(0); 499 } 500 501 static int CreateVectorFromPetscVec(Ceed ceed, Vec p, CeedVector *v) { 502 PetscErrorCode ierr; 503 PetscInt m; 504 505 PetscFunctionBeginUser; 506 ierr = VecGetLocalSize(p, &m); CHKERRQ(ierr); 507 ierr = CeedVectorCreate(ceed, m, v); CHKERRQ(ierr); 508 PetscFunctionReturn(0); 509 } 510 511 static int VectorPlacePetscVec(CeedVector c, Vec p) { 512 PetscErrorCode ierr; 513 PetscInt mceed,mpetsc; 514 PetscScalar *a; 515 516 PetscFunctionBeginUser; 517 ierr = CeedVectorGetLength(c, &mceed); CHKERRQ(ierr); 518 ierr = VecGetLocalSize(p, &mpetsc); CHKERRQ(ierr); 519 if (mceed != mpetsc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, 520 "Cannot place PETSc Vec of length %D in CeedVector of length %D", 521 mpetsc, mceed); 522 ierr = VecGetArray(p, &a); CHKERRQ(ierr); 523 CeedVectorSetArray(c, CEED_MEM_HOST, CEED_USE_POINTER, a); 524 PetscFunctionReturn(0); 525 } 526 527 static PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 528 PetscBool insertEssential, Vec Qloc, PetscReal time, Vec faceGeomFVM, 529 Vec cellGeomFVM, Vec gradFVM) { 530 PetscErrorCode ierr; 531 Vec Qbc; 532 533 PetscFunctionBegin; 534 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 535 ierr = VecAXPY(Qloc, 1., Qbc); CHKERRQ(ierr); 536 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 537 PetscFunctionReturn(0); 538 } 539 540 // This is the RHS of the ODE, given as u_t = G(t,u) 541 // This function takes in a state vector Q and writes into G 542 static PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *userData) { 543 PetscErrorCode ierr; 544 User user = *(User *)userData; 545 PetscScalar *q, *g; 546 Vec Qloc, Gloc; 547 548 // Global-to-local 549 PetscFunctionBeginUser; 550 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 551 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 552 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 553 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 554 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 555 NULL, NULL, NULL); CHKERRQ(ierr); 556 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 557 558 // Ceed Vectors 559 ierr = VecGetArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 560 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 561 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); 562 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 563 564 // Apply CEED operator 565 CeedOperatorApply(user->op_rhs, user->qceed, user->gceed, 566 CEED_REQUEST_IMMEDIATE); 567 568 // Restore vectors 569 ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 570 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 571 572 ierr = VecZeroEntries(G); CHKERRQ(ierr); 573 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 574 575 // Inverse of the lumped mass matrix 576 ierr = VecPointwiseMult(G, G, user->M); // M is Minv 577 CHKERRQ(ierr); 578 579 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 580 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 584 static PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Qdot, Vec G, 585 void *userData) { 586 PetscErrorCode ierr; 587 User user = *(User *)userData; 588 const PetscScalar *q, *qdot; 589 PetscScalar *g; 590 Vec Qloc, Qdotloc, Gloc; 591 592 // Global-to-local 593 PetscFunctionBeginUser; 594 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 595 ierr = DMGetLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 596 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 597 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 598 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 599 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 600 NULL, NULL, NULL); CHKERRQ(ierr); 601 ierr = VecZeroEntries(Qdotloc); CHKERRQ(ierr); 602 ierr = DMGlobalToLocal(user->dm, Qdot, INSERT_VALUES, Qdotloc); CHKERRQ(ierr); 603 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 604 605 // Ceed Vectors 606 ierr = VecGetArrayRead(Qloc, &q); CHKERRQ(ierr); 607 ierr = VecGetArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 608 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 609 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, 610 (PetscScalar *)q); 611 CeedVectorSetArray(user->qdotceed, CEED_MEM_HOST, CEED_USE_POINTER, 612 (PetscScalar *)qdot); 613 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 614 615 // Apply CEED operator 616 CeedOperatorApply(user->op_ifunction, user->qceed, user->gceed, 617 CEED_REQUEST_IMMEDIATE); 618 619 // Restore vectors 620 ierr = VecRestoreArrayRead(Qloc, &q); CHKERRQ(ierr); 621 ierr = VecRestoreArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 622 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 623 624 ierr = VecZeroEntries(G); CHKERRQ(ierr); 625 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 626 627 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 628 ierr = DMRestoreLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 629 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 630 PetscFunctionReturn(0); 631 } 632 633 // User provided TS Monitor 634 static PetscErrorCode TSMonitor_NS(TS ts, PetscInt stepno, PetscReal time, 635 Vec Q, void *ctx) { 636 User user = ctx; 637 Vec Qloc; 638 char filepath[PETSC_MAX_PATH_LEN]; 639 PetscViewer viewer; 640 PetscErrorCode ierr; 641 642 // Set up output 643 PetscFunctionBeginUser; 644 // Print every 'outputfreq' steps 645 if (stepno % user->outputfreq != 0) 646 PetscFunctionReturn(0); 647 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 648 ierr = PetscObjectSetName((PetscObject)Qloc, "StateVec"); CHKERRQ(ierr); 649 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 650 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 651 652 // Output 653 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-%03D.vtu", 654 user->outputfolder, stepno + user->contsteps); 655 CHKERRQ(ierr); 656 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), filepath, 657 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 658 ierr = VecView(Qloc, viewer); CHKERRQ(ierr); 659 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 660 if (user->dmviz) { 661 Vec Qrefined, Qrefined_loc; 662 char filepath_refined[PETSC_MAX_PATH_LEN]; 663 PetscViewer viewer_refined; 664 665 ierr = DMGetGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 666 ierr = DMGetLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 667 ierr = PetscObjectSetName((PetscObject)Qrefined_loc, "Refined"); 668 CHKERRQ(ierr); 669 ierr = MatInterpolate(user->interpviz, Q, Qrefined); CHKERRQ(ierr); 670 ierr = VecZeroEntries(Qrefined_loc); CHKERRQ(ierr); 671 ierr = DMGlobalToLocal(user->dmviz, Qrefined, INSERT_VALUES, Qrefined_loc); 672 CHKERRQ(ierr); 673 ierr = PetscSNPrintf(filepath_refined, sizeof filepath_refined, 674 "%s/nsrefined-%03D.vtu", 675 user->outputfolder, stepno + user->contsteps); 676 CHKERRQ(ierr); 677 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Qrefined), 678 filepath_refined, 679 FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 680 ierr = VecView(Qrefined_loc, viewer_refined); CHKERRQ(ierr); 681 ierr = DMRestoreLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 682 ierr = DMRestoreGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 683 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 684 } 685 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 686 687 // Save data in a binary file for continuation of simulations 688 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 689 user->outputfolder); CHKERRQ(ierr); 690 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 691 CHKERRQ(ierr); 692 ierr = VecView(Q, viewer); CHKERRQ(ierr); 693 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 694 695 // Save time stamp 696 // Dimensionalize time back 697 time /= user->units->second; 698 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 699 user->outputfolder); CHKERRQ(ierr); 700 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 701 CHKERRQ(ierr); 702 #if PETSC_VERSION_GE(3,13,0) 703 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 704 #else 705 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 706 #endif 707 CHKERRQ(ierr); 708 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 709 710 PetscFunctionReturn(0); 711 } 712 713 static PetscErrorCode ICs_FixMultiplicity(CeedOperator op_ics, 714 CeedVector xcorners, CeedVector q0ceed, DM dm, Vec Qloc, Vec Q, 715 CeedElemRestriction restrictq, SetupContext ctxSetup, CeedScalar time) { 716 PetscErrorCode ierr; 717 CeedVector multlvec; 718 Vec Multiplicity, MultiplicityLoc; 719 720 ctxSetup->time = time; 721 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 722 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 723 CeedOperatorApply(op_ics, xcorners, q0ceed, CEED_REQUEST_IMMEDIATE); 724 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 725 ierr = DMLocalToGlobal(dm, Qloc, ADD_VALUES, Q); CHKERRQ(ierr); 726 727 // Fix multiplicity for output of ICs 728 ierr = DMGetLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 729 CeedElemRestrictionCreateVector(restrictq, &multlvec, NULL); 730 ierr = VectorPlacePetscVec(multlvec, MultiplicityLoc); CHKERRQ(ierr); 731 CeedElemRestrictionGetMultiplicity(restrictq, multlvec); 732 CeedVectorDestroy(&multlvec); 733 ierr = DMGetGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 734 ierr = VecZeroEntries(Multiplicity); CHKERRQ(ierr); 735 ierr = DMLocalToGlobal(dm, MultiplicityLoc, ADD_VALUES, Multiplicity); 736 CHKERRQ(ierr); 737 ierr = VecPointwiseDivide(Q, Q, Multiplicity); CHKERRQ(ierr); 738 ierr = VecPointwiseDivide(Qloc, Qloc, MultiplicityLoc); CHKERRQ(ierr); 739 ierr = DMRestoreLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 740 ierr = DMRestoreGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 741 742 PetscFunctionReturn(0); 743 } 744 745 static PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, 746 CeedElemRestriction restrictq, CeedBasis basisq, 747 CeedElemRestriction restrictqdi, CeedVector qdata, Vec M) { 748 PetscErrorCode ierr; 749 CeedQFunction qf_mass; 750 CeedOperator op_mass; 751 CeedVector mceed; 752 Vec Mloc; 753 CeedInt ncompq, qdatasize; 754 755 PetscFunctionBeginUser; 756 CeedElemRestrictionGetNumComponents(restrictq, &ncompq); 757 CeedElemRestrictionGetNumComponents(restrictqdi, &qdatasize); 758 // Create the Q-function that defines the action of the mass operator 759 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 760 CeedQFunctionAddInput(qf_mass, "q", ncompq, CEED_EVAL_INTERP); 761 CeedQFunctionAddInput(qf_mass, "qdata", qdatasize, CEED_EVAL_NONE); 762 CeedQFunctionAddOutput(qf_mass, "v", ncompq, CEED_EVAL_INTERP); 763 764 // Create the mass operator 765 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 766 CeedOperatorSetField(op_mass, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 767 CeedOperatorSetField(op_mass, "qdata", restrictqdi, 768 CEED_BASIS_COLLOCATED, qdata); 769 CeedOperatorSetField(op_mass, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 770 771 ierr = DMGetLocalVector(dm, &Mloc); CHKERRQ(ierr); 772 ierr = VecZeroEntries(Mloc); CHKERRQ(ierr); 773 CeedElemRestrictionCreateVector(restrictq, &mceed, NULL); 774 ierr = VectorPlacePetscVec(mceed, Mloc); CHKERRQ(ierr); 775 776 { 777 // Compute a lumped mass matrix 778 CeedVector onesvec; 779 CeedElemRestrictionCreateVector(restrictq, &onesvec, NULL); 780 CeedVectorSetValue(onesvec, 1.0); 781 CeedOperatorApply(op_mass, onesvec, mceed, CEED_REQUEST_IMMEDIATE); 782 CeedVectorDestroy(&onesvec); 783 CeedOperatorDestroy(&op_mass); 784 CeedVectorDestroy(&mceed); 785 } 786 CeedQFunctionDestroy(&qf_mass); 787 788 ierr = VecZeroEntries(M); CHKERRQ(ierr); 789 ierr = DMLocalToGlobal(dm, Mloc, ADD_VALUES, M); CHKERRQ(ierr); 790 ierr = DMRestoreLocalVector(dm, &Mloc); CHKERRQ(ierr); 791 792 // Invert diagonally lumped mass vector for RHS function 793 ierr = VecReciprocal(M); CHKERRQ(ierr); 794 PetscFunctionReturn(0); 795 } 796 797 static PetscErrorCode SetUpDM(DM dm, problemData *problem, PetscInt degree, 798 SimpleBC bc, void *ctxSetup) { 799 PetscErrorCode ierr; 800 801 PetscFunctionBeginUser; 802 { 803 // Configure the finite element space and boundary conditions 804 PetscFE fe; 805 PetscInt ncompq = 5; 806 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, ncompq, 807 PETSC_FALSE, degree, PETSC_DECIDE, 808 &fe); CHKERRQ(ierr); 809 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 810 ierr = DMAddField(dm,NULL,(PetscObject)fe); CHKERRQ(ierr); 811 ierr = DMCreateDS(dm); CHKERRQ(ierr); 812 { 813 PetscInt comps[1] = {1}; 814 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 815 1, comps, (void(*)(void))NULL, bc->nslip[0], 816 bc->slips[0], ctxSetup); CHKERRQ(ierr); 817 comps[0] = 2; 818 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 819 1, comps, (void(*)(void))NULL, bc->nslip[1], 820 bc->slips[1], ctxSetup); CHKERRQ(ierr); 821 comps[0] = 3; 822 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 823 1, comps, (void(*)(void))NULL, bc->nslip[2], 824 bc->slips[2], ctxSetup); CHKERRQ(ierr); 825 } 826 if (bc->userbc == PETSC_TRUE) { 827 for (PetscInt c = 0; c < 3; c++) { 828 for (PetscInt s = 0; s < bc->nslip[c]; s++) { 829 for (PetscInt w = 0; w < bc->nwall; w++) { 830 if (bc->slips[c][s] == bc->walls[w]) 831 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, 832 "Boundary condition already set on face %D!\n", bc->walls[w]); 833 834 } 835 } 836 } 837 } 838 // Wall boundary conditions are zero energy density and zero flux for 839 // velocity in advection/advection2d, and zero velocity and zero flux 840 // for mass density and energy density in density_current 841 { 842 if (problem->bc == Exact_Advection || problem->bc == Exact_Advection2d) { 843 PetscInt comps[1] = {4}; 844 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 845 1, comps, (void(*)(void))problem->bc, 846 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 847 } else if (problem->bc == Exact_DC) { 848 PetscInt comps[3] = {1, 2, 3}; 849 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 850 3, comps, (void(*)(void))problem->bc, 851 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 852 } else 853 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, 854 "Undefined boundary conditions for this problem"); 855 } 856 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 857 CHKERRQ(ierr); 858 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 859 } 860 { 861 // Empty name for conserved field (because there is only one field) 862 PetscSection section; 863 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 864 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 865 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 866 CHKERRQ(ierr); 867 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 868 CHKERRQ(ierr); 869 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 870 CHKERRQ(ierr); 871 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 872 CHKERRQ(ierr); 873 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 874 CHKERRQ(ierr); 875 } 876 PetscFunctionReturn(0); 877 } 878 879 int main(int argc, char **argv) { 880 PetscInt ierr; 881 MPI_Comm comm; 882 DM dm, dmcoord, dmviz; 883 Mat interpviz; 884 TS ts; 885 TSAdapt adapt; 886 User user; 887 Units units; 888 char ceedresource[4096] = "/cpu/self"; 889 PetscInt localNelemVol, lnodes, gnodes, steps; 890 const PetscInt ncompq = 5; 891 PetscMPIInt rank; 892 PetscScalar ftime; 893 Vec Q, Qloc, Xloc; 894 Ceed ceed; 895 CeedInt numP, numQ; 896 CeedVector xcorners, qdata, q0ceed; 897 CeedBasis basisx, basisxc, basisq; 898 CeedElemRestriction restrictx, restrictq, restrictqdi; 899 CeedQFunction qf_setupVol, qf_ics, qf_rhsVol, qf_ifunctionVol; 900 CeedOperator op_setupVol, op_ics; 901 CeedScalar Rd; 902 CeedMemType memtyperequested; 903 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 904 kgpersquaredms, Joulepercubicm; 905 problemType problemChoice; 906 problemData *problem = NULL; 907 StabilizationType stab; 908 testType testChoice; 909 testData *test = NULL; 910 PetscBool implicit; 911 PetscInt viz_refine = 0; 912 struct SimpleBC_ bc = { 913 .nslip = {2, 2, 2}, 914 .slips = {{5, 6}, {3, 4}, {1, 2}} 915 }; 916 double start, cpu_time_used; 917 // Check PETSc CUDA support 918 PetscBool petschavecuda, setmemtyperequest = PETSC_FALSE; 919 // *INDENT-OFF* 920 #ifdef PETSC_HAVE_CUDA 921 petschavecuda = PETSC_TRUE; 922 #else 923 petschavecuda = PETSC_FALSE; 924 #endif 925 // *INDENT-ON* 926 927 // Create the libCEED contexts 928 PetscScalar meter = 1e-2; // 1 meter in scaled length units 929 PetscScalar second = 1e-2; // 1 second in scaled time units 930 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 931 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 932 CeedScalar theta0 = 300.; // K 933 CeedScalar thetaC = -15.; // K 934 CeedScalar P0 = 1.e5; // Pa 935 CeedScalar P_left = 1.5e5; // Pa 936 CeedScalar rho_left = 1.2; // Kg/m^3 937 CeedScalar N = 0.01; // 1/s 938 CeedScalar cv = 717.; // J/(kg K) 939 CeedScalar cp = 1004.; // J/(kg K) 940 CeedScalar g = 9.81; // m/s^2 941 CeedScalar lambda = -2./3.; // - 942 CeedScalar mu = 75.; // Pa s, dynamic viscosity 943 // mu = 75 is not physical for air, but is good for numerical stability 944 CeedScalar k = 0.02638; // W/(m K) 945 CeedScalar CtauS = 0.; // dimensionless 946 CeedScalar strong_form = 0.; // [0,1] 947 PetscScalar lx = 8000.; // m 948 PetscScalar ly = 8000.; // m 949 PetscScalar lz = 4000.; // m 950 CeedScalar rc = 1000.; // m (Radius of bubble) 951 PetscScalar resx = 1000.; // m (resolution in x) 952 PetscScalar resy = 1000.; // m (resolution in y) 953 PetscScalar resz = 1000.; // m (resolution in z) 954 PetscInt outputfreq = 10; // - 955 PetscInt contsteps = 0; // - 956 PetscInt degree = 1; // - 957 PetscInt qextra = 2; // - 958 PetscInt qextraSur = 2; // - 959 PetscReal center[3], dc_axis[3] = {0, 0, 0}; 960 961 ierr = PetscInitialize(&argc, &argv, NULL, help); 962 if (ierr) return ierr; 963 964 // Allocate PETSc context 965 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 966 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 967 968 // Parse command line options 969 comm = PETSC_COMM_WORLD; 970 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 971 NULL); CHKERRQ(ierr); 972 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 973 NULL, ceedresource, ceedresource, 974 sizeof(ceedresource), NULL); CHKERRQ(ierr); 975 testChoice = TEST_NONE; 976 ierr = PetscOptionsEnum("-test", "Run tests", NULL, 977 testTypes, (PetscEnum)testChoice, 978 (PetscEnum *)&testChoice, 979 NULL); CHKERRQ(ierr); 980 test = &testOptions[testChoice]; 981 problemChoice = NS_DENSITY_CURRENT; 982 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 983 problemTypes, (PetscEnum)problemChoice, 984 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 985 problem = &problemOptions[problemChoice]; 986 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 987 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 988 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 989 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 990 NULL, implicit=PETSC_FALSE, &implicit, NULL); 991 CHKERRQ(ierr); 992 if (!implicit && stab != STAB_NONE) { 993 ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); 994 CHKERRQ(ierr); 995 } 996 { 997 PetscInt len; 998 PetscBool flg; 999 ierr = PetscOptionsIntArray("-bc_outflow", 1000 "Use outflow boundary conditions on this list of faces", 1001 NULL, bc.outflow, 1002 (len = sizeof(bc.outflow) / sizeof(bc.outflow[0]), 1003 &len), &flg); CHKERRQ(ierr); 1004 if (flg) { 1005 bc.noutflow = len; 1006 // Using outflow boundaries disables automatic wall/slip boundaries (they must be set explicitly) 1007 bc.nwall = 0; 1008 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1009 } 1010 ierr = PetscOptionsIntArray("-bc_inflow", 1011 "Use inflow boundary conditions on this list of faces", 1012 NULL, bc.inflow, 1013 (len = sizeof(bc.inflow) / sizeof(bc.inflow[0]), 1014 &len), &flg); CHKERRQ(ierr); 1015 if (flg) { 1016 bc.ninflow = len; 1017 // Using inflow boundaries disables automatic wall/slip boundaries (they must be set explicitly) 1018 bc.nwall = 0; 1019 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1020 } 1021 ierr = PetscOptionsIntArray("-bc_wall", 1022 "Use wall boundary conditions on this list of faces", 1023 NULL, bc.walls, 1024 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 1025 &len), &flg); CHKERRQ(ierr); 1026 if (flg) { 1027 bc.nwall = len; 1028 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 1029 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1030 } 1031 for (PetscInt j=0; j<3; j++) { 1032 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 1033 ierr = PetscOptionsIntArray(flags[j], 1034 "Use slip boundary conditions on this list of faces", 1035 NULL, bc.slips[j], 1036 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 1037 &len), &flg); 1038 CHKERRQ(ierr); 1039 if (flg) { 1040 bc.nslip[j] = len; 1041 bc.userbc = PETSC_TRUE; 1042 } 1043 } 1044 } 1045 ierr = PetscOptionsInt("-viz_refine", 1046 "Regular refinement levels for visualization", 1047 NULL, viz_refine, &viz_refine, NULL); 1048 CHKERRQ(ierr); 1049 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1050 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1051 meter = fabs(meter); 1052 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1053 NULL, second, &second, NULL); CHKERRQ(ierr); 1054 second = fabs(second); 1055 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1056 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1057 kilogram = fabs(kilogram); 1058 ierr = PetscOptionsScalar("-units_Kelvin", 1059 "1 Kelvin in scaled temperature units", 1060 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1061 Kelvin = fabs(Kelvin); 1062 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1063 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1064 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1065 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1066 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1067 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1068 ierr = PetscOptionsScalar("-P_left", "Inflow pressure", 1069 NULL, P_left, &P_left, NULL); CHKERRQ(ierr); 1070 ierr = PetscOptionsScalar("-rho_left", "Inflow density", 1071 NULL, rho_left, &rho_left, NULL); CHKERRQ(ierr); 1072 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1073 NULL, N, &N, NULL); CHKERRQ(ierr); 1074 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1075 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1076 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1077 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1078 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1079 NULL, g, &g, NULL); CHKERRQ(ierr); 1080 ierr = PetscOptionsScalar("-lambda", 1081 "Stokes hypothesis second viscosity coefficient", 1082 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1083 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1084 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1085 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1086 NULL, k, &k, NULL); CHKERRQ(ierr); 1087 ierr = PetscOptionsScalar("-CtauS", 1088 "Scale coefficient for tau (nondimensional)", 1089 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1090 if (stab == STAB_NONE && CtauS != 0) { 1091 ierr = PetscPrintf(comm, 1092 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1093 CHKERRQ(ierr); 1094 } 1095 ierr = PetscOptionsScalar("-strong_form", 1096 "Strong (1) or weak/integrated by parts (0) advection residual", 1097 NULL, strong_form, &strong_form, NULL); 1098 CHKERRQ(ierr); 1099 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1100 ierr = PetscPrintf(comm, 1101 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1102 CHKERRQ(ierr); 1103 } 1104 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1105 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1106 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1107 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1108 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1109 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1110 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1111 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1112 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1113 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1114 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1115 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1116 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1117 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1118 PetscInt n = problem->dim; 1119 center[0] = 0.5 * lx; 1120 center[1] = 0.5 * ly; 1121 center[2] = 0.5 * lz; 1122 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1123 NULL, center, &n, NULL); CHKERRQ(ierr); 1124 n = problem->dim; 1125 ierr = PetscOptionsRealArray("-dc_axis", 1126 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1127 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1128 { 1129 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1130 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1131 if (norm > 0) { 1132 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1133 } 1134 } 1135 ierr = PetscOptionsInt("-output_freq", 1136 "Frequency of output, in number of steps", 1137 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1138 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1139 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1140 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1141 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1142 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1143 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1144 ierr = PetscStrncpy(user->outputfolder, ".", 2); CHKERRQ(ierr); 1145 ierr = PetscOptionsString("-of", "Output folder", 1146 NULL, user->outputfolder, user->outputfolder, 1147 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 1148 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1149 ierr = PetscOptionsEnum("-memtype", 1150 "CEED MemType requested", NULL, 1151 memTypes, (PetscEnum)memtyperequested, 1152 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1153 CHKERRQ(ierr); 1154 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1155 1156 // Define derived units 1157 Pascal = kilogram / (meter * PetscSqr(second)); 1158 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1159 mpersquareds = meter / PetscSqr(second); 1160 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1161 kgpercubicm = kilogram / pow(meter,3); 1162 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1163 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1164 1165 // Scale variables to desired units 1166 theta0 *= Kelvin; 1167 thetaC *= Kelvin; 1168 P0 *= Pascal; 1169 P_left *= Pascal; 1170 rho_left *= kgpercubicm; 1171 N *= (1./second); 1172 cv *= JperkgK; 1173 cp *= JperkgK; 1174 Rd = cp - cv; 1175 g *= mpersquareds; 1176 mu *= Pascal * second; 1177 k *= WpermK; 1178 lx = fabs(lx) * meter; 1179 ly = fabs(ly) * meter; 1180 lz = fabs(lz) * meter; 1181 rc = fabs(rc) * meter; 1182 resx = fabs(resx) * meter; 1183 resy = fabs(resy) * meter; 1184 resz = fabs(resz) * meter; 1185 for (int i=0; i<3; i++) center[i] *= meter; 1186 1187 const CeedInt dim = problem->dim, ncompx = problem->dim, 1188 qdatasizeVol = problem->qdatasizeVol; 1189 // Set up the libCEED context 1190 struct SetupContext_ ctxSetup = { 1191 .theta0 = theta0, 1192 .thetaC = thetaC, 1193 .P0 = P0, 1194 .N = N, 1195 .cv = cv, 1196 .cp = cp, 1197 .Rd = Rd, 1198 .g = g, 1199 .rc = rc, 1200 .lx = lx, 1201 .ly = ly, 1202 .lz = lz, 1203 .center[0] = center[0], 1204 .center[1] = center[1], 1205 .center[2] = center[2], 1206 .dc_axis[0] = dc_axis[0], 1207 .dc_axis[1] = dc_axis[1], 1208 .dc_axis[2] = dc_axis[2], 1209 .time = 0, 1210 }; 1211 1212 // Create the mesh 1213 { 1214 const PetscReal scale[3] = {lx, ly, lz}; 1215 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1216 NULL, PETSC_TRUE, &dm); 1217 CHKERRQ(ierr); 1218 } 1219 1220 // Distribute the mesh over processes 1221 { 1222 DM dmDist = NULL; 1223 PetscPartitioner part; 1224 1225 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1226 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1227 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1228 if (dmDist) { 1229 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1230 dm = dmDist; 1231 } 1232 } 1233 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1234 1235 // Setup DM 1236 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1237 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1238 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetup); CHKERRQ(ierr); 1239 1240 // Refine DM for high-order viz 1241 dmviz = NULL; 1242 interpviz = NULL; 1243 if (viz_refine) { 1244 DM dmhierarchy[viz_refine+1]; 1245 1246 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1247 dmhierarchy[0] = dm; 1248 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1249 Mat interp_next; 1250 1251 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1252 CHKERRQ(ierr); 1253 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1254 d = (d + 1) / 2; 1255 if (i + 1 == viz_refine) d = 1; 1256 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 1257 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1258 &interp_next, NULL); CHKERRQ(ierr); 1259 if (!i) interpviz = interp_next; 1260 else { 1261 Mat C; 1262 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1263 PETSC_DECIDE, &C); CHKERRQ(ierr); 1264 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1265 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1266 interpviz = C; 1267 } 1268 } 1269 for (PetscInt i=1; i<viz_refine; i++) { 1270 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1271 } 1272 dmviz = dmhierarchy[viz_refine]; 1273 } 1274 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1275 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1276 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1277 lnodes /= ncompq; 1278 1279 // Initialize CEED 1280 CeedInit(ceedresource, &ceed); 1281 // Set memtype 1282 CeedMemType memtypebackend; 1283 CeedGetPreferredMemType(ceed, &memtypebackend); 1284 // Check memtype compatibility 1285 if (!setmemtyperequest) 1286 memtyperequested = memtypebackend; 1287 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1288 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1289 "PETSc was not built with CUDA. " 1290 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1291 1292 // Set number of 1D nodes and quadrature points 1293 numP = degree + 1; 1294 numQ = numP + qextra; 1295 1296 // Print summary 1297 if (testChoice == TEST_NONE) { 1298 CeedInt gdofs, odofs; 1299 int comm_size; 1300 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1301 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1302 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1303 gnodes = gdofs/ncompq; 1304 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1305 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1306 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1307 const char *usedresource; 1308 CeedGetResource(ceed, &usedresource); 1309 1310 ierr = PetscPrintf(comm, 1311 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1312 " rank(s) : %d\n" 1313 " Problem:\n" 1314 " Problem Name : %s\n" 1315 " Stabilization : %s\n" 1316 " PETSc:\n" 1317 " Box Faces : %s\n" 1318 " libCEED:\n" 1319 " libCEED Backend : %s\n" 1320 " libCEED Backend MemType : %s\n" 1321 " libCEED User Requested MemType : %s\n" 1322 " Mesh:\n" 1323 " Number of 1D Basis Nodes (P) : %d\n" 1324 " Number of 1D Quadrature Points (Q) : %d\n" 1325 " Global DoFs : %D\n" 1326 " Owned DoFs : %D\n" 1327 " DoFs per node : %D\n" 1328 " Global nodes : %D\n" 1329 " Owned nodes : %D\n", 1330 comm_size, problemTypes[problemChoice], 1331 StabilizationTypes[stab], box_faces_str, usedresource, 1332 CeedMemTypes[memtypebackend], 1333 (setmemtyperequest) ? 1334 CeedMemTypes[memtyperequested] : "none", 1335 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1336 CHKERRQ(ierr); 1337 } 1338 1339 // Set up global mass vector 1340 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1341 1342 // Set up libCEED 1343 // CEED Bases 1344 CeedInit(ceedresource, &ceed); 1345 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1346 &basisq); 1347 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1348 &basisx); 1349 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1350 CEED_GAUSS_LOBATTO, &basisxc); 1351 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1352 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1353 CHKERRQ(ierr); 1354 1355 // CEED Restrictions 1356 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1357 qdatasizeVol, &restrictq, &restrictx, 1358 &restrictqdi); CHKERRQ(ierr); 1359 1360 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1361 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1362 1363 // Create the CEED vectors that will be needed in setup 1364 CeedInt NqptsVol; 1365 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1366 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1367 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1368 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1369 1370 // Create the Q-function that builds the quadrature data for the NS operator 1371 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1372 &qf_setupVol); 1373 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1374 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1375 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1376 1377 // Create the Q-function that sets the ICs of the operator 1378 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1379 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1380 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1381 1382 qf_rhsVol = NULL; 1383 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1384 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1385 problem->applyVol_rhs_loc, &qf_rhsVol); 1386 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1387 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1388 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1389 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1390 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1391 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1392 } 1393 1394 qf_ifunctionVol = NULL; 1395 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1396 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1397 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1398 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1399 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1400 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1401 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1402 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1403 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1404 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1405 } 1406 1407 // Create the operator that builds the quadrature data for the NS operator 1408 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1409 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1410 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1411 basisx, CEED_VECTOR_NONE); 1412 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1413 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1414 1415 // Create the operator that sets the ICs 1416 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1417 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1418 CeedOperatorSetField(op_ics, "q0", restrictq, 1419 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1420 1421 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1422 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1423 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1424 1425 if (qf_rhsVol) { // Create the RHS physics operator 1426 CeedOperator op; 1427 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1428 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1429 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1430 CeedOperatorSetField(op, "qdata", restrictqdi, 1431 CEED_BASIS_COLLOCATED, qdata); 1432 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1433 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1434 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1435 user->op_rhs_vol = op; 1436 } 1437 1438 if (qf_ifunctionVol) { // Create the IFunction operator 1439 CeedOperator op; 1440 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1441 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1442 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1443 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1444 CeedOperatorSetField(op, "qdata", restrictqdi, 1445 CEED_BASIS_COLLOCATED, qdata); 1446 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1447 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1448 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1449 user->op_ifunction_vol = op; 1450 } 1451 1452 //--------------------------------------------------------------------------------------// 1453 // In/OutFlow Boundary Conditions 1454 //--------------------------------------------------------------------------------------// 1455 // Set up CEED for the boundaries 1456 CeedInt height = 1; 1457 CeedInt dimSur = dim - height; 1458 CeedInt numP_Sur = degree + 1; 1459 CeedInt numQ_Sur = numP_Sur + qextraSur; 1460 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1461 CeedBasis basisxSur, basisxcSur, basisqSur; 1462 CeedInt NqptsSur; 1463 CeedQFunction qf_setupSur, qf_rhsOut, qf_ifunctionOut, qf_rhsIn, qf_ifunctionIn; 1464 1465 // CEED bases for the boundaries 1466 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, CEED_GAUSS, 1467 &basisqSur); 1468 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1469 &basisxSur); 1470 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1471 CEED_GAUSS_LOBATTO, &basisxcSur); 1472 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1473 1474 // Create the Q-function that builds the quadrature data for the Surface operator 1475 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1476 &qf_setupSur); 1477 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1478 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1479 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1480 1481 // Creat Q-Function for OutFlow BCs 1482 qf_rhsOut = NULL; 1483 if (problem->applyOut_rhs) { // Create the Q-function that defines the action of the RHS operator 1484 CeedQFunctionCreateInterior(ceed, 1, problem->applyOut_rhs, 1485 problem->applyOut_rhs_loc, &qf_rhsOut); 1486 CeedQFunctionAddInput(qf_rhsOut, "q", ncompq, CEED_EVAL_INTERP); 1487 CeedQFunctionAddInput(qf_rhsOut, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1488 CeedQFunctionAddInput(qf_rhsOut, "x", ncompx, CEED_EVAL_INTERP); 1489 CeedQFunctionAddOutput(qf_rhsOut, "v", ncompq, CEED_EVAL_INTERP); 1490 } 1491 qf_ifunctionOut = NULL; 1492 if (problem->applyOut_ifunction) { // Create the Q-function that defines the action of the IFunction 1493 CeedQFunctionCreateInterior(ceed, 1, problem->applyOut_ifunction, 1494 problem->applyOut_ifunction_loc, &qf_ifunctionOut); 1495 CeedQFunctionAddInput(qf_ifunctionOut, "q", ncompq, CEED_EVAL_INTERP); 1496 CeedQFunctionAddInput(qf_ifunctionOut, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1497 CeedQFunctionAddInput(qf_ifunctionOut, "x", ncompx, CEED_EVAL_INTERP); 1498 CeedQFunctionAddOutput(qf_ifunctionOut, "v", ncompq, CEED_EVAL_INTERP); 1499 } 1500 1501 // Creat Q-Function for InFlow BCs 1502 qf_rhsIn = NULL; 1503 if (problem->applyIn_rhs) { // Create the Q-function that defines the action of the RHS operator 1504 CeedQFunctionCreateInterior(ceed, 1, problem->applyIn_rhs, 1505 problem->applyIn_rhs_loc, &qf_rhsIn); 1506 CeedQFunctionAddInput(qf_rhsIn, "q", ncompq, CEED_EVAL_INTERP); 1507 CeedQFunctionAddInput(qf_rhsIn, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1508 CeedQFunctionAddInput(qf_rhsIn, "x", ncompx, CEED_EVAL_INTERP); 1509 CeedQFunctionAddOutput(qf_rhsIn, "v", ncompq, CEED_EVAL_INTERP); 1510 } 1511 qf_ifunctionIn = NULL; 1512 if (problem->applyIn_ifunction) { // Create the Q-function that defines the action of the IFunction 1513 CeedQFunctionCreateInterior(ceed, 1, problem->applyIn_ifunction, 1514 problem->applyIn_ifunction_loc, &qf_ifunctionIn); 1515 CeedQFunctionAddInput(qf_ifunctionIn, "q", ncompq, CEED_EVAL_INTERP); 1516 CeedQFunctionAddInput(qf_ifunctionIn, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1517 CeedQFunctionAddInput(qf_ifunctionIn, "x", ncompx, CEED_EVAL_INTERP); 1518 CeedQFunctionAddOutput(qf_ifunctionIn, "v", ncompq, CEED_EVAL_INTERP); 1519 } 1520 1521 // Create CEED Operator for the whole domain 1522 if (!implicit) 1523 ierr = CreateOperatorForDomain(ceed, dm, user->op_rhs_vol, qf_rhsOut, qf_rhsIn, qf_setupSur, height, 1524 bc.noutflow, bc.outflow, bc.ninflow, bc.inflow, numP_Sur, numQ_Sur, qdatasizeSur, 1525 NqptsSur, basisxSur, basisqSur, &user->op_rhs); 1526 CHKERRQ(ierr); 1527 if (implicit) 1528 ierr = CreateOperatorForDomain(ceed, dm, user->op_ifunction_vol, qf_ifunctionOut, qf_ifunctionIn, qf_setupSur, height, 1529 bc.noutflow, bc.outflow, bc.ninflow, bc.inflow, numP_Sur, numQ_Sur, qdatasizeSur, 1530 NqptsSur, basisxSur, basisqSur, &user->op_ifunction); 1531 CHKERRQ(ierr); 1532 1533 //--------------------------------------------------------------------------------------// 1534 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1535 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1536 CeedScalar ctxIn[5] = {cv, cp, Rd, P_left, rho_left}; 1537 struct Advection2dContext_ ctxAdvection2d = { 1538 .CtauS = CtauS, 1539 .strong_form = strong_form, 1540 .stabilization = stab, 1541 }; 1542 switch (problemChoice) { 1543 case NS_DENSITY_CURRENT: 1544 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1545 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1546 sizeof ctxNS); 1547 if (qf_rhsOut) CeedQFunctionSetContext(qf_rhsOut, &ctxNS, sizeof ctxNS); 1548 if (qf_ifunctionOut) CeedQFunctionSetContext(qf_ifunctionOut, &ctxNS, 1549 sizeof ctxNS); 1550 if (qf_rhsIn) CeedQFunctionSetContext(qf_rhsIn, &ctxIn, sizeof ctxIn); 1551 if (qf_ifunctionIn) CeedQFunctionSetContext(qf_ifunctionIn, &ctxIn, sizeof ctxIn); 1552 break; 1553 case NS_ADVECTION: 1554 case NS_ADVECTION2D: 1555 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1556 sizeof ctxAdvection2d); 1557 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1558 sizeof ctxAdvection2d); 1559 if (qf_rhsOut) CeedQFunctionSetContext(qf_rhsOut, &ctxAdvection2d, 1560 sizeof ctxAdvection2d); 1561 if (qf_ifunctionOut) CeedQFunctionSetContext(qf_ifunctionOut, &ctxAdvection2d, 1562 sizeof ctxAdvection2d); 1563 if (qf_rhsIn) CeedQFunctionSetContext(qf_rhsIn, &ctxIn, sizeof ctxIn); 1564 if (qf_ifunctionIn) CeedQFunctionSetContext(qf_ifunctionIn, &ctxIn, sizeof ctxIn); 1565 } 1566 1567 // Set up PETSc context 1568 // Set up units structure 1569 units->meter = meter; 1570 units->kilogram = kilogram; 1571 units->second = second; 1572 units->Kelvin = Kelvin; 1573 units->Pascal = Pascal; 1574 units->JperkgK = JperkgK; 1575 units->mpersquareds = mpersquareds; 1576 units->WpermK = WpermK; 1577 units->kgpercubicm = kgpercubicm; 1578 units->kgpersquaredms = kgpersquaredms; 1579 units->Joulepercubicm = Joulepercubicm; 1580 1581 // Set up user structure 1582 user->comm = comm; 1583 user->outputfreq = outputfreq; 1584 user->contsteps = contsteps; 1585 user->units = units; 1586 user->dm = dm; 1587 user->dmviz = dmviz; 1588 user->interpviz = interpviz; 1589 user->ceed = ceed; 1590 1591 // Calculate qdata and ICs 1592 // Set up state global and local vectors 1593 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1594 1595 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1596 1597 // Apply Setup Ceed Operators 1598 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1599 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1600 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1601 user->M); CHKERRQ(ierr); 1602 1603 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1604 &ctxSetup, 0.0); CHKERRQ(ierr); 1605 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1606 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1607 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1608 // slower execution. 1609 Vec Qbc; 1610 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1611 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1612 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1613 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1614 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1615 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1616 ierr = PetscObjectComposeFunction((PetscObject)dm, 1617 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1618 CHKERRQ(ierr); 1619 } 1620 1621 MPI_Comm_rank(comm, &rank); 1622 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1623 // Gather initial Q values 1624 // In case of continuation of simulation, set up initial values from binary file 1625 if (contsteps) { // continue from existent solution 1626 PetscViewer viewer; 1627 char filepath[PETSC_MAX_PATH_LEN]; 1628 // Read input 1629 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1630 user->outputfolder); 1631 CHKERRQ(ierr); 1632 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1633 CHKERRQ(ierr); 1634 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1635 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1636 } 1637 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1638 1639 // Create and setup TS 1640 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1641 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1642 if (implicit) { 1643 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1644 if (user->op_ifunction) { 1645 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1646 } else { // Implicit integrators can fall back to using an RHSFunction 1647 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1648 } 1649 } else { 1650 if (!user->op_rhs) SETERRQ(comm,PETSC_ERR_ARG_NULL, 1651 "Problem does not provide RHSFunction"); 1652 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1653 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1654 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1655 } 1656 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1657 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1658 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1659 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1660 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1661 ierr = TSAdaptSetStepLimits(adapt, 1662 1.e-12 * units->second, 1663 1.e2 * units->second); CHKERRQ(ierr); 1664 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1665 if (!contsteps) { // print initial condition 1666 if (testChoice == TEST_NONE) { 1667 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1668 } 1669 } else { // continue from time of last output 1670 PetscReal time; 1671 PetscInt count; 1672 PetscViewer viewer; 1673 char filepath[PETSC_MAX_PATH_LEN]; 1674 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1675 user->outputfolder); CHKERRQ(ierr); 1676 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1677 CHKERRQ(ierr); 1678 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1679 CHKERRQ(ierr); 1680 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1681 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1682 } 1683 if (testChoice == TEST_NONE) { 1684 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1685 } 1686 1687 // Solve 1688 start = MPI_Wtime(); 1689 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1690 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1691 cpu_time_used = MPI_Wtime() - start; 1692 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1693 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1694 comm); CHKERRQ(ierr); 1695 if (testChoice == TEST_NONE) { 1696 ierr = PetscPrintf(PETSC_COMM_WORLD, 1697 "Time taken for solution (sec): %g\n", 1698 (double)cpu_time_used); CHKERRQ(ierr); 1699 } 1700 1701 // Get error 1702 if (problem->non_zero_time && testChoice == TEST_NONE) { 1703 Vec Qexact, Qexactloc; 1704 PetscReal norm; 1705 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1706 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1707 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1708 1709 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1710 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1711 1712 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1713 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1714 CeedVectorDestroy(&q0ceed); 1715 ierr = PetscPrintf(PETSC_COMM_WORLD, 1716 "Max Error: %g\n", 1717 (double)norm); CHKERRQ(ierr); 1718 // Clean up vectors 1719 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1720 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1721 } 1722 1723 // Output Statistics 1724 ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr); 1725 if (testChoice == TEST_NONE) { 1726 ierr = PetscPrintf(PETSC_COMM_WORLD, 1727 "Time integrator took %D time steps to reach final time %g\n", 1728 steps, (double)ftime); CHKERRQ(ierr); 1729 } 1730 // Output numerical values from command line 1731 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1732 1733 // compare reference solution values with current run 1734 if (testChoice != TEST_NONE) { 1735 PetscViewer viewer; 1736 // Read reference file 1737 Vec Qref; 1738 PetscReal error, Qrefnorm; 1739 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1740 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1741 CHKERRQ(ierr); 1742 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1743 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1744 1745 // Compute error with respect to reference solution 1746 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1747 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1748 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1749 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1750 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1751 // Check error 1752 if (error > test->testtol) { 1753 ierr = PetscPrintf(PETSC_COMM_WORLD, 1754 "Test failed with error norm %g\n", 1755 (double)error); CHKERRQ(ierr); 1756 } 1757 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1758 } 1759 1760 // Clean up libCEED 1761 CeedVectorDestroy(&qdata); 1762 CeedVectorDestroy(&user->qceed); 1763 CeedVectorDestroy(&user->qdotceed); 1764 CeedVectorDestroy(&user->gceed); 1765 CeedVectorDestroy(&xcorners); 1766 CeedBasisDestroy(&basisq); 1767 CeedBasisDestroy(&basisx); 1768 CeedBasisDestroy(&basisxc); 1769 CeedElemRestrictionDestroy(&restrictq); 1770 CeedElemRestrictionDestroy(&restrictx); 1771 CeedElemRestrictionDestroy(&restrictqdi); 1772 CeedQFunctionDestroy(&qf_setupVol); 1773 CeedQFunctionDestroy(&qf_ics); 1774 CeedQFunctionDestroy(&qf_rhsVol); 1775 CeedQFunctionDestroy(&qf_ifunctionVol); 1776 CeedOperatorDestroy(&op_setupVol); 1777 CeedOperatorDestroy(&op_ics); 1778 CeedOperatorDestroy(&user->op_rhs_vol); 1779 CeedOperatorDestroy(&user->op_ifunction_vol); 1780 CeedDestroy(&ceed); 1781 CeedBasisDestroy(&basisqSur); 1782 CeedBasisDestroy(&basisxSur); 1783 CeedBasisDestroy(&basisxcSur); 1784 CeedQFunctionDestroy(&qf_setupSur); 1785 CeedQFunctionDestroy(&qf_rhsOut); 1786 CeedQFunctionDestroy(&qf_ifunctionOut); 1787 CeedQFunctionDestroy(&qf_rhsIn); 1788 CeedQFunctionDestroy(&qf_ifunctionIn); 1789 CeedOperatorDestroy(&user->op_rhs); 1790 CeedOperatorDestroy(&user->op_ifunction); 1791 1792 // Clean up PETSc 1793 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1794 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1795 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1796 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1797 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1798 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1799 ierr = PetscFree(units); CHKERRQ(ierr); 1800 ierr = PetscFree(user); CHKERRQ(ierr); 1801 return PetscFinalize(); 1802 } 1803 1804