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