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