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