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