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