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 for (CeedInt i=0; i<nFace; i++) { 454 ierr = GetRestrictionForDomain(ceed, dm, height, domainLabel, i+1, numP_Sur, 455 numQ_Sur, qdatasizeSur, &restrictqSur[i], 456 &restrictxSur[i], &restrictqdiSur[i]); 457 CHKERRQ(ierr); 458 // Create the CEED vectors that will be needed in Boundary setup 459 CeedElemRestrictionGetNumElements(restrictqSur[i], &localNelemSur[i]); 460 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, 461 &qdataSur[i]); 462 // Create the operator that builds the quadrature data for the Boundary operator 463 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur[i]); 464 CeedOperatorSetField(op_setupSur[i], "dx", restrictxSur[i], basisxSur, 465 CEED_VECTOR_ACTIVE); 466 CeedOperatorSetField(op_setupSur[i], "weight", CEED_ELEMRESTRICTION_NONE, 467 basisxSur, CEED_VECTOR_NONE); 468 CeedOperatorSetField(op_setupSur[i], "qdataSur", restrictqdiSur[i], 469 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 470 // Create Boundary operator 471 CeedOperatorCreate(ceed, qf_applySur, NULL, NULL, &op_applySur[i]); 472 CeedOperatorSetField(op_applySur[i], "q", restrictqSur[i], basisqSur, 473 CEED_VECTOR_ACTIVE); 474 CeedOperatorSetField(op_applySur[i], "qdataSur", restrictqdiSur[i], 475 CEED_BASIS_COLLOCATED, qdataSur[i]); 476 CeedOperatorSetField(op_applySur[i], "x", restrictxSur[i], basisxSur, 477 xcorners); 478 CeedOperatorSetField(op_applySur[i], "v", restrictqSur[i], basisqSur, 479 CEED_VECTOR_ACTIVE); 480 // Apply CEED operator for Boundary setup 481 CeedOperatorApply(op_setupSur[i], xcorners, qdataSur[i], 482 CEED_REQUEST_IMMEDIATE); 483 // --Apply Sub-Operator for the Boundary 484 CeedCompositeOperatorAddSub(*op_apply, op_applySur[i]); 485 } 486 CeedVectorDestroy(&xcorners); 487 } 488 489 if (problemChoice == NS_EULER_VORTEX) { 490 // All faces are slip, except for 5 and 6 491 bc->nwall = 0; 492 bc->nslip[0] = 0; 493 bc->nslip[1] = 2; bc->slips[1][0] = 3; bc->slips[1][1] = 4; 494 bc->nslip[2] = 2; bc->slips[2][0] = 1; bc->slips[2][1] = 2; 495 496 // Create CEED Operator for each boundary face 497 PetscInt localNelemSur[2]; 498 CeedVector qdataSur[2]; 499 CeedQFunction qf_apply; 500 CeedOperator op_setupSur[2], op_applySur[2]; 501 CeedElemRestriction restrictxSur[2], restrictqSur[2], restrictqdiSur[2]; 502 for (CeedInt i=0; i<2; i++) { 503 if (i == 0) qf_apply = qf_applyIn; 504 if (i == 1) qf_apply = qf_applyOut; 505 ierr = GetRestrictionForDomain(ceed, dm, height, domainLabel, i+5, numP_Sur, 506 numQ_Sur, qdatasizeSur, &restrictqSur[i], 507 &restrictxSur[i], &restrictqdiSur[i]); 508 CHKERRQ(ierr); 509 // Create the CEED vectors that will be needed in Boundary setup 510 CeedElemRestrictionGetNumElements(restrictqSur[i], &localNelemSur[i]); 511 CeedVectorCreate(ceed, qdatasizeSur*localNelemSur[i]*NqptsSur, 512 &qdataSur[i]); 513 // Create the operator that builds the quadrature data for the Boundary operator 514 CeedOperatorCreate(ceed, qf_setupSur, NULL, NULL, &op_setupSur[i]); 515 CeedOperatorSetField(op_setupSur[i], "dx", restrictxSur[i], basisxSur, 516 CEED_VECTOR_ACTIVE); 517 CeedOperatorSetField(op_setupSur[i], "weight", CEED_ELEMRESTRICTION_NONE, 518 basisxSur, CEED_VECTOR_NONE); 519 CeedOperatorSetField(op_setupSur[i], "qdataSur", restrictqdiSur[i], 520 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 521 // Create Boundary operator 522 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_applySur[i]); 523 CeedOperatorSetField(op_applySur[i], "q", restrictqSur[i], basisqSur, 524 CEED_VECTOR_ACTIVE); 525 CeedOperatorSetField(op_applySur[i], "qdataSur", restrictqdiSur[i], 526 CEED_BASIS_COLLOCATED, qdataSur[i]); 527 CeedOperatorSetField(op_applySur[i], "x", restrictxSur[i], basisxSur, 528 xcorners); 529 CeedOperatorSetField(op_applySur[i], "v", restrictqSur[i], basisqSur, 530 CEED_VECTOR_ACTIVE); 531 // Apply CEED operator for Boundary setup 532 CeedOperatorApply(op_setupSur[i], xcorners, qdataSur[i], 533 CEED_REQUEST_IMMEDIATE); 534 // --Apply Sub-Operator for the Boundary 535 CeedCompositeOperatorAddSub(*op_apply, op_applySur[i]); 536 } 537 CeedVectorDestroy(&xcorners); 538 539 } 540 PetscFunctionReturn(0); 541 } 542 543 static int CreateVectorFromPetscVec(Ceed ceed, Vec p, CeedVector *v) { 544 PetscErrorCode ierr; 545 PetscInt m; 546 547 PetscFunctionBeginUser; 548 ierr = VecGetLocalSize(p, &m); CHKERRQ(ierr); 549 ierr = CeedVectorCreate(ceed, m, v); CHKERRQ(ierr); 550 PetscFunctionReturn(0); 551 } 552 553 static int VectorPlacePetscVec(CeedVector c, Vec p) { 554 PetscErrorCode ierr; 555 PetscInt mceed, mpetsc; 556 PetscScalar *a; 557 558 PetscFunctionBeginUser; 559 ierr = CeedVectorGetLength(c, &mceed); CHKERRQ(ierr); 560 ierr = VecGetLocalSize(p, &mpetsc); CHKERRQ(ierr); 561 if (mceed != mpetsc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 562 "Cannot place PETSc Vec of length %D in CeedVector of length %D", 563 mpetsc, mceed); 564 ierr = VecGetArray(p, &a); CHKERRQ(ierr); 565 CeedVectorSetArray(c, CEED_MEM_HOST, CEED_USE_POINTER, a); 566 PetscFunctionReturn(0); 567 } 568 569 static PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 570 PetscBool insertEssential, Vec Qloc, PetscReal time, Vec faceGeomFVM, 571 Vec cellGeomFVM, Vec gradFVM) { 572 PetscErrorCode ierr; 573 Vec Qbc; 574 575 PetscFunctionBegin; 576 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 577 ierr = VecAXPY(Qloc, 1., Qbc); CHKERRQ(ierr); 578 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 579 PetscFunctionReturn(0); 580 } 581 582 // This is the RHS of the ODE, given as u_t = G(t,u) 583 // This function takes in a state vector Q and writes into G 584 static PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *userData) { 585 PetscErrorCode ierr; 586 User user = *(User *)userData; 587 PetscScalar *q, *g; 588 Vec Qloc, Gloc; 589 590 // Global-to-local 591 PetscFunctionBeginUser; 592 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 593 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 594 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 595 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 596 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 597 NULL, NULL, NULL); CHKERRQ(ierr); 598 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 599 600 // Ceed Vectors 601 ierr = VecGetArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 602 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 603 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); 604 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 605 606 // Apply CEED operator 607 CeedOperatorApply(user->op_rhs, user->qceed, user->gceed, 608 CEED_REQUEST_IMMEDIATE); 609 610 // Restore vectors 611 ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); 612 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 613 614 ierr = VecZeroEntries(G); CHKERRQ(ierr); 615 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 616 617 // Inverse of the lumped mass matrix 618 ierr = VecPointwiseMult(G, G, user->M); // M is Minv 619 CHKERRQ(ierr); 620 621 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 622 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 623 PetscFunctionReturn(0); 624 } 625 626 static PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Qdot, Vec G, 627 void *userData) { 628 PetscErrorCode ierr; 629 User user = *(User *)userData; 630 const PetscScalar *q, *qdot; 631 PetscScalar *g; 632 Vec Qloc, Qdotloc, Gloc; 633 634 // Global-to-local 635 PetscFunctionBeginUser; 636 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 637 ierr = DMGetLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 638 ierr = DMGetLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 639 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 640 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 641 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, Qloc, 0.0, 642 NULL, NULL, NULL); CHKERRQ(ierr); 643 ierr = VecZeroEntries(Qdotloc); CHKERRQ(ierr); 644 ierr = DMGlobalToLocal(user->dm, Qdot, INSERT_VALUES, Qdotloc); CHKERRQ(ierr); 645 ierr = VecZeroEntries(Gloc); CHKERRQ(ierr); 646 647 // Ceed Vectors 648 ierr = VecGetArrayRead(Qloc, &q); CHKERRQ(ierr); 649 ierr = VecGetArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 650 ierr = VecGetArray(Gloc, &g); CHKERRQ(ierr); 651 CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, 652 (PetscScalar *)q); 653 CeedVectorSetArray(user->qdotceed, CEED_MEM_HOST, CEED_USE_POINTER, 654 (PetscScalar *)qdot); 655 CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); 656 657 // Apply CEED operator 658 CeedOperatorApply(user->op_ifunction, user->qceed, user->gceed, 659 CEED_REQUEST_IMMEDIATE); 660 661 // Restore vectors 662 ierr = VecRestoreArrayRead(Qloc, &q); CHKERRQ(ierr); 663 ierr = VecRestoreArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); 664 ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); 665 666 ierr = VecZeroEntries(G); CHKERRQ(ierr); 667 ierr = DMLocalToGlobal(user->dm, Gloc, ADD_VALUES, G); CHKERRQ(ierr); 668 669 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 670 ierr = DMRestoreLocalVector(user->dm, &Qdotloc); CHKERRQ(ierr); 671 ierr = DMRestoreLocalVector(user->dm, &Gloc); CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 // User provided TS Monitor 676 static PetscErrorCode TSMonitor_NS(TS ts, PetscInt stepno, PetscReal time, 677 Vec Q, void *ctx) { 678 User user = ctx; 679 Vec Qloc; 680 char filepath[PETSC_MAX_PATH_LEN]; 681 PetscViewer viewer; 682 PetscErrorCode ierr; 683 684 // Set up output 685 PetscFunctionBeginUser; 686 // Print every 'outputfreq' steps 687 if (stepno % user->outputfreq != 0) 688 PetscFunctionReturn(0); 689 ierr = DMGetLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 690 ierr = PetscObjectSetName((PetscObject)Qloc, "StateVec"); CHKERRQ(ierr); 691 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 692 ierr = DMGlobalToLocal(user->dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 693 694 // Output 695 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-%03D.vtu", 696 user->outputdir, stepno + user->contsteps); 697 CHKERRQ(ierr); 698 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Q), filepath, 699 FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); 700 ierr = VecView(Qloc, viewer); CHKERRQ(ierr); 701 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 702 if (user->dmviz) { 703 Vec Qrefined, Qrefined_loc; 704 char filepath_refined[PETSC_MAX_PATH_LEN]; 705 PetscViewer viewer_refined; 706 707 ierr = DMGetGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 708 ierr = DMGetLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 709 ierr = PetscObjectSetName((PetscObject)Qrefined_loc, "Refined"); 710 CHKERRQ(ierr); 711 ierr = MatInterpolate(user->interpviz, Q, Qrefined); CHKERRQ(ierr); 712 ierr = VecZeroEntries(Qrefined_loc); CHKERRQ(ierr); 713 ierr = DMGlobalToLocal(user->dmviz, Qrefined, INSERT_VALUES, Qrefined_loc); 714 CHKERRQ(ierr); 715 ierr = PetscSNPrintf(filepath_refined, sizeof filepath_refined, 716 "%s/nsrefined-%03D.vtu", 717 user->outputdir, stepno + user->contsteps); 718 CHKERRQ(ierr); 719 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)Qrefined), 720 filepath_refined, 721 FILE_MODE_WRITE, &viewer_refined); CHKERRQ(ierr); 722 ierr = VecView(Qrefined_loc, viewer_refined); CHKERRQ(ierr); 723 ierr = DMRestoreLocalVector(user->dmviz, &Qrefined_loc); CHKERRQ(ierr); 724 ierr = DMRestoreGlobalVector(user->dmviz, &Qrefined); CHKERRQ(ierr); 725 ierr = PetscViewerDestroy(&viewer_refined); CHKERRQ(ierr); 726 } 727 ierr = DMRestoreLocalVector(user->dm, &Qloc); CHKERRQ(ierr); 728 729 // Save data in a binary file for continuation of simulations 730 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 731 user->outputdir); CHKERRQ(ierr); 732 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 733 CHKERRQ(ierr); 734 ierr = VecView(Q, viewer); CHKERRQ(ierr); 735 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 736 737 // Save time stamp 738 // Dimensionalize time back 739 time /= user->units->second; 740 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 741 user->outputdir); CHKERRQ(ierr); 742 ierr = PetscViewerBinaryOpen(user->comm, filepath, FILE_MODE_WRITE, &viewer); 743 CHKERRQ(ierr); 744 #if PETSC_VERSION_GE(3,13,0) 745 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL); 746 #else 747 ierr = PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL, true); 748 #endif 749 CHKERRQ(ierr); 750 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 751 752 PetscFunctionReturn(0); 753 } 754 755 static PetscErrorCode ICs_FixMultiplicity(CeedOperator op_ics, 756 CeedVector xcorners, CeedVector q0ceed, DM dm, Vec Qloc, Vec Q, 757 CeedElemRestriction restrictq, CeedQFunctionContext ctxSetup, CeedScalar time) { 758 PetscErrorCode ierr; 759 CeedVector multlvec; 760 Vec Multiplicity, MultiplicityLoc; 761 762 SetupContext ctxSetupData; 763 CeedQFunctionContextGetData(ctxSetup, CEED_MEM_HOST, (void **)&ctxSetupData); 764 ctxSetupData->time = time; 765 CeedQFunctionContextRestoreData(ctxSetup, (void **)&ctxSetupData); 766 767 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 768 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 769 CeedOperatorApply(op_ics, xcorners, q0ceed, CEED_REQUEST_IMMEDIATE); 770 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 771 ierr = DMLocalToGlobal(dm, Qloc, ADD_VALUES, Q); CHKERRQ(ierr); 772 773 // Fix multiplicity for output of ICs 774 ierr = DMGetLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 775 CeedElemRestrictionCreateVector(restrictq, &multlvec, NULL); 776 ierr = VectorPlacePetscVec(multlvec, MultiplicityLoc); CHKERRQ(ierr); 777 CeedElemRestrictionGetMultiplicity(restrictq, multlvec); 778 CeedVectorDestroy(&multlvec); 779 ierr = DMGetGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 780 ierr = VecZeroEntries(Multiplicity); CHKERRQ(ierr); 781 ierr = DMLocalToGlobal(dm, MultiplicityLoc, ADD_VALUES, Multiplicity); 782 CHKERRQ(ierr); 783 ierr = VecPointwiseDivide(Q, Q, Multiplicity); CHKERRQ(ierr); 784 ierr = VecPointwiseDivide(Qloc, Qloc, MultiplicityLoc); CHKERRQ(ierr); 785 ierr = DMRestoreLocalVector(dm, &MultiplicityLoc); CHKERRQ(ierr); 786 ierr = DMRestoreGlobalVector(dm, &Multiplicity); CHKERRQ(ierr); 787 788 PetscFunctionReturn(0); 789 } 790 791 static PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, 792 CeedElemRestriction restrictq, CeedBasis basisq, 793 CeedElemRestriction restrictqdi, CeedVector qdata, Vec M) { 794 PetscErrorCode ierr; 795 CeedQFunction qf_mass; 796 CeedOperator op_mass; 797 CeedVector mceed; 798 Vec Mloc; 799 CeedInt ncompq, qdatasize; 800 801 PetscFunctionBeginUser; 802 CeedElemRestrictionGetNumComponents(restrictq, &ncompq); 803 CeedElemRestrictionGetNumComponents(restrictqdi, &qdatasize); 804 // Create the Q-function that defines the action of the mass operator 805 CeedQFunctionCreateInterior(ceed, 1, Mass, Mass_loc, &qf_mass); 806 CeedQFunctionAddInput(qf_mass, "q", ncompq, CEED_EVAL_INTERP); 807 CeedQFunctionAddInput(qf_mass, "qdata", qdatasize, CEED_EVAL_NONE); 808 CeedQFunctionAddOutput(qf_mass, "v", ncompq, CEED_EVAL_INTERP); 809 810 // Create the mass operator 811 CeedOperatorCreate(ceed, qf_mass, NULL, NULL, &op_mass); 812 CeedOperatorSetField(op_mass, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 813 CeedOperatorSetField(op_mass, "qdata", restrictqdi, 814 CEED_BASIS_COLLOCATED, qdata); 815 CeedOperatorSetField(op_mass, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 816 817 ierr = DMGetLocalVector(dm, &Mloc); CHKERRQ(ierr); 818 ierr = VecZeroEntries(Mloc); CHKERRQ(ierr); 819 CeedElemRestrictionCreateVector(restrictq, &mceed, NULL); 820 ierr = VectorPlacePetscVec(mceed, Mloc); CHKERRQ(ierr); 821 822 { 823 // Compute a lumped mass matrix 824 CeedVector onesvec; 825 CeedElemRestrictionCreateVector(restrictq, &onesvec, NULL); 826 CeedVectorSetValue(onesvec, 1.0); 827 CeedOperatorApply(op_mass, onesvec, mceed, CEED_REQUEST_IMMEDIATE); 828 CeedVectorDestroy(&onesvec); 829 CeedOperatorDestroy(&op_mass); 830 CeedVectorDestroy(&mceed); 831 } 832 CeedQFunctionDestroy(&qf_mass); 833 834 ierr = VecZeroEntries(M); CHKERRQ(ierr); 835 ierr = DMLocalToGlobal(dm, Mloc, ADD_VALUES, M); CHKERRQ(ierr); 836 ierr = DMRestoreLocalVector(dm, &Mloc); CHKERRQ(ierr); 837 838 // Invert diagonally lumped mass vector for RHS function 839 ierr = VecReciprocal(M); CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 843 static PetscErrorCode SetUpDM(DM dm, problemData *problem, PetscInt degree, 844 SimpleBC bc, void *ctxSetupData) { 845 PetscErrorCode ierr; 846 847 PetscFunctionBeginUser; 848 { 849 // Configure the finite element space and boundary conditions 850 PetscFE fe; 851 PetscInt ncompq = 5; 852 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, ncompq, 853 PETSC_FALSE, degree, PETSC_DECIDE, 854 &fe); CHKERRQ(ierr); 855 ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); 856 ierr = DMAddField(dm, NULL,(PetscObject)fe); CHKERRQ(ierr); 857 ierr = DMCreateDS(dm); CHKERRQ(ierr); 858 { 859 PetscInt comps[1] = {1}; 860 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 861 1, comps, (void(*)(void))NULL, NULL, bc->nslip[0], 862 bc->slips[0], ctxSetupData); CHKERRQ(ierr); 863 comps[0] = 2; 864 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 865 1, comps, (void(*)(void))NULL, NULL, bc->nslip[1], 866 bc->slips[1], ctxSetupData); CHKERRQ(ierr); 867 comps[0] = 3; 868 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 869 1, comps, (void(*)(void))NULL, NULL, bc->nslip[2], 870 bc->slips[2], ctxSetupData); CHKERRQ(ierr); 871 } 872 if (bc->userbc == PETSC_TRUE) { 873 for (PetscInt c = 0; c < 3; c++) { 874 for (PetscInt s = 0; s < bc->nslip[c]; s++) { 875 for (PetscInt w = 0; w < bc->nwall; w++) { 876 if (bc->slips[c][s] == bc->walls[w]) 877 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, 878 "Boundary condition already set on face %D!\n", 879 bc->walls[w]); 880 881 } 882 } 883 } 884 } 885 // Wall boundary conditions are zero energy density and zero flux for 886 // velocity in advection/advection2d, and zero velocity and zero flux 887 // for mass density and energy density in density_current 888 { 889 if (problem->bc == Exact_Advection || problem->bc == Exact_Advection2d) { 890 PetscInt comps[1] = {4}; 891 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 892 1, comps, (void(*)(void))problem->bc, NULL, 893 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 894 } else if (problem->bc == Exact_DC || problem->bc == Exact_Euler) { 895 PetscInt comps[3] = {1, 2, 3}; 896 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 897 3, comps, (void(*)(void))problem->bc, NULL, 898 bc->nwall, bc->walls, ctxSetupData); CHKERRQ(ierr); 899 } else 900 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, 901 "Undefined boundary conditions for this problem"); 902 } 903 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 904 CHKERRQ(ierr); 905 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 906 } 907 { 908 // Empty name for conserved field (because there is only one field) 909 PetscSection section; 910 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 911 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 912 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 913 CHKERRQ(ierr); 914 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 915 CHKERRQ(ierr); 916 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 917 CHKERRQ(ierr); 918 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 919 CHKERRQ(ierr); 920 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 921 CHKERRQ(ierr); 922 } 923 PetscFunctionReturn(0); 924 } 925 926 int main(int argc, char **argv) { 927 PetscInt ierr; 928 MPI_Comm comm; 929 DM dm, dmcoord, dmviz; 930 Mat interpviz; 931 TS ts; 932 TSAdapt adapt; 933 User user; 934 Units units; 935 char ceedresource[4096] = "/cpu/self"; 936 PetscInt localNelemVol, lnodes, gnodes, steps; 937 const PetscInt ncompq = 5; 938 PetscMPIInt rank; 939 PetscScalar ftime; 940 Vec Q, Qloc, Xloc; 941 Ceed ceed; 942 CeedInt numP, numQ; 943 CeedVector xcorners, qdata, q0ceed; 944 CeedBasis basisx, basisxc, basisq; 945 CeedElemRestriction restrictx, restrictq, restrictqdi; 946 CeedQFunction qf_setupVol, qf_ics, qf_rhsVol, qf_ifunctionVol; 947 CeedQFunctionContext ctxSetup, ctxNS, ctxAdvection2d, ctxSurface; 948 CeedOperator op_setupVol, op_ics; 949 CeedScalar Rd; 950 CeedMemType memtyperequested; 951 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 952 kgpersquaredms, Joulepercubicm, Joule; 953 problemType problemChoice; 954 problemData *problem = NULL; 955 WindType wind_type; 956 StabilizationType stab; 957 PetscBool implicit; 958 PetscInt viz_refine = 0; 959 struct SimpleBC_ bc = { 960 .nslip = {2, 2, 2}, 961 .slips = {{5, 6}, {3, 4}, {1, 2}} 962 }; 963 double start, cpu_time_used; 964 // Test variables 965 PetscBool test; 966 PetscScalar testtol = 0.; 967 char filepath[PETSC_MAX_PATH_LEN]; 968 // Check PETSc CUDA support 969 PetscBool petschavecuda, setmemtyperequest = PETSC_FALSE; 970 // *INDENT-OFF* 971 #ifdef PETSC_HAVE_CUDA 972 petschavecuda = PETSC_TRUE; 973 #else 974 petschavecuda = PETSC_FALSE; 975 #endif 976 // *INDENT-ON* 977 978 // Create the libCEED contexts 979 PetscScalar meter = 1e-2; // 1 meter in scaled length units 980 PetscScalar second = 1e-2; // 1 second in scaled time units 981 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 982 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 983 CeedScalar theta0 = 300.; // K 984 CeedScalar thetaC = -15.; // K 985 CeedScalar P0 = 1.e5; // Pa 986 CeedScalar E_wind = 1.e6; // J 987 CeedScalar N = 0.01; // 1/s 988 CeedScalar cv = 717.; // J/(kg K) 989 CeedScalar cp = 1004.; // J/(kg K) 990 CeedScalar vortex_strength = 5.; // - 991 CeedScalar rho_enter = 1.2; // Kg/m^3 992 CeedScalar g = 9.81; // m/s^2 993 CeedScalar lambda = -2./3.; // - 994 CeedScalar mu = 75.; // Pa s, dynamic viscosity 995 // mu = 75 is not physical for air, but is good for numerical stability 996 CeedScalar k = 0.02638; // W/(m K) 997 CeedScalar CtauS = 0.; // dimensionless 998 CeedScalar strong_form = 0.; // [0,1] 999 PetscScalar lx = 8000.; // m 1000 PetscScalar ly = 8000.; // m 1001 PetscScalar lz = 4000.; // m 1002 CeedScalar rc = 1000.; // m (Radius of bubble) 1003 PetscScalar resx = 1000.; // m (resolution in x) 1004 PetscScalar resy = 1000.; // m (resolution in y) 1005 PetscScalar resz = 1000.; // m (resolution in z) 1006 PetscInt outputfreq = 10; // - 1007 PetscInt contsteps = 0; // - 1008 PetscInt degree = 1; // - 1009 PetscInt qextra = 2; // - 1010 PetscInt qextraSur = 2; // - 1011 PetscReal center[3], dc_axis[3] = {0, 0, 0}, wind[3] = {1., 0, 0}, 1012 u_enter[3] = {2., 0, 0}; 1013 1014 ierr = PetscInitialize(&argc, &argv, NULL, help); 1015 if (ierr) return ierr; 1016 1017 // Allocate PETSc context 1018 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 1019 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 1020 1021 // Parse command line options 1022 comm = PETSC_COMM_WORLD; 1023 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 1024 NULL); CHKERRQ(ierr); 1025 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 1026 NULL, ceedresource, ceedresource, 1027 sizeof(ceedresource), NULL); CHKERRQ(ierr); 1028 ierr = PetscOptionsBool("-test", "Run in test mode", 1029 NULL, test=PETSC_FALSE, &test, NULL); CHKERRQ(ierr); 1030 ierr = PetscOptionsScalar("-compare_final_state_atol", 1031 "Test absolute tolerance", 1032 NULL, testtol, &testtol, NULL); CHKERRQ(ierr); 1033 ierr = PetscOptionsString("-compare_final_state_filename", "Test filename", 1034 NULL, filepath, filepath, 1035 sizeof(filepath), NULL); CHKERRQ(ierr); 1036 problemChoice = NS_DENSITY_CURRENT; 1037 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 1038 problemTypes, (PetscEnum)problemChoice, 1039 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 1040 problem = &problemOptions[problemChoice]; 1041 ierr = PetscOptionsEnum("-problem_advection_wind", "Wind type in Advection", 1042 NULL, WindTypes, 1043 (PetscEnum)(wind_type = ADVECTION_WIND_ROTATION), 1044 (PetscEnum *)&wind_type, NULL); CHKERRQ(ierr); 1045 PetscInt n = problem->dim; 1046 PetscBool userWind; 1047 ierr = PetscOptionsRealArray("-problem_advection_wind_translation", 1048 "Constant wind vector", 1049 NULL, wind, &n, &userWind); CHKERRQ(ierr); 1050 if (wind_type == ADVECTION_WIND_ROTATION && userWind) { 1051 ierr = PetscPrintf(comm, 1052 "Warning! Use -problem_advection_wind_translation only with -problem_advection_wind translation\n"); 1053 CHKERRQ(ierr); 1054 } 1055 if (wind_type == ADVECTION_WIND_TRANSLATION 1056 && (problemChoice == NS_DENSITY_CURRENT || 1057 problemChoice == NS_EULER_VORTEX)) { 1058 SETERRQ(comm, PETSC_ERR_ARG_INCOMP, 1059 "-problem_advection_wind translation is not defined for -problem density_current or -problem euler_vortex"); 1060 } 1061 ierr = PetscOptionsRealArray("-problem_euler_vortex_velocity", 1062 "Incoming velocity vector", 1063 NULL, u_enter, &n, NULL); CHKERRQ(ierr); 1064 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 1065 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 1066 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 1067 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 1068 NULL, implicit=PETSC_FALSE, &implicit, NULL); 1069 CHKERRQ(ierr); 1070 if (!implicit && stab != STAB_NONE) { 1071 ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); 1072 CHKERRQ(ierr); 1073 } 1074 { 1075 PetscInt len; 1076 PetscBool flg; 1077 ierr = PetscOptionsIntArray("-bc_wall", 1078 "Use wall boundary conditions on this list of faces", 1079 NULL, bc.walls, 1080 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 1081 &len), &flg); CHKERRQ(ierr); 1082 if (flg) { 1083 bc.nwall = len; 1084 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 1085 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1086 } 1087 for (PetscInt j=0; j<3; j++) { 1088 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 1089 ierr = PetscOptionsIntArray(flags[j], 1090 "Use slip boundary conditions on this list of faces", 1091 NULL, bc.slips[j], 1092 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 1093 &len), &flg); 1094 CHKERRQ(ierr); 1095 if (flg) { 1096 bc.nslip[j] = len; 1097 bc.userbc = PETSC_TRUE; 1098 } 1099 } 1100 } 1101 ierr = PetscOptionsInt("-viz_refine", 1102 "Regular refinement levels for visualization", 1103 NULL, viz_refine, &viz_refine, NULL); 1104 CHKERRQ(ierr); 1105 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1106 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1107 meter = fabs(meter); 1108 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1109 NULL, second, &second, NULL); CHKERRQ(ierr); 1110 second = fabs(second); 1111 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1112 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1113 kilogram = fabs(kilogram); 1114 ierr = PetscOptionsScalar("-units_Kelvin", 1115 "1 Kelvin in scaled temperature units", 1116 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1117 Kelvin = fabs(Kelvin); 1118 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1119 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1120 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1121 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1122 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1123 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1124 ierr = PetscOptionsScalar("-E_wind", "Total energy of inflow wind", 1125 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 1126 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1127 NULL, N, &N, NULL); CHKERRQ(ierr); 1128 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1129 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1130 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1131 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1132 PetscBool userVortex; 1133 ierr = PetscOptionsScalar("-vortex_strength", "Strength of Vortex", 1134 NULL, vortex_strength, &vortex_strength, &userVortex); 1135 CHKERRQ(ierr); 1136 if (problemChoice != NS_EULER_VORTEX && userVortex) { 1137 ierr = PetscPrintf(comm, 1138 "Warning! Use -vortex_strength only with -problem euler_vortex\n"); 1139 CHKERRQ(ierr); 1140 } 1141 ierr = PetscOptionsScalar("-problem_euler_vortex_rho", "Incoming density", 1142 NULL, rho_enter, &rho_enter, NULL); 1143 CHKERRQ(ierr); 1144 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1145 NULL, g, &g, NULL); CHKERRQ(ierr); 1146 ierr = PetscOptionsScalar("-lambda", 1147 "Stokes hypothesis second viscosity coefficient", 1148 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1149 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1150 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1151 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1152 NULL, k, &k, NULL); CHKERRQ(ierr); 1153 ierr = PetscOptionsScalar("-CtauS", 1154 "Scale coefficient for tau (nondimensional)", 1155 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1156 if (stab == STAB_NONE && CtauS != 0) { 1157 ierr = PetscPrintf(comm, 1158 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1159 CHKERRQ(ierr); 1160 } 1161 ierr = PetscOptionsScalar("-strong_form", 1162 "Strong (1) or weak/integrated by parts (0) advection residual", 1163 NULL, strong_form, &strong_form, NULL); 1164 CHKERRQ(ierr); 1165 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1166 ierr = PetscPrintf(comm, 1167 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1168 CHKERRQ(ierr); 1169 } 1170 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1171 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1172 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1173 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1174 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1175 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1176 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1177 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1178 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1179 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1180 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1181 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1182 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1183 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1184 n = problem->dim; 1185 center[0] = 0.5 * lx; 1186 center[1] = 0.5 * ly; 1187 center[2] = 0.5 * lz; 1188 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1189 NULL, center, &n, NULL); CHKERRQ(ierr); 1190 n = problem->dim; 1191 ierr = PetscOptionsRealArray("-dc_axis", 1192 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1193 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1194 { 1195 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1196 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1197 if (norm > 0) { 1198 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1199 } 1200 } 1201 ierr = PetscOptionsInt("-output_freq", 1202 "Frequency of output, in number of steps", 1203 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1204 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1205 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1206 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1207 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1208 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1209 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1210 PetscBool userQextraSur; 1211 ierr = PetscOptionsInt("-qextra_boundary", 1212 "Number of extra quadrature points on in/outflow faces", 1213 NULL, qextraSur, &qextraSur, &userQextraSur); 1214 CHKERRQ(ierr); 1215 if ((wind_type == ADVECTION_WIND_ROTATION 1216 || problemChoice == NS_DENSITY_CURRENT) && userQextraSur) { 1217 ierr = PetscPrintf(comm, 1218 "Warning! Use -qextra_boundary only with -problem_advection_wind translation\n"); 1219 CHKERRQ(ierr); 1220 } 1221 ierr = PetscStrncpy(user->outputdir, ".", 2); CHKERRQ(ierr); 1222 ierr = PetscOptionsString("-output_dir", "Output directory", 1223 NULL, user->outputdir, user->outputdir, 1224 sizeof(user->outputdir), NULL); CHKERRQ(ierr); 1225 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1226 ierr = PetscOptionsEnum("-memtype", 1227 "CEED MemType requested", NULL, 1228 memTypes, (PetscEnum)memtyperequested, 1229 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1230 CHKERRQ(ierr); 1231 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1232 1233 // Define derived units 1234 Pascal = kilogram / (meter * PetscSqr(second)); 1235 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1236 mpersquareds = meter / PetscSqr(second); 1237 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1238 kgpercubicm = kilogram / pow(meter,3); 1239 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1240 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1241 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 1242 1243 // Scale variables to desired units 1244 theta0 *= Kelvin; 1245 thetaC *= Kelvin; 1246 P0 *= Pascal; 1247 E_wind *= Joule; 1248 rho_enter *= kgpercubicm; 1249 N *= (1./second); 1250 cv *= JperkgK; 1251 cp *= JperkgK; 1252 Rd = cp - cv; 1253 g *= mpersquareds; 1254 mu *= Pascal * second; 1255 k *= WpermK; 1256 lx = fabs(lx) * meter; 1257 ly = fabs(ly) * meter; 1258 lz = fabs(lz) * meter; 1259 rc = fabs(rc) * meter; 1260 resx = fabs(resx) * meter; 1261 resy = fabs(resy) * meter; 1262 resz = fabs(resz) * meter; 1263 for (int i=0; i<3; i++) center[i] *= meter; 1264 1265 const CeedInt dim = problem->dim, ncompx = problem->dim, 1266 qdatasizeVol = problem->qdatasizeVol; 1267 // Set up the libCEED context 1268 struct SetupContext_ ctxSetupData = { 1269 .theta0 = theta0, 1270 .thetaC = thetaC, 1271 .P0 = P0, 1272 .N = N, 1273 .cv = cv, 1274 .cp = cp, 1275 .Rd = Rd, 1276 .g = g, 1277 .rc = rc, 1278 .lx = lx, 1279 .ly = ly, 1280 .lz = lz, 1281 .center[0] = center[0], 1282 .center[1] = center[1], 1283 .center[2] = center[2], 1284 .dc_axis[0] = dc_axis[0], 1285 .dc_axis[1] = dc_axis[1], 1286 .dc_axis[2] = dc_axis[2], 1287 .wind[0] = wind[0], 1288 .wind[1] = wind[1], 1289 .wind[2] = wind[2], 1290 .time = 0, 1291 .vortex_strength = vortex_strength, 1292 .wind_type = wind_type, 1293 }; 1294 1295 // Create the mesh 1296 { 1297 const PetscReal scale[3] = {lx, ly, lz}; 1298 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1299 NULL, PETSC_TRUE, &dm); 1300 CHKERRQ(ierr); 1301 } 1302 1303 // Distribute the mesh over processes 1304 { 1305 DM dmDist = NULL; 1306 PetscPartitioner part; 1307 1308 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1309 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1310 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1311 if (dmDist) { 1312 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1313 dm = dmDist; 1314 } 1315 } 1316 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1317 1318 // Setup DM 1319 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1320 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1321 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetupData); 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, qf_applyIn, qf_applyOut; 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 qf_applySur = NULL; 1567 if (problem->applySur) { 1568 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1569 problem->applySur_loc, &qf_applySur); 1570 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1571 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1572 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1573 CeedQFunctionAddOutput(qf_applySur, "v", ncompq, CEED_EVAL_INTERP); 1574 } 1575 qf_applyIn = NULL; 1576 if (problem->applyIn) { 1577 CeedQFunctionCreateInterior(ceed, 1, problem->applyIn, 1578 problem->applyIn_loc, &qf_applyIn); 1579 CeedQFunctionAddInput(qf_applyIn, "q", ncompq, CEED_EVAL_INTERP); 1580 CeedQFunctionAddInput(qf_applyIn, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1581 CeedQFunctionAddInput(qf_applyIn, "x", ncompx, CEED_EVAL_INTERP); 1582 CeedQFunctionAddOutput(qf_applyIn, "v", ncompq, CEED_EVAL_INTERP); 1583 } 1584 qf_applyOut = NULL; 1585 if (problem->applyOut) { 1586 CeedQFunctionCreateInterior(ceed, 1, problem->applyOut, 1587 problem->applyOut_loc, &qf_applyOut); 1588 CeedQFunctionAddInput(qf_applyOut, "q", ncompq, CEED_EVAL_INTERP); 1589 CeedQFunctionAddInput(qf_applyOut, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1590 CeedQFunctionAddInput(qf_applyOut, "x", ncompx, CEED_EVAL_INTERP); 1591 CeedQFunctionAddOutput(qf_applyOut, "v", ncompq, CEED_EVAL_INTERP); 1592 } 1593 1594 // Create CEED Operator for the whole domain 1595 if (!implicit) 1596 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1597 user->op_rhs_vol,qf_applySur, 1598 qf_setupSur, qf_applyIn, qf_applyOut, 1599 height, numP_Sur, numQ_Sur, 1600 qdatasizeSur, NqptsSur, basisxSur, 1601 basisqSur, &user->op_rhs); 1602 CHKERRQ(ierr); 1603 if (implicit) 1604 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1605 user->op_ifunction_vol, qf_applySur, 1606 qf_setupSur, NULL, NULL, 1607 height, numP_Sur, numQ_Sur, 1608 qdatasizeSur, NqptsSur, basisxSur, 1609 basisqSur, &user->op_ifunction); 1610 CHKERRQ(ierr); 1611 // Set up contex for QFunctions 1612 CeedQFunctionContextCreate(ceed, &ctxSetup); 1613 CeedQFunctionContextSetData(ctxSetup, CEED_MEM_HOST, CEED_USE_POINTER, 1614 sizeof ctxSetupData, &ctxSetupData); 1615 CeedQFunctionSetContext(qf_ics, ctxSetup); 1616 1617 CeedScalar ctxNSData[8] = {lambda, mu, k, cv, cp, g, Rd}; 1618 CeedQFunctionContextCreate(ceed, &ctxNS); 1619 CeedQFunctionContextSetData(ctxNS, CEED_MEM_HOST, CEED_USE_POINTER, 1620 sizeof ctxNSData, &ctxNSData); 1621 1622 struct Advection2dContext_ ctxAdvection2dData = { 1623 .CtauS = CtauS, 1624 .strong_form = strong_form, 1625 .stabilization = stab, 1626 }; 1627 CeedQFunctionContextCreate(ceed, &ctxAdvection2d); 1628 CeedQFunctionContextSetData(ctxAdvection2d, CEED_MEM_HOST, CEED_USE_POINTER, 1629 sizeof ctxAdvection2dData, &ctxAdvection2dData); 1630 1631 struct SurfaceContext_ ctxSurfaceData = { 1632 .E_wind = E_wind, 1633 .strong_form = strong_form, 1634 .implicit = implicit, 1635 }; 1636 struct EulerContext_ ctxEuler = { 1637 .cv = cv, 1638 .cp = cp, 1639 .rho_enter = rho_enter, 1640 .u_enter[0] = u_enter[0], 1641 .u_enter[1] = u_enter[1], 1642 .u_enter[2] = u_enter[2], 1643 .implicit = implicit, 1644 }; 1645 CeedQFunctionContextCreate(ceed, &ctxSurface); 1646 CeedQFunctionContextSetData(ctxSurface, CEED_MEM_HOST, CEED_USE_POINTER, 1647 sizeof ctxSurfaceData, &ctxSurfaceData); 1648 1649 switch (problemChoice) { 1650 case NS_DENSITY_CURRENT: 1651 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1652 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1653 sizeof ctxNS); 1654 break; 1655 case NS_ADVECTION: 1656 case NS_ADVECTION2D: 1657 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1658 sizeof ctxAdvection2d); 1659 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1660 sizeof ctxAdvection2d); 1661 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &ctxSurface, 1662 sizeof ctxSurface); 1663 case NS_EULER_VORTEX: 1664 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxEuler, 1665 sizeof ctxEuler); 1666 if (qf_applyOut) CeedQFunctionSetContext(qf_applyOut, &ctxEuler, 1667 sizeof ctxEuler); 1668 if (qf_applyIn) CeedQFunctionSetContext(qf_applyIn, &ctxEuler, 1669 sizeof ctxEuler); 1670 } 1671 1672 // Set up PETSc context 1673 // Set up units structure 1674 units->meter = meter; 1675 units->kilogram = kilogram; 1676 units->second = second; 1677 units->Kelvin = Kelvin; 1678 units->Pascal = Pascal; 1679 units->JperkgK = JperkgK; 1680 units->mpersquareds = mpersquareds; 1681 units->WpermK = WpermK; 1682 units->kgpercubicm = kgpercubicm; 1683 units->kgpersquaredms = kgpersquaredms; 1684 units->Joulepercubicm = Joulepercubicm; 1685 units->Joule = Joule; 1686 1687 // Set up user structure 1688 user->comm = comm; 1689 user->outputfreq = outputfreq; 1690 user->contsteps = contsteps; 1691 user->units = units; 1692 user->dm = dm; 1693 user->dmviz = dmviz; 1694 user->interpviz = interpviz; 1695 user->ceed = ceed; 1696 1697 // Calculate qdata and ICs 1698 // Set up state global and local vectors 1699 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1700 1701 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1702 1703 // Apply Setup Ceed Operators 1704 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1705 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1706 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1707 user->M); CHKERRQ(ierr); 1708 1709 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1710 ctxSetup, 0.0); CHKERRQ(ierr); 1711 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1712 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1713 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1714 // slower execution. 1715 Vec Qbc; 1716 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1717 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1718 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1719 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1720 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1721 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1722 ierr = PetscObjectComposeFunction((PetscObject)dm, 1723 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1724 CHKERRQ(ierr); 1725 } 1726 1727 MPI_Comm_rank(comm, &rank); 1728 if (!rank) {ierr = PetscMkdir(user->outputdir); CHKERRQ(ierr);} 1729 // Gather initial Q values 1730 // In case of continuation of simulation, set up initial values from binary file 1731 if (contsteps) { // continue from existent solution 1732 PetscViewer viewer; 1733 char filepath[PETSC_MAX_PATH_LEN]; 1734 // Read input 1735 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1736 user->outputdir); 1737 CHKERRQ(ierr); 1738 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1739 CHKERRQ(ierr); 1740 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1741 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1742 } 1743 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1744 1745 // Create and setup TS 1746 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1747 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1748 if (implicit) { 1749 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1750 if (user->op_ifunction) { 1751 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1752 } else { // Implicit integrators can fall back to using an RHSFunction 1753 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1754 } 1755 } else { 1756 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1757 "Problem does not provide RHSFunction"); 1758 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1759 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1760 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1761 } 1762 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1763 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1764 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1765 if (test) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1766 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1767 ierr = TSAdaptSetStepLimits(adapt, 1768 1.e-12 * units->second, 1769 1.e2 * units->second); CHKERRQ(ierr); 1770 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1771 if (!contsteps) { // print initial condition 1772 if (!test) { 1773 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1774 } 1775 } else { // continue from time of last output 1776 PetscReal time; 1777 PetscInt count; 1778 PetscViewer viewer; 1779 char filepath[PETSC_MAX_PATH_LEN]; 1780 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1781 user->outputdir); CHKERRQ(ierr); 1782 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1783 CHKERRQ(ierr); 1784 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1785 CHKERRQ(ierr); 1786 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1787 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1788 } 1789 if (!test) { 1790 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1791 } 1792 1793 // Solve 1794 start = MPI_Wtime(); 1795 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1796 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1797 cpu_time_used = MPI_Wtime() - start; 1798 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1799 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1800 comm); CHKERRQ(ierr); 1801 if (!test) { 1802 ierr = PetscPrintf(PETSC_COMM_WORLD, 1803 "Time taken for solution (sec): %g\n", 1804 (double)cpu_time_used); CHKERRQ(ierr); 1805 } 1806 1807 // Get error 1808 if (problem->non_zero_time && !test) { 1809 Vec Qexact, Qexactloc; 1810 PetscReal norm; 1811 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1812 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1813 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1814 1815 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1816 restrictq, ctxSetup, ftime); CHKERRQ(ierr); 1817 1818 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1819 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1820 CeedVectorDestroy(&q0ceed); 1821 ierr = PetscPrintf(PETSC_COMM_WORLD, 1822 "Max Error: %g\n", 1823 (double)norm); CHKERRQ(ierr); 1824 // Clean up vectors 1825 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1826 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1827 } 1828 1829 // Output Statistics 1830 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1831 if (!test) { 1832 ierr = PetscPrintf(PETSC_COMM_WORLD, 1833 "Time integrator took %D time steps to reach final time %g\n", 1834 steps, (double)ftime); CHKERRQ(ierr); 1835 } 1836 // Output numerical values from command line 1837 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1838 1839 // Compare reference solution values with current test run for CI 1840 if (test) { 1841 PetscViewer viewer; 1842 // Read reference file 1843 Vec Qref; 1844 PetscReal error, Qrefnorm; 1845 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1846 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1847 CHKERRQ(ierr); 1848 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1849 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1850 1851 // Compute error with respect to reference solution 1852 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1853 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1854 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1855 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1856 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1857 // Check error 1858 if (error > testtol) { 1859 ierr = PetscPrintf(PETSC_COMM_WORLD, 1860 "Test failed with error norm %g\n", 1861 (double)error); CHKERRQ(ierr); 1862 } 1863 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1864 } 1865 1866 // Clean up libCEED 1867 CeedVectorDestroy(&qdata); 1868 CeedVectorDestroy(&user->qceed); 1869 CeedVectorDestroy(&user->qdotceed); 1870 CeedVectorDestroy(&user->gceed); 1871 CeedVectorDestroy(&xcorners); 1872 CeedBasisDestroy(&basisq); 1873 CeedBasisDestroy(&basisx); 1874 CeedBasisDestroy(&basisxc); 1875 CeedElemRestrictionDestroy(&restrictq); 1876 CeedElemRestrictionDestroy(&restrictx); 1877 CeedElemRestrictionDestroy(&restrictqdi); 1878 CeedQFunctionDestroy(&qf_setupVol); 1879 CeedQFunctionDestroy(&qf_ics); 1880 CeedQFunctionDestroy(&qf_rhsVol); 1881 CeedQFunctionDestroy(&qf_ifunctionVol); 1882 CeedQFunctionContextDestroy(&ctxSetup); 1883 CeedQFunctionContextDestroy(&ctxNS); 1884 CeedQFunctionContextDestroy(&ctxAdvection2d); 1885 CeedQFunctionContextDestroy(&ctxSurface); 1886 CeedOperatorDestroy(&op_setupVol); 1887 CeedOperatorDestroy(&op_ics); 1888 CeedOperatorDestroy(&user->op_rhs_vol); 1889 CeedOperatorDestroy(&user->op_ifunction_vol); 1890 CeedDestroy(&ceed); 1891 CeedBasisDestroy(&basisqSur); 1892 CeedBasisDestroy(&basisxSur); 1893 CeedBasisDestroy(&basisxcSur); 1894 CeedQFunctionDestroy(&qf_setupSur); 1895 CeedQFunctionDestroy(&qf_applySur); 1896 CeedOperatorDestroy(&user->op_rhs); 1897 CeedOperatorDestroy(&user->op_ifunction); 1898 1899 // Clean up PETSc 1900 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1901 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1902 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1903 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1904 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1905 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1906 ierr = PetscFree(units); CHKERRQ(ierr); 1907 ierr = PetscFree(user); CHKERRQ(ierr); 1908 return PetscFinalize(); 1909 } 1910