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