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_left = 1.5e5; // Pa 938 CeedScalar rho_left = 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 if (wind_type == ADVECTION_WIND_TRANSLATION) { 992 PetscInt n = problem->dim; 993 ierr = PetscOptionsRealArray("-problem_advection_wind_translation", "Constant wind vector", 994 NULL, wind, &n, NULL); CHKERRQ(ierr); 995 } 996 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 997 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 998 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 999 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 1000 NULL, implicit=PETSC_FALSE, &implicit, NULL); 1001 CHKERRQ(ierr); 1002 if (!implicit && stab != STAB_NONE) { 1003 ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); 1004 CHKERRQ(ierr); 1005 } 1006 { 1007 PetscInt len; 1008 PetscBool flg; 1009 ierr = PetscOptionsIntArray("-bc_outflow", 1010 "Use outflow boundary conditions on this list of faces", 1011 NULL, bc.outflow, 1012 (len = sizeof(bc.outflow) / sizeof(bc.outflow[0]), 1013 &len), &flg); CHKERRQ(ierr); 1014 if (flg) { 1015 bc.noutflow = len; 1016 // Using outflow boundaries disables automatic wall/slip boundaries (they must be set explicitly) 1017 bc.nwall = 0; 1018 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1019 } 1020 ierr = PetscOptionsIntArray("-bc_inflow", 1021 "Use inflow boundary conditions on this list of faces", 1022 NULL, bc.inflow, 1023 (len = sizeof(bc.inflow) / sizeof(bc.inflow[0]), 1024 &len), &flg); CHKERRQ(ierr); 1025 if (flg) { 1026 bc.ninflow = len; 1027 // Using inflow boundaries disables automatic wall/slip boundaries (they must be set explicitly) 1028 bc.nwall = 0; 1029 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1030 } 1031 ierr = PetscOptionsIntArray("-bc_wall", 1032 "Use wall boundary conditions on this list of faces", 1033 NULL, bc.walls, 1034 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 1035 &len), &flg); CHKERRQ(ierr); 1036 if (flg) { 1037 bc.nwall = len; 1038 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 1039 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1040 } 1041 for (PetscInt j=0; j<3; j++) { 1042 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 1043 ierr = PetscOptionsIntArray(flags[j], 1044 "Use slip boundary conditions on this list of faces", 1045 NULL, bc.slips[j], 1046 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 1047 &len), &flg); 1048 CHKERRQ(ierr); 1049 if (flg) { 1050 bc.nslip[j] = len; 1051 bc.userbc = PETSC_TRUE; 1052 } 1053 } 1054 } 1055 ierr = PetscOptionsInt("-viz_refine", 1056 "Regular refinement levels for visualization", 1057 NULL, viz_refine, &viz_refine, NULL); 1058 CHKERRQ(ierr); 1059 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1060 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1061 meter = fabs(meter); 1062 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1063 NULL, second, &second, NULL); CHKERRQ(ierr); 1064 second = fabs(second); 1065 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1066 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1067 kilogram = fabs(kilogram); 1068 ierr = PetscOptionsScalar("-units_Kelvin", 1069 "1 Kelvin in scaled temperature units", 1070 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1071 Kelvin = fabs(Kelvin); 1072 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1073 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1074 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1075 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1076 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1077 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1078 ierr = PetscOptionsScalar("-P_left", "Inflow pressure", 1079 NULL, P_left, &P_left, NULL); CHKERRQ(ierr); 1080 ierr = PetscOptionsScalar("-rho_left", "Inflow density", 1081 NULL, rho_left, &rho_left, NULL); CHKERRQ(ierr); 1082 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1083 NULL, N, &N, NULL); CHKERRQ(ierr); 1084 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1085 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1086 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1087 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1088 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1089 NULL, g, &g, NULL); CHKERRQ(ierr); 1090 ierr = PetscOptionsScalar("-lambda", 1091 "Stokes hypothesis second viscosity coefficient", 1092 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1093 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1094 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1095 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1096 NULL, k, &k, NULL); CHKERRQ(ierr); 1097 ierr = PetscOptionsScalar("-CtauS", 1098 "Scale coefficient for tau (nondimensional)", 1099 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1100 if (stab == STAB_NONE && CtauS != 0) { 1101 ierr = PetscPrintf(comm, 1102 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1103 CHKERRQ(ierr); 1104 } 1105 ierr = PetscOptionsScalar("-strong_form", 1106 "Strong (1) or weak/integrated by parts (0) advection residual", 1107 NULL, strong_form, &strong_form, NULL); 1108 CHKERRQ(ierr); 1109 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1110 ierr = PetscPrintf(comm, 1111 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1112 CHKERRQ(ierr); 1113 } 1114 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1115 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1116 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1117 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1118 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1119 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1120 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1121 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1122 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1123 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1124 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1125 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1126 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1127 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1128 PetscInt n = problem->dim; 1129 center[0] = 0.5 * lx; 1130 center[1] = 0.5 * ly; 1131 center[2] = 0.5 * lz; 1132 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1133 NULL, center, &n, NULL); CHKERRQ(ierr); 1134 n = problem->dim; 1135 ierr = PetscOptionsRealArray("-dc_axis", 1136 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1137 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1138 { 1139 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1140 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1141 if (norm > 0) { 1142 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1143 } 1144 } 1145 ierr = PetscOptionsInt("-output_freq", 1146 "Frequency of output, in number of steps", 1147 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1148 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1149 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1150 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1151 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1152 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1153 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1154 ierr = PetscStrncpy(user->outputfolder, ".", 2); CHKERRQ(ierr); 1155 ierr = PetscOptionsString("-of", "Output folder", 1156 NULL, user->outputfolder, user->outputfolder, 1157 sizeof(user->outputfolder), NULL); CHKERRQ(ierr); 1158 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1159 ierr = PetscOptionsEnum("-memtype", 1160 "CEED MemType requested", NULL, 1161 memTypes, (PetscEnum)memtyperequested, 1162 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1163 CHKERRQ(ierr); 1164 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1165 1166 // Setup BCs for Rotation or Translation wind types in Advection (3d) 1167 if (problemChoice == NS_ADVECTION) { 1168 switch (wind_type) { 1169 case ADVECTION_WIND_ROTATION: 1170 // No in/out-flow 1171 bc.ninflow = bc.noutflow = 0; 1172 break; 1173 case ADVECTION_WIND_TRANSLATION: 1174 // Face 6 is inflow and Face 5 is outflow 1175 bc.ninflow = bc.noutflow = 1; 1176 bc.inflow[0] = 6; bc.outflow[0] = 5; 1177 // Faces 3 and 4 are slip 1178 bc.nslip[0] = bc.nslip[2] = 0; bc.nslip[1] = 2; 1179 bc.slips[1][0] = 3; bc.slips[1][1] = 4; 1180 // Faces 1 and 2 are wall 1181 bc.nwall = 2; 1182 bc.walls[0] = 1; bc.walls[1] = 2; 1183 break; 1184 } 1185 } 1186 // Setup BCs for Rotation or Translation wind types in Advection (2d) 1187 if (problemChoice == NS_ADVECTION2D) { 1188 switch (wind_type) { 1189 case ADVECTION_WIND_ROTATION: 1190 // No in/out-flow 1191 bc.ninflow = bc.noutflow = 0; 1192 break; 1193 case ADVECTION_WIND_TRANSLATION: 1194 // Face 4 is inflow and Face 2 is outflow 1195 bc.ninflow = bc.noutflow = 1; 1196 bc.inflow[0] = 4; bc.outflow[0] = 2; 1197 // Face 3 is slip 1198 bc.nslip[0] = bc.nslip[2] = 0; bc.nslip[1] = 1; 1199 bc.slips[1][0] = 3; 1200 // Face 1 is wall 1201 bc.nwall = 1; 1202 bc.walls[0] = 1; 1203 break; 1204 } 1205 } 1206 1207 // Define derived units 1208 Pascal = kilogram / (meter * PetscSqr(second)); 1209 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1210 mpersquareds = meter / PetscSqr(second); 1211 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1212 kgpercubicm = kilogram / pow(meter,3); 1213 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1214 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1215 1216 // Scale variables to desired units 1217 theta0 *= Kelvin; 1218 thetaC *= Kelvin; 1219 P0 *= Pascal; 1220 P_left *= Pascal; 1221 rho_left *= kgpercubicm; 1222 N *= (1./second); 1223 cv *= JperkgK; 1224 cp *= JperkgK; 1225 Rd = cp - cv; 1226 g *= mpersquareds; 1227 mu *= Pascal * second; 1228 k *= WpermK; 1229 lx = fabs(lx) * meter; 1230 ly = fabs(ly) * meter; 1231 lz = fabs(lz) * meter; 1232 rc = fabs(rc) * meter; 1233 resx = fabs(resx) * meter; 1234 resy = fabs(resy) * meter; 1235 resz = fabs(resz) * meter; 1236 for (int i=0; i<3; i++) center[i] *= meter; 1237 1238 const CeedInt dim = problem->dim, ncompx = problem->dim, 1239 qdatasizeVol = problem->qdatasizeVol; 1240 // Set up the libCEED context 1241 struct SetupContext_ ctxSetup = { 1242 .theta0 = theta0, 1243 .thetaC = thetaC, 1244 .P0 = P0, 1245 .N = N, 1246 .cv = cv, 1247 .cp = cp, 1248 .Rd = Rd, 1249 .g = g, 1250 .rc = rc, 1251 .lx = lx, 1252 .ly = ly, 1253 .lz = lz, 1254 .center[0] = center[0], 1255 .center[1] = center[1], 1256 .center[2] = center[2], 1257 .dc_axis[0] = dc_axis[0], 1258 .dc_axis[1] = dc_axis[1], 1259 .dc_axis[2] = dc_axis[2], 1260 .wind[0] = wind[0], 1261 .wind[1] = wind[1], 1262 .wind[2] = wind[2], 1263 .time = 0, 1264 .wind_type = wind_type, 1265 }; 1266 1267 // Create the mesh 1268 { 1269 const PetscReal scale[3] = {lx, ly, lz}; 1270 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1271 NULL, PETSC_TRUE, &dm); 1272 CHKERRQ(ierr); 1273 } 1274 1275 // Distribute the mesh over processes 1276 { 1277 DM dmDist = NULL; 1278 PetscPartitioner part; 1279 1280 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1281 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1282 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1283 if (dmDist) { 1284 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1285 dm = dmDist; 1286 } 1287 } 1288 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1289 1290 // Setup DM 1291 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1292 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1293 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetup); CHKERRQ(ierr); 1294 1295 // Refine DM for high-order viz 1296 dmviz = NULL; 1297 interpviz = NULL; 1298 if (viz_refine) { 1299 DM dmhierarchy[viz_refine+1]; 1300 1301 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1302 dmhierarchy[0] = dm; 1303 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1304 Mat interp_next; 1305 1306 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1307 CHKERRQ(ierr); 1308 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1309 d = (d + 1) / 2; 1310 if (i + 1 == viz_refine) d = 1; 1311 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetup); CHKERRQ(ierr); 1312 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1313 &interp_next, NULL); CHKERRQ(ierr); 1314 if (!i) interpviz = interp_next; 1315 else { 1316 Mat C; 1317 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1318 PETSC_DECIDE, &C); CHKERRQ(ierr); 1319 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1320 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1321 interpviz = C; 1322 } 1323 } 1324 for (PetscInt i=1; i<viz_refine; i++) { 1325 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1326 } 1327 dmviz = dmhierarchy[viz_refine]; 1328 } 1329 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1330 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1331 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1332 lnodes /= ncompq; 1333 1334 // Initialize CEED 1335 CeedInit(ceedresource, &ceed); 1336 // Set memtype 1337 CeedMemType memtypebackend; 1338 CeedGetPreferredMemType(ceed, &memtypebackend); 1339 // Check memtype compatibility 1340 if (!setmemtyperequest) 1341 memtyperequested = memtypebackend; 1342 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1343 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1344 "PETSc was not built with CUDA. " 1345 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1346 1347 // Set number of 1D nodes and quadrature points 1348 numP = degree + 1; 1349 numQ = numP + qextra; 1350 1351 // Print summary 1352 if (testChoice == TEST_NONE) { 1353 CeedInt gdofs, odofs; 1354 int comm_size; 1355 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1356 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1357 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1358 gnodes = gdofs/ncompq; 1359 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1360 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1361 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1362 const char *usedresource; 1363 CeedGetResource(ceed, &usedresource); 1364 1365 ierr = PetscPrintf(comm, 1366 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1367 " rank(s) : %d\n" 1368 " Problem:\n" 1369 " Problem Name : %s\n" 1370 " Stabilization : %s\n" 1371 " PETSc:\n" 1372 " Box Faces : %s\n" 1373 " libCEED:\n" 1374 " libCEED Backend : %s\n" 1375 " libCEED Backend MemType : %s\n" 1376 " libCEED User Requested MemType : %s\n" 1377 " Mesh:\n" 1378 " Number of 1D Basis Nodes (P) : %d\n" 1379 " Number of 1D Quadrature Points (Q) : %d\n" 1380 " Global DoFs : %D\n" 1381 " Owned DoFs : %D\n" 1382 " DoFs per node : %D\n" 1383 " Global nodes : %D\n" 1384 " Owned nodes : %D\n", 1385 comm_size, problemTypes[problemChoice], 1386 StabilizationTypes[stab], box_faces_str, usedresource, 1387 CeedMemTypes[memtypebackend], 1388 (setmemtyperequest) ? 1389 CeedMemTypes[memtyperequested] : "none", 1390 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1391 CHKERRQ(ierr); 1392 } 1393 1394 // Set up global mass vector 1395 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1396 1397 // Set up libCEED 1398 // CEED Bases 1399 CeedInit(ceedresource, &ceed); 1400 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1401 &basisq); 1402 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1403 &basisx); 1404 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1405 CEED_GAUSS_LOBATTO, &basisxc); 1406 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1407 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1408 CHKERRQ(ierr); 1409 1410 // CEED Restrictions 1411 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1412 qdatasizeVol, &restrictq, &restrictx, 1413 &restrictqdi); CHKERRQ(ierr); 1414 1415 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1416 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1417 1418 // Create the CEED vectors that will be needed in setup 1419 CeedInt NqptsVol; 1420 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1421 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1422 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1423 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1424 1425 // Create the Q-function that builds the quadrature data for the NS operator 1426 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1427 &qf_setupVol); 1428 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1429 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1430 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1431 1432 // Create the Q-function that sets the ICs of the operator 1433 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1434 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1435 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1436 1437 qf_rhsVol = NULL; 1438 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1439 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1440 problem->applyVol_rhs_loc, &qf_rhsVol); 1441 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1442 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1443 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1444 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1445 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1446 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1447 } 1448 1449 qf_ifunctionVol = NULL; 1450 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1451 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1452 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1453 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1454 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1455 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1456 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1457 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1458 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1459 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1460 } 1461 1462 // Create the operator that builds the quadrature data for the NS operator 1463 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1464 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1465 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1466 basisx, CEED_VECTOR_NONE); 1467 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1468 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1469 1470 // Create the operator that sets the ICs 1471 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1472 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1473 CeedOperatorSetField(op_ics, "q0", restrictq, 1474 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1475 1476 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1477 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1478 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1479 1480 if (qf_rhsVol) { // Create the RHS physics operator 1481 CeedOperator op; 1482 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1483 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1484 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1485 CeedOperatorSetField(op, "qdata", restrictqdi, 1486 CEED_BASIS_COLLOCATED, qdata); 1487 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1488 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1489 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1490 user->op_rhs_vol = op; 1491 } 1492 1493 if (qf_ifunctionVol) { // Create the IFunction operator 1494 CeedOperator op; 1495 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1496 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1497 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1498 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1499 CeedOperatorSetField(op, "qdata", restrictqdi, 1500 CEED_BASIS_COLLOCATED, qdata); 1501 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1502 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1503 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1504 user->op_ifunction_vol = op; 1505 } 1506 1507 //--------------------------------------------------------------------------------------// 1508 // In/OutFlow Boundary Conditions 1509 //--------------------------------------------------------------------------------------// 1510 // Set up CEED for the boundaries 1511 CeedInt height = 1; 1512 CeedInt dimSur = dim - height; 1513 CeedInt numP_Sur = degree + 1; 1514 CeedInt numQ_Sur = numP_Sur + qextraSur; 1515 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1516 CeedBasis basisxSur, basisxcSur, basisqSur; 1517 CeedInt NqptsSur; 1518 CeedQFunction qf_setupSur, qf_rhsOut, qf_ifunctionOut, qf_rhsIn, qf_ifunctionIn; 1519 1520 // CEED bases for the boundaries 1521 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, CEED_GAUSS, 1522 &basisqSur); 1523 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1524 &basisxSur); 1525 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1526 CEED_GAUSS_LOBATTO, &basisxcSur); 1527 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1528 1529 // Create the Q-function that builds the quadrature data for the Surface operator 1530 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1531 &qf_setupSur); 1532 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1533 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1534 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1535 1536 // Creat Q-Function for OutFlow BCs 1537 qf_rhsOut = NULL; 1538 if (problem->applyOut_rhs) { // Create the Q-function that defines the action of the RHS operator 1539 CeedQFunctionCreateInterior(ceed, 1, problem->applyOut_rhs, 1540 problem->applyOut_rhs_loc, &qf_rhsOut); 1541 CeedQFunctionAddInput(qf_rhsOut, "q", ncompq, CEED_EVAL_INTERP); 1542 CeedQFunctionAddInput(qf_rhsOut, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1543 CeedQFunctionAddInput(qf_rhsOut, "x", ncompx, CEED_EVAL_INTERP); 1544 CeedQFunctionAddOutput(qf_rhsOut, "v", ncompq, CEED_EVAL_INTERP); 1545 } 1546 qf_ifunctionOut = NULL; 1547 if (problem->applyOut_ifunction) { // Create the Q-function that defines the action of the IFunction 1548 CeedQFunctionCreateInterior(ceed, 1, problem->applyOut_ifunction, 1549 problem->applyOut_ifunction_loc, &qf_ifunctionOut); 1550 CeedQFunctionAddInput(qf_ifunctionOut, "q", ncompq, CEED_EVAL_INTERP); 1551 CeedQFunctionAddInput(qf_ifunctionOut, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1552 CeedQFunctionAddInput(qf_ifunctionOut, "x", ncompx, CEED_EVAL_INTERP); 1553 CeedQFunctionAddOutput(qf_ifunctionOut, "v", ncompq, CEED_EVAL_INTERP); 1554 } 1555 1556 // Creat Q-Function for InFlow BCs 1557 qf_rhsIn = NULL; 1558 if (problem->applyIn_rhs) { // Create the Q-function that defines the action of the RHS operator 1559 CeedQFunctionCreateInterior(ceed, 1, problem->applyIn_rhs, 1560 problem->applyIn_rhs_loc, &qf_rhsIn); 1561 CeedQFunctionAddInput(qf_rhsIn, "q", ncompq, CEED_EVAL_INTERP); 1562 CeedQFunctionAddInput(qf_rhsIn, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1563 CeedQFunctionAddInput(qf_rhsIn, "x", ncompx, CEED_EVAL_INTERP); 1564 CeedQFunctionAddOutput(qf_rhsIn, "v", ncompq, CEED_EVAL_INTERP); 1565 } 1566 qf_ifunctionIn = NULL; 1567 if (problem->applyIn_ifunction) { // Create the Q-function that defines the action of the IFunction 1568 CeedQFunctionCreateInterior(ceed, 1, problem->applyIn_ifunction, 1569 problem->applyIn_ifunction_loc, &qf_ifunctionIn); 1570 CeedQFunctionAddInput(qf_ifunctionIn, "q", ncompq, CEED_EVAL_INTERP); 1571 CeedQFunctionAddInput(qf_ifunctionIn, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1572 CeedQFunctionAddInput(qf_ifunctionIn, "x", ncompx, CEED_EVAL_INTERP); 1573 CeedQFunctionAddOutput(qf_ifunctionIn, "v", ncompq, CEED_EVAL_INTERP); 1574 } 1575 1576 // Create CEED Operator for the whole domain 1577 if (!implicit) 1578 ierr = CreateOperatorForDomain(ceed, dm, user->op_rhs_vol, qf_rhsOut, qf_rhsIn, qf_setupSur, height, 1579 bc.noutflow, bc.outflow, bc.ninflow, bc.inflow, numP_Sur, numQ_Sur, qdatasizeSur, 1580 NqptsSur, basisxSur, basisqSur, &user->op_rhs); 1581 CHKERRQ(ierr); 1582 if (implicit) 1583 ierr = CreateOperatorForDomain(ceed, dm, user->op_ifunction_vol, qf_ifunctionOut, qf_ifunctionIn, qf_setupSur, height, 1584 bc.noutflow, bc.outflow, bc.ninflow, bc.inflow, numP_Sur, numQ_Sur, qdatasizeSur, 1585 NqptsSur, basisxSur, basisqSur, &user->op_ifunction); 1586 CHKERRQ(ierr); 1587 1588 //--------------------------------------------------------------------------------------// 1589 CeedQFunctionSetContext(qf_ics, &ctxSetup, sizeof ctxSetup); 1590 CeedScalar ctxNS[8] = {lambda, mu, k, cv, cp, g, Rd}; 1591 CeedScalar ctxIn[5] = {cv, cp, Rd, P_left, rho_left}; 1592 struct Advection2dContext_ ctxAdvection2d = { 1593 .CtauS = CtauS, 1594 .strong_form = strong_form, 1595 .stabilization = stab, 1596 }; 1597 switch (problemChoice) { 1598 case NS_DENSITY_CURRENT: 1599 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1600 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1601 sizeof ctxNS); 1602 break; 1603 case NS_ADVECTION: 1604 case NS_ADVECTION2D: 1605 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1606 sizeof ctxAdvection2d); 1607 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1608 sizeof ctxAdvection2d); 1609 if (qf_rhsOut) CeedQFunctionSetContext(qf_rhsOut, &ctxAdvection2d, 1610 sizeof ctxAdvection2d); 1611 if (qf_ifunctionOut) CeedQFunctionSetContext(qf_ifunctionOut, &ctxAdvection2d, 1612 sizeof ctxAdvection2d); 1613 if (qf_rhsIn) CeedQFunctionSetContext(qf_rhsIn, &ctxIn, sizeof ctxIn); 1614 if (qf_ifunctionIn) CeedQFunctionSetContext(qf_ifunctionIn, &ctxIn, sizeof ctxIn); 1615 } 1616 1617 // Set up PETSc context 1618 // Set up units structure 1619 units->meter = meter; 1620 units->kilogram = kilogram; 1621 units->second = second; 1622 units->Kelvin = Kelvin; 1623 units->Pascal = Pascal; 1624 units->JperkgK = JperkgK; 1625 units->mpersquareds = mpersquareds; 1626 units->WpermK = WpermK; 1627 units->kgpercubicm = kgpercubicm; 1628 units->kgpersquaredms = kgpersquaredms; 1629 units->Joulepercubicm = Joulepercubicm; 1630 1631 // Set up user structure 1632 user->comm = comm; 1633 user->outputfreq = outputfreq; 1634 user->contsteps = contsteps; 1635 user->units = units; 1636 user->dm = dm; 1637 user->dmviz = dmviz; 1638 user->interpviz = interpviz; 1639 user->ceed = ceed; 1640 1641 // Calculate qdata and ICs 1642 // Set up state global and local vectors 1643 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1644 1645 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1646 1647 // Apply Setup Ceed Operators 1648 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1649 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1650 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1651 user->M); CHKERRQ(ierr); 1652 1653 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1654 &ctxSetup, 0.0); CHKERRQ(ierr); 1655 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1656 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1657 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1658 // slower execution. 1659 Vec Qbc; 1660 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1661 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1662 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1663 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1664 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1665 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1666 ierr = PetscObjectComposeFunction((PetscObject)dm, 1667 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1668 CHKERRQ(ierr); 1669 } 1670 1671 MPI_Comm_rank(comm, &rank); 1672 if (!rank) {ierr = PetscMkdir(user->outputfolder); CHKERRQ(ierr);} 1673 // Gather initial Q values 1674 // In case of continuation of simulation, set up initial values from binary file 1675 if (contsteps) { // continue from existent solution 1676 PetscViewer viewer; 1677 char filepath[PETSC_MAX_PATH_LEN]; 1678 // Read input 1679 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1680 user->outputfolder); 1681 CHKERRQ(ierr); 1682 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1683 CHKERRQ(ierr); 1684 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1685 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1686 } 1687 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1688 1689 // Create and setup TS 1690 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1691 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1692 if (implicit) { 1693 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1694 if (user->op_ifunction) { 1695 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1696 } else { // Implicit integrators can fall back to using an RHSFunction 1697 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1698 } 1699 } else { 1700 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1701 "Problem does not provide RHSFunction"); 1702 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1703 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1704 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1705 } 1706 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1707 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1708 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1709 if (testChoice != TEST_NONE) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1710 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1711 ierr = TSAdaptSetStepLimits(adapt, 1712 1.e-12 * units->second, 1713 1.e2 * units->second); CHKERRQ(ierr); 1714 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1715 if (!contsteps) { // print initial condition 1716 if (testChoice == TEST_NONE) { 1717 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1718 } 1719 } else { // continue from time of last output 1720 PetscReal time; 1721 PetscInt count; 1722 PetscViewer viewer; 1723 char filepath[PETSC_MAX_PATH_LEN]; 1724 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1725 user->outputfolder); CHKERRQ(ierr); 1726 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1727 CHKERRQ(ierr); 1728 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1729 CHKERRQ(ierr); 1730 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1731 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1732 } 1733 if (testChoice == TEST_NONE) { 1734 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1735 } 1736 1737 // Solve 1738 start = MPI_Wtime(); 1739 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1740 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1741 cpu_time_used = MPI_Wtime() - start; 1742 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1743 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1744 comm); CHKERRQ(ierr); 1745 if (testChoice == TEST_NONE) { 1746 ierr = PetscPrintf(PETSC_COMM_WORLD, 1747 "Time taken for solution (sec): %g\n", 1748 (double)cpu_time_used); CHKERRQ(ierr); 1749 } 1750 1751 // Get error 1752 if (problem->non_zero_time && testChoice == TEST_NONE) { 1753 Vec Qexact, Qexactloc; 1754 PetscReal norm; 1755 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1756 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1757 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1758 1759 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1760 restrictq, &ctxSetup, ftime); CHKERRQ(ierr); 1761 1762 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1763 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1764 CeedVectorDestroy(&q0ceed); 1765 ierr = PetscPrintf(PETSC_COMM_WORLD, 1766 "Max Error: %g\n", 1767 (double)norm); CHKERRQ(ierr); 1768 // Clean up vectors 1769 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1770 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1771 } 1772 1773 // Output Statistics 1774 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1775 if (testChoice == TEST_NONE) { 1776 ierr = PetscPrintf(PETSC_COMM_WORLD, 1777 "Time integrator took %D time steps to reach final time %g\n", 1778 steps, (double)ftime); CHKERRQ(ierr); 1779 } 1780 // Output numerical values from command line 1781 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1782 1783 // compare reference solution values with current run 1784 if (testChoice != TEST_NONE) { 1785 PetscViewer viewer; 1786 // Read reference file 1787 Vec Qref; 1788 PetscReal error, Qrefnorm; 1789 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1790 ierr = PetscViewerBinaryOpen(comm, test->filepath, FILE_MODE_READ, &viewer); 1791 CHKERRQ(ierr); 1792 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1793 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1794 1795 // Compute error with respect to reference solution 1796 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1797 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1798 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1799 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1800 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1801 // Check error 1802 if (error > test->testtol) { 1803 ierr = PetscPrintf(PETSC_COMM_WORLD, 1804 "Test failed with error norm %g\n", 1805 (double)error); CHKERRQ(ierr); 1806 } 1807 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1808 } 1809 1810 // Clean up libCEED 1811 CeedVectorDestroy(&qdata); 1812 CeedVectorDestroy(&user->qceed); 1813 CeedVectorDestroy(&user->qdotceed); 1814 CeedVectorDestroy(&user->gceed); 1815 CeedVectorDestroy(&xcorners); 1816 CeedBasisDestroy(&basisq); 1817 CeedBasisDestroy(&basisx); 1818 CeedBasisDestroy(&basisxc); 1819 CeedElemRestrictionDestroy(&restrictq); 1820 CeedElemRestrictionDestroy(&restrictx); 1821 CeedElemRestrictionDestroy(&restrictqdi); 1822 CeedQFunctionDestroy(&qf_setupVol); 1823 CeedQFunctionDestroy(&qf_ics); 1824 CeedQFunctionDestroy(&qf_rhsVol); 1825 CeedQFunctionDestroy(&qf_ifunctionVol); 1826 CeedOperatorDestroy(&op_setupVol); 1827 CeedOperatorDestroy(&op_ics); 1828 CeedOperatorDestroy(&user->op_rhs_vol); 1829 CeedOperatorDestroy(&user->op_ifunction_vol); 1830 CeedDestroy(&ceed); 1831 CeedBasisDestroy(&basisqSur); 1832 CeedBasisDestroy(&basisxSur); 1833 CeedBasisDestroy(&basisxcSur); 1834 CeedQFunctionDestroy(&qf_setupSur); 1835 CeedQFunctionDestroy(&qf_rhsOut); 1836 CeedQFunctionDestroy(&qf_ifunctionOut); 1837 CeedQFunctionDestroy(&qf_rhsIn); 1838 CeedQFunctionDestroy(&qf_ifunctionIn); 1839 CeedOperatorDestroy(&user->op_rhs); 1840 CeedOperatorDestroy(&user->op_ifunction); 1841 1842 // Clean up PETSc 1843 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1844 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1845 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1846 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1847 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1848 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1849 ierr = PetscFree(units); CHKERRQ(ierr); 1850 ierr = PetscFree(user); CHKERRQ(ierr); 1851 return PetscFinalize(); 1852 } 1853 1854