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.; // 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 ctxEulerData->stabilization = stab; 1635 user->ctxEulerData = ctxEulerData; 1636 1637 CeedQFunctionContextCreate(ceed, &ctxEuler); 1638 CeedQFunctionContextSetData(ctxEuler, CEED_MEM_HOST, CEED_USE_POINTER, 1639 sizeof *ctxEulerData, ctxEulerData); 1640 1641 switch (problemChoice) { 1642 case NS_DENSITY_CURRENT: 1643 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxNS); 1644 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxNS); 1645 break; 1646 case NS_ADVECTION: 1647 case NS_ADVECTION2D: 1648 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxAdvection2d); 1649 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxAdvection2d); 1650 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, ctxSurface); 1651 case NS_EULER_VORTEX: 1652 if (qf_ics) CeedQFunctionSetContext(qf_ics, ctxEuler); 1653 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxEuler); 1654 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, ctxEuler); 1655 } 1656 1657 // Set up PETSc context 1658 // Set up units structure 1659 units->meter = meter; 1660 units->kilogram = kilogram; 1661 units->second = second; 1662 units->Kelvin = Kelvin; 1663 units->Pascal = Pascal; 1664 units->JperkgK = JperkgK; 1665 units->mpersquareds = mpersquareds; 1666 units->WpermK = WpermK; 1667 units->kgpercubicm = kgpercubicm; 1668 units->kgpersquaredms = kgpersquaredms; 1669 units->Joulepercubicm = Joulepercubicm; 1670 units->Joule = Joule; 1671 1672 // Set up user structure 1673 user->comm = comm; 1674 user->outputfreq = outputfreq; 1675 user->contsteps = contsteps; 1676 user->units = units; 1677 user->dm = dm; 1678 user->dmviz = dmviz; 1679 user->interpviz = interpviz; 1680 user->ceed = ceed; 1681 1682 // Calculate qdata and ICs 1683 // Set up state global and local vectors 1684 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1685 1686 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1687 1688 // Apply Setup Ceed Operators 1689 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1690 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1691 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1692 user->M); CHKERRQ(ierr); 1693 1694 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1695 ctxSetup, 0.0); CHKERRQ(ierr); 1696 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1697 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1698 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1699 // slower execution. 1700 Vec Qbc; 1701 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1702 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1703 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1704 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1705 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1706 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1707 ierr = PetscObjectComposeFunction((PetscObject)dm, 1708 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1709 CHKERRQ(ierr); 1710 } 1711 1712 MPI_Comm_rank(comm, &rank); 1713 if (!rank) {ierr = PetscMkdir(user->outputdir); CHKERRQ(ierr);} 1714 // Gather initial Q values 1715 // In case of continuation of simulation, set up initial values from binary file 1716 if (contsteps) { // continue from existent solution 1717 PetscViewer viewer; 1718 char filepath[PETSC_MAX_PATH_LEN]; 1719 // Read input 1720 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1721 user->outputdir); 1722 CHKERRQ(ierr); 1723 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1724 CHKERRQ(ierr); 1725 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1726 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1727 } 1728 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1729 1730 // Create and setup TS 1731 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1732 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1733 if (implicit) { 1734 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1735 if (user->op_ifunction) { 1736 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1737 } else { // Implicit integrators can fall back to using an RHSFunction 1738 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1739 } 1740 } else { 1741 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1742 "Problem does not provide RHSFunction"); 1743 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1744 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1745 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1746 } 1747 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1748 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1749 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1750 if (test) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1751 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1752 ierr = TSAdaptSetStepLimits(adapt, 1753 1.e-12 * units->second, 1754 1.e2 * units->second); CHKERRQ(ierr); 1755 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1756 if (!contsteps) { // print initial condition 1757 if (!test) { 1758 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1759 } 1760 } else { // continue from time of last output 1761 PetscReal time; 1762 PetscInt count; 1763 PetscViewer viewer; 1764 char filepath[PETSC_MAX_PATH_LEN]; 1765 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1766 user->outputdir); CHKERRQ(ierr); 1767 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1768 CHKERRQ(ierr); 1769 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1770 CHKERRQ(ierr); 1771 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1772 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1773 } 1774 if (!test) { 1775 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1776 } 1777 1778 // Solve 1779 start = MPI_Wtime(); 1780 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1781 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1782 cpu_time_used = MPI_Wtime() - start; 1783 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1784 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1785 comm); CHKERRQ(ierr); 1786 if (!test) { 1787 ierr = PetscPrintf(PETSC_COMM_WORLD, 1788 "Time taken for solution (sec): %g\n", 1789 (double)cpu_time_used); CHKERRQ(ierr); 1790 } 1791 1792 // Get error 1793 if (problem->non_zero_time && !test) { 1794 Vec Qexact, Qexactloc; 1795 PetscReal rel_error, norm_error, norm_exact; 1796 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1797 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1798 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1799 1800 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1801 restrictq, ctxSetup, ftime); CHKERRQ(ierr); 1802 ierr = VecNorm(Qexact, NORM_1, &norm_exact); CHKERRQ(ierr); 1803 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1804 ierr = VecNorm(Q, NORM_1, &norm_error); CHKERRQ(ierr); 1805 rel_error = norm_error / norm_exact; 1806 CeedVectorDestroy(&q0ceed); 1807 ierr = PetscPrintf(PETSC_COMM_WORLD, 1808 "Relative Error: %g\n", 1809 (double)rel_error); CHKERRQ(ierr); 1810 // Clean up vectors 1811 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1812 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1813 } 1814 1815 // Output Statistics 1816 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1817 if (!test) { 1818 ierr = PetscPrintf(PETSC_COMM_WORLD, 1819 "Time integrator took %D time steps to reach final time %g\n", 1820 steps, (double)ftime); CHKERRQ(ierr); 1821 } 1822 1823 1824 // Compare reference solution values with current test run for CI 1825 if (test) { 1826 PetscViewer viewer; 1827 // Read reference file 1828 Vec Qref; 1829 PetscReal error, Qrefnorm; 1830 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1831 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1832 CHKERRQ(ierr); 1833 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1834 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1835 1836 // Compute error with respect to reference solution 1837 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1838 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1839 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1840 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1841 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1842 // Check error 1843 if (error > testtol) { 1844 ierr = PetscPrintf(PETSC_COMM_WORLD, 1845 "Test failed with error norm %g\n", 1846 (double)error); CHKERRQ(ierr); 1847 } 1848 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1849 } 1850 1851 // Clean up libCEED 1852 CeedVectorDestroy(&qdata); 1853 CeedVectorDestroy(&user->qceed); 1854 CeedVectorDestroy(&user->qdotceed); 1855 CeedVectorDestroy(&user->gceed); 1856 CeedVectorDestroy(&xcorners); 1857 CeedBasisDestroy(&basisq); 1858 CeedBasisDestroy(&basisx); 1859 CeedBasisDestroy(&basisxc); 1860 CeedElemRestrictionDestroy(&restrictq); 1861 CeedElemRestrictionDestroy(&restrictx); 1862 CeedElemRestrictionDestroy(&restrictqdi); 1863 CeedQFunctionDestroy(&qf_setupVol); 1864 CeedQFunctionDestroy(&qf_ics); 1865 CeedQFunctionDestroy(&qf_rhsVol); 1866 CeedQFunctionDestroy(&qf_ifunctionVol); 1867 CeedQFunctionContextDestroy(&ctxSetup); 1868 CeedQFunctionContextDestroy(&ctxNS); 1869 CeedQFunctionContextDestroy(&ctxAdvection2d); 1870 CeedQFunctionContextDestroy(&ctxSurface); 1871 CeedQFunctionContextDestroy(&ctxEuler); 1872 CeedOperatorDestroy(&op_setupVol); 1873 CeedOperatorDestroy(&op_ics); 1874 CeedOperatorDestroy(&user->op_rhs_vol); 1875 CeedOperatorDestroy(&user->op_ifunction_vol); 1876 CeedDestroy(&ceed); 1877 CeedBasisDestroy(&basisqSur); 1878 CeedBasisDestroy(&basisxSur); 1879 CeedBasisDestroy(&basisxcSur); 1880 CeedQFunctionDestroy(&qf_setupSur); 1881 CeedQFunctionDestroy(&qf_applySur); 1882 CeedOperatorDestroy(&user->op_rhs); 1883 CeedOperatorDestroy(&user->op_ifunction); 1884 1885 // Clean up PETSc 1886 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1887 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1888 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1889 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1890 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1891 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1892 ierr = PetscFree(units); CHKERRQ(ierr); 1893 ierr = PetscFree(user); CHKERRQ(ierr); 1894 ierr = PetscFree(ctxEulerData); CHKERRQ(ierr); 1895 return PetscFinalize(); 1896 } 1897