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) { 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, ctxSetupData); 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 .wind_type = wind_type, 1292 }; 1293 1294 // Create the mesh 1295 { 1296 const PetscReal scale[3] = {lx, ly, lz}; 1297 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1298 NULL, PETSC_TRUE, &dm); 1299 CHKERRQ(ierr); 1300 } 1301 1302 // Distribute the mesh over processes 1303 { 1304 DM dmDist = NULL; 1305 PetscPartitioner part; 1306 1307 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1308 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1309 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1310 if (dmDist) { 1311 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1312 dm = dmDist; 1313 } 1314 } 1315 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1316 1317 // Setup DM 1318 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1319 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1320 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetupData); 1321 CHKERRQ(ierr); 1322 1323 // Refine DM for high-order viz 1324 dmviz = NULL; 1325 interpviz = NULL; 1326 if (viz_refine) { 1327 DM dmhierarchy[viz_refine+1]; 1328 1329 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1330 dmhierarchy[0] = dm; 1331 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1332 Mat interp_next; 1333 1334 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1335 CHKERRQ(ierr); 1336 ierr = DMClearDS(dmhierarchy[i+1]); CHKERRQ(ierr); 1337 ierr = DMClearFields(dmhierarchy[i+1]); CHKERRQ(ierr); 1338 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1339 d = (d + 1) / 2; 1340 if (i + 1 == viz_refine) d = 1; 1341 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetupData); 1342 CHKERRQ(ierr); 1343 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1344 &interp_next, NULL); CHKERRQ(ierr); 1345 if (!i) interpviz = interp_next; 1346 else { 1347 Mat C; 1348 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1349 PETSC_DECIDE, &C); CHKERRQ(ierr); 1350 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1351 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1352 interpviz = C; 1353 } 1354 } 1355 for (PetscInt i=1; i<viz_refine; i++) { 1356 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1357 } 1358 dmviz = dmhierarchy[viz_refine]; 1359 } 1360 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1361 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1362 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1363 lnodes /= ncompq; 1364 1365 // Initialize CEED 1366 CeedInit(ceedresource, &ceed); 1367 // Set memtype 1368 CeedMemType memtypebackend; 1369 CeedGetPreferredMemType(ceed, &memtypebackend); 1370 // Check memtype compatibility 1371 if (!setmemtyperequest) 1372 memtyperequested = memtypebackend; 1373 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1374 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1375 "PETSc was not built with CUDA. " 1376 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1377 1378 // Set number of 1D nodes and quadrature points 1379 numP = degree + 1; 1380 numQ = numP + qextra; 1381 1382 // Print summary 1383 if (!test) { 1384 CeedInt gdofs, odofs; 1385 int comm_size; 1386 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1387 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1388 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1389 gnodes = gdofs/ncompq; 1390 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1391 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1392 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1393 const char *usedresource; 1394 CeedGetResource(ceed, &usedresource); 1395 1396 ierr = PetscPrintf(comm, 1397 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1398 " rank(s) : %d\n" 1399 " Problem:\n" 1400 " Problem Name : %s\n" 1401 " Stabilization : %s\n" 1402 " PETSc:\n" 1403 " Box Faces : %s\n" 1404 " libCEED:\n" 1405 " libCEED Backend : %s\n" 1406 " libCEED Backend MemType : %s\n" 1407 " libCEED User Requested MemType : %s\n" 1408 " Mesh:\n" 1409 " Number of 1D Basis Nodes (P) : %d\n" 1410 " Number of 1D Quadrature Points (Q) : %d\n" 1411 " Global DoFs : %D\n" 1412 " Owned DoFs : %D\n" 1413 " DoFs per node : %D\n" 1414 " Global nodes : %D\n" 1415 " Owned nodes : %D\n", 1416 comm_size, problemTypes[problemChoice], 1417 StabilizationTypes[stab], box_faces_str, usedresource, 1418 CeedMemTypes[memtypebackend], 1419 (setmemtyperequest) ? 1420 CeedMemTypes[memtyperequested] : "none", 1421 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1422 CHKERRQ(ierr); 1423 } 1424 1425 // Set up global mass vector 1426 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1427 1428 // Set up libCEED 1429 // CEED Bases 1430 CeedInit(ceedresource, &ceed); 1431 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1432 &basisq); 1433 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1434 &basisx); 1435 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1436 CEED_GAUSS_LOBATTO, &basisxc); 1437 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1438 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1439 CHKERRQ(ierr); 1440 1441 // CEED Restrictions 1442 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1443 qdatasizeVol, &restrictq, &restrictx, 1444 &restrictqdi); CHKERRQ(ierr); 1445 1446 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1447 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1448 1449 // Create the CEED vectors that will be needed in setup 1450 CeedInt NqptsVol; 1451 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1452 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1453 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1454 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1455 1456 // Create the Q-function that builds the quadrature data for the NS operator 1457 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1458 &qf_setupVol); 1459 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1460 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1461 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1462 1463 // Create the Q-function that sets the ICs of the operator 1464 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1465 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1466 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1467 1468 qf_rhsVol = NULL; 1469 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1470 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1471 problem->applyVol_rhs_loc, &qf_rhsVol); 1472 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1473 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1474 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1475 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1476 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1477 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1478 } 1479 1480 qf_ifunctionVol = NULL; 1481 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1482 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1483 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1484 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1485 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1486 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1487 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1488 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1489 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1490 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1491 } 1492 1493 // Create the operator that builds the quadrature data for the NS operator 1494 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1495 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1496 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1497 basisx, CEED_VECTOR_NONE); 1498 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1499 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1500 1501 // Create the operator that sets the ICs 1502 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1503 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1504 CeedOperatorSetField(op_ics, "q0", restrictq, 1505 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1506 1507 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1508 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1509 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1510 1511 if (qf_rhsVol) { // Create the RHS physics operator 1512 CeedOperator op; 1513 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1514 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1515 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1516 CeedOperatorSetField(op, "qdata", restrictqdi, 1517 CEED_BASIS_COLLOCATED, qdata); 1518 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1519 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1520 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1521 user->op_rhs_vol = op; 1522 } 1523 1524 if (qf_ifunctionVol) { // Create the IFunction operator 1525 CeedOperator op; 1526 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1527 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1528 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1529 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1530 CeedOperatorSetField(op, "qdata", restrictqdi, 1531 CEED_BASIS_COLLOCATED, qdata); 1532 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1533 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1534 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1535 user->op_ifunction_vol = op; 1536 } 1537 1538 // Set up CEED for the boundaries 1539 CeedInt height = 1; 1540 CeedInt dimSur = dim - height; 1541 CeedInt numP_Sur = degree + 1; 1542 CeedInt numQ_Sur = numP_Sur + qextraSur; 1543 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1544 CeedBasis basisxSur, basisxcSur, basisqSur; 1545 CeedInt NqptsSur; 1546 CeedQFunction qf_setupSur, qf_applySur; 1547 1548 // CEED bases for the boundaries 1549 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, 1550 CEED_GAUSS, 1551 &basisqSur); 1552 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1553 &basisxSur); 1554 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1555 CEED_GAUSS_LOBATTO, &basisxcSur); 1556 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1557 1558 // Create the Q-function that builds the quadrature data for the Surface operator 1559 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1560 &qf_setupSur); 1561 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1562 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1563 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1564 1565 // Creat Q-Function for Boundaries 1566 // -- Defined for Advection(2d) test cases 1567 qf_applySur = NULL; 1568 if (problem->applySur) { 1569 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1570 problem->applySur_loc, &qf_applySur); 1571 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1572 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1573 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1574 CeedQFunctionAddOutput(qf_applySur, "v", ncompq, CEED_EVAL_INTERP); 1575 } 1576 1577 // Create CEED Operator for the whole domain 1578 if (!implicit) 1579 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1580 user->op_rhs_vol, qf_applySur, 1581 qf_setupSur, height, numP_Sur, numQ_Sur, 1582 qdatasizeSur, NqptsSur, basisxSur, 1583 basisqSur, &user->op_rhs); 1584 CHKERRQ(ierr); 1585 if (implicit) 1586 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1587 user->op_ifunction_vol, qf_applySur, 1588 qf_setupSur, height, numP_Sur, numQ_Sur, 1589 qdatasizeSur, NqptsSur, basisxSur, 1590 basisqSur, &user->op_ifunction); 1591 CHKERRQ(ierr); 1592 // Set up contex for QFunctions 1593 CeedQFunctionContextCreate(ceed, &ctxSetup); 1594 CeedQFunctionContextSetData(ctxSetup, CEED_MEM_HOST, CEED_USE_POINTER, 1595 sizeof ctxSetupData, &ctxSetupData); 1596 if (qf_ics && problemChoice != NS_EULER_VORTEX) 1597 CeedQFunctionSetContext(qf_ics, ctxSetup); 1598 1599 CeedScalar ctxNSData[8] = {lambda, mu, k, cv, cp, g, Rd}; 1600 CeedQFunctionContextCreate(ceed, &ctxNS); 1601 CeedQFunctionContextSetData(ctxNS, CEED_MEM_HOST, CEED_USE_POINTER, 1602 sizeof ctxNSData, &ctxNSData); 1603 1604 struct Advection2dContext_ ctxAdvection2dData = { 1605 .CtauS = CtauS, 1606 .strong_form = strong_form, 1607 .stabilization = stab, 1608 }; 1609 CeedQFunctionContextCreate(ceed, &ctxAdvection2d); 1610 CeedQFunctionContextSetData(ctxAdvection2d, CEED_MEM_HOST, CEED_USE_POINTER, 1611 sizeof ctxAdvection2dData, &ctxAdvection2dData); 1612 1613 struct SurfaceContext_ ctxSurfaceData = { 1614 .E_wind = E_wind, 1615 .strong_form = strong_form, 1616 .implicit = implicit, 1617 }; 1618 CeedQFunctionContextCreate(ceed, &ctxSurface); 1619 CeedQFunctionContextSetData(ctxSurface, CEED_MEM_HOST, CEED_USE_POINTER, 1620 sizeof ctxSurfaceData, &ctxSurfaceData); 1621 1622 // Set up ctxEulerData structure 1623 ctxEulerData->time = 0.; 1624 ctxEulerData->currentTime = 0.; 1625 ctxEulerData->euler_test = euler_test; 1626 ctxEulerData->center[0] = center[0]; 1627 ctxEulerData->center[1] = center[1]; 1628 ctxEulerData->center[2] = center[2]; 1629 ctxEulerData->vortex_strength = vortex_strength; 1630 ctxEulerData->etv_mean_velocity[0] = etv_mean_velocity[0]; 1631 ctxEulerData->etv_mean_velocity[1] = etv_mean_velocity[1]; 1632 ctxEulerData->etv_mean_velocity[2] = etv_mean_velocity[2]; 1633 ctxEulerData->T_inlet = T_inlet; 1634 ctxEulerData->P_inlet = P_inlet; 1635 ctxEulerData->stabilization = stab; 1636 ctxEulerData->implicit = implicit; 1637 user->ctxEulerData = ctxEulerData; 1638 1639 CeedQFunctionContextCreate(ceed, &ctxEuler); 1640 CeedQFunctionContextSetData(ctxEuler, CEED_MEM_HOST, CEED_USE_POINTER, 1641 sizeof *ctxEulerData, ctxEulerData); 1642 1643 switch (problemChoice) { 1644 case NS_DENSITY_CURRENT: 1645 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxNS); 1646 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxNS); 1647 break; 1648 case NS_ADVECTION: 1649 case NS_ADVECTION2D: 1650 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxAdvection2d); 1651 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxAdvection2d); 1652 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, ctxSurface); 1653 case NS_EULER_VORTEX: 1654 if (qf_ics) CeedQFunctionSetContext(qf_ics, ctxEuler); 1655 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxEuler); 1656 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxEuler); 1657 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, ctxEuler); 1658 } 1659 1660 // Set up PETSc context 1661 // Set up units structure 1662 units->meter = meter; 1663 units->kilogram = kilogram; 1664 units->second = second; 1665 units->Kelvin = Kelvin; 1666 units->Pascal = Pascal; 1667 units->JperkgK = JperkgK; 1668 units->mpersquareds = mpersquareds; 1669 units->WpermK = WpermK; 1670 units->kgpercubicm = kgpercubicm; 1671 units->kgpersquaredms = kgpersquaredms; 1672 units->Joulepercubicm = Joulepercubicm; 1673 units->Joule = Joule; 1674 1675 // Set up user structure 1676 user->comm = comm; 1677 user->outputfreq = outputfreq; 1678 user->contsteps = contsteps; 1679 user->units = units; 1680 user->dm = dm; 1681 user->dmviz = dmviz; 1682 user->interpviz = interpviz; 1683 user->ceed = ceed; 1684 1685 // Calculate qdata and ICs 1686 // Set up state global and local vectors 1687 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1688 1689 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1690 1691 // Apply Setup Ceed Operators 1692 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1693 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1694 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1695 user->M); CHKERRQ(ierr); 1696 1697 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1698 ctxSetup, 0.0); CHKERRQ(ierr); 1699 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1700 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1701 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1702 // slower execution. 1703 Vec Qbc; 1704 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1705 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1706 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1707 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1708 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1709 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1710 ierr = PetscObjectComposeFunction((PetscObject)dm, 1711 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1712 CHKERRQ(ierr); 1713 } 1714 1715 MPI_Comm_rank(comm, &rank); 1716 if (!rank) {ierr = PetscMkdir(user->outputdir); CHKERRQ(ierr);} 1717 // Gather initial Q values 1718 // In case of continuation of simulation, set up initial values from binary file 1719 if (contsteps) { // continue from existent solution 1720 PetscViewer viewer; 1721 char filepath[PETSC_MAX_PATH_LEN]; 1722 // Read input 1723 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1724 user->outputdir); 1725 CHKERRQ(ierr); 1726 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1727 CHKERRQ(ierr); 1728 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1729 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1730 } 1731 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1732 1733 // Create and setup TS 1734 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1735 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1736 if (implicit) { 1737 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1738 if (user->op_ifunction) { 1739 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1740 } else { // Implicit integrators can fall back to using an RHSFunction 1741 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1742 } 1743 } else { 1744 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1745 "Problem does not provide RHSFunction"); 1746 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1747 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1748 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1749 } 1750 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1751 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1752 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1753 if (test) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1754 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1755 ierr = TSAdaptSetStepLimits(adapt, 1756 1.e-12 * units->second, 1757 1.e2 * units->second); CHKERRQ(ierr); 1758 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1759 if (!contsteps) { // print initial condition 1760 if (!test) { 1761 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1762 } 1763 } else { // continue from time of last output 1764 PetscReal time; 1765 PetscInt count; 1766 PetscViewer viewer; 1767 char filepath[PETSC_MAX_PATH_LEN]; 1768 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1769 user->outputdir); CHKERRQ(ierr); 1770 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1771 CHKERRQ(ierr); 1772 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1773 CHKERRQ(ierr); 1774 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1775 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1776 } 1777 if (!test) { 1778 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1779 } 1780 1781 // Solve 1782 start = MPI_Wtime(); 1783 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1784 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1785 cpu_time_used = MPI_Wtime() - start; 1786 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1787 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1788 comm); CHKERRQ(ierr); 1789 if (!test) { 1790 ierr = PetscPrintf(PETSC_COMM_WORLD, 1791 "Time taken for solution (sec): %g\n", 1792 (double)cpu_time_used); CHKERRQ(ierr); 1793 } 1794 1795 // Get error 1796 if (problem->non_zero_time && !test) { 1797 Vec Qexact, Qexactloc; 1798 PetscReal rel_error, norm_error, norm_exact; 1799 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1800 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1801 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1802 1803 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1804 restrictq, ctxSetup, ftime); CHKERRQ(ierr); 1805 ierr = VecNorm(Qexact, NORM_1, &norm_exact); CHKERRQ(ierr); 1806 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1807 ierr = VecNorm(Q, NORM_1, &norm_error); CHKERRQ(ierr); 1808 rel_error = norm_error / norm_exact; 1809 CeedVectorDestroy(&q0ceed); 1810 ierr = PetscPrintf(PETSC_COMM_WORLD, 1811 "Relative Error: %g\n", 1812 (double)rel_error); CHKERRQ(ierr); 1813 // Clean up vectors 1814 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1815 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1816 } 1817 1818 // Output Statistics 1819 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1820 if (!test) { 1821 ierr = PetscPrintf(PETSC_COMM_WORLD, 1822 "Time integrator took %D time steps to reach final time %g\n", 1823 steps, (double)ftime); CHKERRQ(ierr); 1824 } 1825 // Output numerical values from command line 1826 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1827 1828 // Compare reference solution values with current test run for CI 1829 if (test) { 1830 PetscViewer viewer; 1831 // Read reference file 1832 Vec Qref; 1833 PetscReal error, Qrefnorm; 1834 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1835 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1836 CHKERRQ(ierr); 1837 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1838 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1839 1840 // Compute error with respect to reference solution 1841 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1842 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1843 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1844 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1845 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1846 // Check error 1847 if (error > testtol) { 1848 ierr = PetscPrintf(PETSC_COMM_WORLD, 1849 "Test failed with error norm %g\n", 1850 (double)error); CHKERRQ(ierr); 1851 } 1852 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1853 } 1854 1855 // Clean up libCEED 1856 CeedVectorDestroy(&qdata); 1857 CeedVectorDestroy(&user->qceed); 1858 CeedVectorDestroy(&user->qdotceed); 1859 CeedVectorDestroy(&user->gceed); 1860 CeedVectorDestroy(&xcorners); 1861 CeedBasisDestroy(&basisq); 1862 CeedBasisDestroy(&basisx); 1863 CeedBasisDestroy(&basisxc); 1864 CeedElemRestrictionDestroy(&restrictq); 1865 CeedElemRestrictionDestroy(&restrictx); 1866 CeedElemRestrictionDestroy(&restrictqdi); 1867 CeedQFunctionDestroy(&qf_setupVol); 1868 CeedQFunctionDestroy(&qf_ics); 1869 CeedQFunctionDestroy(&qf_rhsVol); 1870 CeedQFunctionDestroy(&qf_ifunctionVol); 1871 CeedQFunctionContextDestroy(&ctxSetup); 1872 CeedQFunctionContextDestroy(&ctxNS); 1873 CeedQFunctionContextDestroy(&ctxAdvection2d); 1874 CeedQFunctionContextDestroy(&ctxSurface); 1875 CeedQFunctionContextDestroy(&ctxEuler); 1876 CeedOperatorDestroy(&op_setupVol); 1877 CeedOperatorDestroy(&op_ics); 1878 CeedOperatorDestroy(&user->op_rhs_vol); 1879 CeedOperatorDestroy(&user->op_ifunction_vol); 1880 CeedDestroy(&ceed); 1881 CeedBasisDestroy(&basisqSur); 1882 CeedBasisDestroy(&basisxSur); 1883 CeedBasisDestroy(&basisxcSur); 1884 CeedQFunctionDestroy(&qf_setupSur); 1885 CeedQFunctionDestroy(&qf_applySur); 1886 CeedOperatorDestroy(&user->op_rhs); 1887 CeedOperatorDestroy(&user->op_ifunction); 1888 1889 // Clean up PETSc 1890 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1891 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1892 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1893 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1894 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1895 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1896 ierr = PetscFree(units); CHKERRQ(ierr); 1897 ierr = PetscFree(user); CHKERRQ(ierr); 1898 ierr = PetscFree(ctxEulerData); CHKERRQ(ierr); 1899 return PetscFinalize(); 1900 } 1901