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