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