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