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