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