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 { 1057 PetscInt len; 1058 PetscBool flg; 1059 ierr = PetscOptionsIntArray("-bc_wall", 1060 "Use wall boundary conditions on this list of faces", 1061 NULL, bc.walls, 1062 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 1063 &len), &flg); CHKERRQ(ierr); 1064 if (flg) { 1065 bc.nwall = len; 1066 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 1067 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1068 } 1069 for (PetscInt j=0; j<3; j++) { 1070 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 1071 ierr = PetscOptionsIntArray(flags[j], 1072 "Use slip boundary conditions on this list of faces", 1073 NULL, bc.slips[j], 1074 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 1075 &len), &flg); 1076 CHKERRQ(ierr); 1077 if (flg) { 1078 bc.nslip[j] = len; 1079 bc.userbc = PETSC_TRUE; 1080 } 1081 } 1082 } 1083 ierr = PetscOptionsInt("-viz_refine", 1084 "Regular refinement levels for visualization", 1085 NULL, viz_refine, &viz_refine, NULL); 1086 CHKERRQ(ierr); 1087 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1088 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1089 meter = fabs(meter); 1090 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1091 NULL, second, &second, NULL); CHKERRQ(ierr); 1092 second = fabs(second); 1093 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1094 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1095 kilogram = fabs(kilogram); 1096 ierr = PetscOptionsScalar("-units_Kelvin", 1097 "1 Kelvin in scaled temperature units", 1098 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1099 Kelvin = fabs(Kelvin); 1100 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1101 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1102 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1103 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1104 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1105 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1106 ierr = PetscOptionsScalar("-E_wind", "Total energy of inflow wind", 1107 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 1108 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1109 NULL, N, &N, NULL); CHKERRQ(ierr); 1110 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1111 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1112 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1113 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1114 PetscBool userVortex; 1115 ierr = PetscOptionsScalar("-vortex_strength", "Strength of Vortex", 1116 NULL, vortex_strength, &vortex_strength, &userVortex); 1117 CHKERRQ(ierr); 1118 if (problemChoice != NS_EULER_VORTEX && userVortex) { 1119 ierr = PetscPrintf(comm, 1120 "Warning! Use -vortex_strength only with -problem euler_vortex\n"); 1121 CHKERRQ(ierr); 1122 } 1123 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1124 NULL, g, &g, NULL); CHKERRQ(ierr); 1125 ierr = PetscOptionsScalar("-lambda", 1126 "Stokes hypothesis second viscosity coefficient", 1127 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1128 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1129 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1130 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1131 NULL, k, &k, NULL); CHKERRQ(ierr); 1132 ierr = PetscOptionsScalar("-CtauS", 1133 "Scale coefficient for tau (nondimensional)", 1134 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1135 if (stab == STAB_NONE && CtauS != 0) { 1136 ierr = PetscPrintf(comm, 1137 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1138 CHKERRQ(ierr); 1139 } 1140 ierr = PetscOptionsScalar("-strong_form", 1141 "Strong (1) or weak/integrated by parts (0) advection residual", 1142 NULL, strong_form, &strong_form, NULL); 1143 CHKERRQ(ierr); 1144 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1145 ierr = PetscPrintf(comm, 1146 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1147 CHKERRQ(ierr); 1148 } 1149 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1150 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1151 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1152 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1153 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1154 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1155 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1156 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1157 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1158 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1159 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1160 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1161 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1162 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1163 n = problem->dim; 1164 center[0] = 0.5 * lx; 1165 center[1] = 0.5 * ly; 1166 center[2] = 0.5 * lz; 1167 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1168 NULL, center, &n, NULL); CHKERRQ(ierr); 1169 n = problem->dim; 1170 ierr = PetscOptionsRealArray("-dc_axis", 1171 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1172 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1173 { 1174 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1175 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1176 if (norm > 0) { 1177 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1178 } 1179 } 1180 ierr = PetscOptionsInt("-output_freq", 1181 "Frequency of output, in number of steps", 1182 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1183 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1184 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1185 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1186 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1187 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1188 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1189 PetscBool userQextraSur; 1190 ierr = PetscOptionsInt("-qextra_boundary", 1191 "Number of extra quadrature points on in/outflow faces", 1192 NULL, qextraSur, &qextraSur, &userQextraSur); 1193 CHKERRQ(ierr); 1194 if ((wind_type == ADVECTION_WIND_ROTATION 1195 || problemChoice == NS_DENSITY_CURRENT) && userQextraSur) { 1196 ierr = PetscPrintf(comm, 1197 "Warning! Use -qextra_boundary only with -problem_advection_wind translation\n"); 1198 CHKERRQ(ierr); 1199 } 1200 ierr = PetscStrncpy(user->outputdir, ".", 2); CHKERRQ(ierr); 1201 ierr = PetscOptionsString("-output_dir", "Output directory", 1202 NULL, user->outputdir, user->outputdir, 1203 sizeof(user->outputdir), NULL); CHKERRQ(ierr); 1204 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1205 ierr = PetscOptionsEnum("-memtype", 1206 "CEED MemType requested", NULL, 1207 memTypes, (PetscEnum)memtyperequested, 1208 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1209 CHKERRQ(ierr); 1210 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1211 1212 // Define derived units 1213 Pascal = kilogram / (meter * PetscSqr(second)); 1214 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1215 mpersquareds = meter / PetscSqr(second); 1216 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1217 kgpercubicm = kilogram / pow(meter,3); 1218 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1219 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1220 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 1221 1222 // Scale variables to desired units 1223 theta0 *= Kelvin; 1224 thetaC *= Kelvin; 1225 P0 *= Pascal; 1226 E_wind *= Joule; 1227 N *= (1./second); 1228 cv *= JperkgK; 1229 cp *= JperkgK; 1230 Rd = cp - cv; 1231 g *= mpersquareds; 1232 mu *= Pascal * second; 1233 k *= WpermK; 1234 lx = fabs(lx) * meter; 1235 ly = fabs(ly) * meter; 1236 lz = fabs(lz) * meter; 1237 rc = fabs(rc) * meter; 1238 resx = fabs(resx) * meter; 1239 resy = fabs(resy) * meter; 1240 resz = fabs(resz) * meter; 1241 for (int i=0; i<3; i++) center[i] *= meter; 1242 1243 const CeedInt dim = problem->dim, ncompx = problem->dim, 1244 qdatasizeVol = problem->qdatasizeVol; 1245 // Set up the libCEED context 1246 struct SetupContext_ ctxSetupData = { 1247 .theta0 = theta0, 1248 .thetaC = thetaC, 1249 .P0 = P0, 1250 .N = N, 1251 .cv = cv, 1252 .cp = cp, 1253 .Rd = Rd, 1254 .g = g, 1255 .rc = rc, 1256 .lx = lx, 1257 .ly = ly, 1258 .lz = lz, 1259 .center[0] = center[0], 1260 .center[1] = center[1], 1261 .center[2] = center[2], 1262 .dc_axis[0] = dc_axis[0], 1263 .dc_axis[1] = dc_axis[1], 1264 .dc_axis[2] = dc_axis[2], 1265 .wind[0] = wind[0], 1266 .wind[1] = wind[1], 1267 .wind[2] = wind[2], 1268 .time = 0, 1269 .wind_type = wind_type, 1270 }; 1271 1272 // Create the mesh 1273 { 1274 const PetscReal scale[3] = {lx, ly, lz}; 1275 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1276 NULL, PETSC_TRUE, &dm); 1277 CHKERRQ(ierr); 1278 } 1279 1280 // Distribute the mesh over processes 1281 { 1282 DM dmDist = NULL; 1283 PetscPartitioner part; 1284 1285 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1286 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1287 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1288 if (dmDist) { 1289 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1290 dm = dmDist; 1291 } 1292 } 1293 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1294 1295 // Setup DM 1296 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1297 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1298 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetupData); 1299 CHKERRQ(ierr); 1300 1301 // Refine DM for high-order viz 1302 dmviz = NULL; 1303 interpviz = NULL; 1304 if (viz_refine) { 1305 DM dmhierarchy[viz_refine+1]; 1306 1307 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1308 dmhierarchy[0] = dm; 1309 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1310 Mat interp_next; 1311 1312 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1313 CHKERRQ(ierr); 1314 ierr = DMClearDS(dmhierarchy[i+1]); CHKERRQ(ierr); 1315 ierr = DMClearFields(dmhierarchy[i+1]); CHKERRQ(ierr); 1316 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1317 d = (d + 1) / 2; 1318 if (i + 1 == viz_refine) d = 1; 1319 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetupData); 1320 CHKERRQ(ierr); 1321 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1322 &interp_next, NULL); CHKERRQ(ierr); 1323 if (!i) interpviz = interp_next; 1324 else { 1325 Mat C; 1326 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1327 PETSC_DECIDE, &C); CHKERRQ(ierr); 1328 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1329 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1330 interpviz = C; 1331 } 1332 } 1333 for (PetscInt i=1; i<viz_refine; i++) { 1334 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1335 } 1336 dmviz = dmhierarchy[viz_refine]; 1337 } 1338 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1339 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1340 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1341 lnodes /= ncompq; 1342 1343 // Initialize CEED 1344 CeedInit(ceedresource, &ceed); 1345 // Set memtype 1346 CeedMemType memtypebackend; 1347 CeedGetPreferredMemType(ceed, &memtypebackend); 1348 // Check memtype compatibility 1349 if (!setmemtyperequest) 1350 memtyperequested = memtypebackend; 1351 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1352 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1353 "PETSc was not built with CUDA. " 1354 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1355 1356 // Set number of 1D nodes and quadrature points 1357 numP = degree + 1; 1358 numQ = numP + qextra; 1359 1360 // Print summary 1361 if (!test) { 1362 CeedInt gdofs, odofs; 1363 int comm_size; 1364 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1365 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1366 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1367 gnodes = gdofs/ncompq; 1368 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1369 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1370 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1371 const char *usedresource; 1372 CeedGetResource(ceed, &usedresource); 1373 1374 ierr = PetscPrintf(comm, 1375 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1376 " rank(s) : %d\n" 1377 " Problem:\n" 1378 " Problem Name : %s\n" 1379 " Stabilization : %s\n" 1380 " PETSc:\n" 1381 " Box Faces : %s\n" 1382 " libCEED:\n" 1383 " libCEED Backend : %s\n" 1384 " libCEED Backend MemType : %s\n" 1385 " libCEED User Requested MemType : %s\n" 1386 " Mesh:\n" 1387 " Number of 1D Basis Nodes (P) : %d\n" 1388 " Number of 1D Quadrature Points (Q) : %d\n" 1389 " Global DoFs : %D\n" 1390 " Owned DoFs : %D\n" 1391 " DoFs per node : %D\n" 1392 " Global nodes : %D\n" 1393 " Owned nodes : %D\n", 1394 comm_size, problemTypes[problemChoice], 1395 StabilizationTypes[stab], box_faces_str, usedresource, 1396 CeedMemTypes[memtypebackend], 1397 (setmemtyperequest) ? 1398 CeedMemTypes[memtyperequested] : "none", 1399 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1400 CHKERRQ(ierr); 1401 } 1402 1403 // Set up global mass vector 1404 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1405 1406 // Set up libCEED 1407 // CEED Bases 1408 CeedInit(ceedresource, &ceed); 1409 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1410 &basisq); 1411 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1412 &basisx); 1413 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1414 CEED_GAUSS_LOBATTO, &basisxc); 1415 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1416 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1417 CHKERRQ(ierr); 1418 1419 // CEED Restrictions 1420 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1421 qdatasizeVol, &restrictq, &restrictx, 1422 &restrictqdi); CHKERRQ(ierr); 1423 1424 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1425 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1426 1427 // Create the CEED vectors that will be needed in setup 1428 CeedInt NqptsVol; 1429 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1430 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1431 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1432 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1433 1434 // Create the Q-function that builds the quadrature data for the NS operator 1435 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1436 &qf_setupVol); 1437 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1438 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1439 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1440 1441 // Create the Q-function that sets the ICs of the operator 1442 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1443 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1444 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1445 1446 qf_rhsVol = NULL; 1447 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1448 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1449 problem->applyVol_rhs_loc, &qf_rhsVol); 1450 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1451 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1452 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1453 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1454 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1455 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1456 } 1457 1458 qf_ifunctionVol = NULL; 1459 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1460 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1461 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1462 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1463 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1464 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1465 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1466 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1467 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1468 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1469 } 1470 1471 // Create the operator that builds the quadrature data for the NS operator 1472 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1473 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1474 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1475 basisx, CEED_VECTOR_NONE); 1476 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1477 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1478 1479 // Create the operator that sets the ICs 1480 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1481 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1482 CeedOperatorSetField(op_ics, "q0", restrictq, 1483 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1484 1485 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1486 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1487 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1488 1489 if (qf_rhsVol) { // Create the RHS physics operator 1490 CeedOperator op; 1491 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1492 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1493 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1494 CeedOperatorSetField(op, "qdata", restrictqdi, 1495 CEED_BASIS_COLLOCATED, qdata); 1496 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1497 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1498 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1499 user->op_rhs_vol = op; 1500 } 1501 1502 if (qf_ifunctionVol) { // Create the IFunction operator 1503 CeedOperator op; 1504 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1505 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1506 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1507 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1508 CeedOperatorSetField(op, "qdata", restrictqdi, 1509 CEED_BASIS_COLLOCATED, qdata); 1510 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1511 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1512 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1513 user->op_ifunction_vol = op; 1514 } 1515 1516 // Set up CEED for the boundaries 1517 CeedInt height = 1; 1518 CeedInt dimSur = dim - height; 1519 CeedInt numP_Sur = degree + 1; 1520 CeedInt numQ_Sur = numP_Sur + qextraSur; 1521 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1522 CeedBasis basisxSur, basisxcSur, basisqSur; 1523 CeedInt NqptsSur; 1524 CeedQFunction qf_setupSur, qf_applySur; 1525 1526 // CEED bases for the boundaries 1527 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, 1528 CEED_GAUSS, 1529 &basisqSur); 1530 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1531 &basisxSur); 1532 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1533 CEED_GAUSS_LOBATTO, &basisxcSur); 1534 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1535 1536 // Create the Q-function that builds the quadrature data for the Surface operator 1537 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1538 &qf_setupSur); 1539 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1540 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1541 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1542 1543 // Creat Q-Function for Boundaries 1544 // -- Defined for Advection(2d) test cases 1545 qf_applySur = NULL; 1546 if (problem->applySur) { 1547 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1548 problem->applySur_loc, &qf_applySur); 1549 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1550 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1551 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1552 CeedQFunctionAddOutput(qf_applySur, "v", ncompq, CEED_EVAL_INTERP); 1553 } 1554 1555 // Create CEED Operator for the whole domain 1556 if (!implicit) 1557 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1558 user->op_rhs_vol, qf_applySur, 1559 qf_setupSur, height, numP_Sur, numQ_Sur, 1560 qdatasizeSur, NqptsSur, basisxSur, 1561 basisqSur, &user->op_rhs); 1562 CHKERRQ(ierr); 1563 if (implicit) 1564 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1565 user->op_ifunction_vol, qf_applySur, 1566 qf_setupSur, height, numP_Sur, numQ_Sur, 1567 qdatasizeSur, NqptsSur, basisxSur, 1568 basisqSur, &user->op_ifunction); 1569 CHKERRQ(ierr); 1570 // Set up contex for QFunctions 1571 CeedQFunctionContextCreate(ceed, &ctxSetup); 1572 CeedQFunctionContextSetData(ctxSetup, CEED_MEM_HOST, CEED_USE_POINTER, 1573 sizeof ctxSetupData, &ctxSetupData); 1574 if (qf_ics && problemChoice != NS_EULER_VORTEX) 1575 CeedQFunctionSetContext(qf_ics, ctxSetup); 1576 1577 CeedScalar ctxNSData[8] = {lambda, mu, k, cv, cp, g, Rd}; 1578 CeedQFunctionContextCreate(ceed, &ctxNS); 1579 CeedQFunctionContextSetData(ctxNS, CEED_MEM_HOST, CEED_USE_POINTER, 1580 sizeof ctxNSData, &ctxNSData); 1581 1582 struct Advection2dContext_ ctxAdvection2dData = { 1583 .CtauS = CtauS, 1584 .strong_form = strong_form, 1585 .stabilization = stab, 1586 }; 1587 CeedQFunctionContextCreate(ceed, &ctxAdvection2d); 1588 CeedQFunctionContextSetData(ctxAdvection2d, CEED_MEM_HOST, CEED_USE_POINTER, 1589 sizeof ctxAdvection2dData, &ctxAdvection2dData); 1590 1591 struct SurfaceContext_ ctxSurfaceData = { 1592 .E_wind = E_wind, 1593 .strong_form = strong_form, 1594 .implicit = implicit, 1595 }; 1596 CeedQFunctionContextCreate(ceed, &ctxSurface); 1597 CeedQFunctionContextSetData(ctxSurface, CEED_MEM_HOST, CEED_USE_POINTER, 1598 sizeof ctxSurfaceData, &ctxSurfaceData); 1599 1600 // Set up ctxEulerData structure 1601 ctxEulerData->time = 0.; 1602 ctxEulerData->currentTime = 0.; 1603 ctxEulerData->euler_test = euler_test; 1604 ctxEulerData->center[0] = center[0]; 1605 ctxEulerData->center[1] = center[1]; 1606 ctxEulerData->center[2] = center[2]; 1607 ctxEulerData->vortex_strength = vortex_strength; 1608 ctxEulerData->etv_mean_velocity[0] = etv_mean_velocity[0]; 1609 ctxEulerData->etv_mean_velocity[1] = etv_mean_velocity[1]; 1610 ctxEulerData->etv_mean_velocity[2] = etv_mean_velocity[2]; 1611 ctxEulerData->stabilization = stab; 1612 ctxEulerData->implicit = implicit; 1613 user->ctxEulerData = ctxEulerData; 1614 1615 CeedQFunctionContextCreate(ceed, &ctxEuler); 1616 CeedQFunctionContextSetData(ctxEuler, CEED_MEM_HOST, CEED_USE_POINTER, 1617 sizeof *ctxEulerData, ctxEulerData); 1618 1619 switch (problemChoice) { 1620 case NS_DENSITY_CURRENT: 1621 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxNS); 1622 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxNS); 1623 break; 1624 case NS_ADVECTION: 1625 case NS_ADVECTION2D: 1626 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxAdvection2d); 1627 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxAdvection2d); 1628 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, ctxSurface); 1629 break; 1630 case NS_EULER_VORTEX: 1631 if (qf_ics) CeedQFunctionSetContext(qf_ics, ctxEuler); 1632 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, ctxEuler); 1633 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, ctxEuler); 1634 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, ctxEuler); 1635 } 1636 1637 // Set up PETSc context 1638 // Set up units structure 1639 units->meter = meter; 1640 units->kilogram = kilogram; 1641 units->second = second; 1642 units->Kelvin = Kelvin; 1643 units->Pascal = Pascal; 1644 units->JperkgK = JperkgK; 1645 units->mpersquareds = mpersquareds; 1646 units->WpermK = WpermK; 1647 units->kgpercubicm = kgpercubicm; 1648 units->kgpersquaredms = kgpersquaredms; 1649 units->Joulepercubicm = Joulepercubicm; 1650 units->Joule = Joule; 1651 1652 // Set up user structure 1653 user->comm = comm; 1654 user->outputfreq = outputfreq; 1655 user->contsteps = contsteps; 1656 user->units = units; 1657 user->dm = dm; 1658 user->dmviz = dmviz; 1659 user->interpviz = interpviz; 1660 user->ceed = ceed; 1661 1662 // Calculate qdata and ICs 1663 // Set up state global and local vectors 1664 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1665 1666 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1667 1668 // Apply Setup Ceed Operators 1669 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1670 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1671 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1672 user->M); CHKERRQ(ierr); 1673 1674 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1675 ctxSetup, 0.0); CHKERRQ(ierr); 1676 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1677 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1678 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1679 // slower execution. 1680 Vec Qbc; 1681 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1682 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1683 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1684 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1685 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1686 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1687 ierr = PetscObjectComposeFunction((PetscObject)dm, 1688 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1689 CHKERRQ(ierr); 1690 } 1691 1692 MPI_Comm_rank(comm, &rank); 1693 if (!rank) {ierr = PetscMkdir(user->outputdir); CHKERRQ(ierr);} 1694 // Gather initial Q values 1695 // In case of continuation of simulation, set up initial values from binary file 1696 if (contsteps) { // continue from existent solution 1697 PetscViewer viewer; 1698 char filepath[PETSC_MAX_PATH_LEN]; 1699 // Read input 1700 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1701 user->outputdir); 1702 CHKERRQ(ierr); 1703 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1704 CHKERRQ(ierr); 1705 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1706 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1707 } 1708 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1709 1710 // Create and setup TS 1711 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1712 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1713 if (implicit) { 1714 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1715 if (user->op_ifunction) { 1716 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1717 } else { // Implicit integrators can fall back to using an RHSFunction 1718 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1719 } 1720 } else { 1721 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1722 "Problem does not provide RHSFunction"); 1723 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1724 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1725 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1726 } 1727 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1728 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1729 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1730 if (test) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1731 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1732 ierr = TSAdaptSetStepLimits(adapt, 1733 1.e-12 * units->second, 1734 1.e2 * units->second); CHKERRQ(ierr); 1735 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1736 if (!contsteps) { // print initial condition 1737 if (!test) { 1738 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1739 } 1740 } else { // continue from time of last output 1741 PetscReal time; 1742 PetscInt count; 1743 PetscViewer viewer; 1744 char filepath[PETSC_MAX_PATH_LEN]; 1745 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1746 user->outputdir); CHKERRQ(ierr); 1747 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1748 CHKERRQ(ierr); 1749 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1750 CHKERRQ(ierr); 1751 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1752 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1753 } 1754 if (!test) { 1755 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1756 } 1757 1758 // Solve 1759 start = MPI_Wtime(); 1760 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1761 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1762 cpu_time_used = MPI_Wtime() - start; 1763 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1764 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1765 comm); CHKERRQ(ierr); 1766 if (!test) { 1767 ierr = PetscPrintf(PETSC_COMM_WORLD, 1768 "Time taken for solution (sec): %g\n", 1769 (double)cpu_time_used); CHKERRQ(ierr); 1770 } 1771 1772 // Get error 1773 if (problem->non_zero_time && !test) { 1774 Vec Qexact, Qexactloc; 1775 PetscReal rel_error, norm_error, norm_exact; 1776 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1777 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1778 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1779 1780 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1781 restrictq, ctxSetup, ftime); CHKERRQ(ierr); 1782 ierr = VecNorm(Qexact, NORM_1, &norm_exact); CHKERRQ(ierr); 1783 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1784 ierr = VecNorm(Q, NORM_1, &norm_error); CHKERRQ(ierr); 1785 rel_error = norm_error / norm_exact; 1786 CeedVectorDestroy(&q0ceed); 1787 ierr = PetscPrintf(PETSC_COMM_WORLD, 1788 "Relative Error: %g\n", 1789 (double)rel_error); CHKERRQ(ierr); 1790 // Clean up vectors 1791 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1792 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1793 } 1794 1795 // Output Statistics 1796 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1797 if (!test) { 1798 ierr = PetscPrintf(PETSC_COMM_WORLD, 1799 "Time integrator took %D time steps to reach final time %g\n", 1800 steps, (double)ftime); CHKERRQ(ierr); 1801 } 1802 // Output numerical values from command line 1803 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1804 1805 // Compare reference solution values with current test run for CI 1806 if (test) { 1807 PetscViewer viewer; 1808 // Read reference file 1809 Vec Qref; 1810 PetscReal error, Qrefnorm; 1811 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1812 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1813 CHKERRQ(ierr); 1814 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1815 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1816 1817 // Compute error with respect to reference solution 1818 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1819 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1820 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1821 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1822 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1823 // Check error 1824 if (error > testtol) { 1825 ierr = PetscPrintf(PETSC_COMM_WORLD, 1826 "Test failed with error norm %g\n", 1827 (double)error); CHKERRQ(ierr); 1828 } 1829 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1830 } 1831 1832 // Clean up libCEED 1833 CeedVectorDestroy(&qdata); 1834 CeedVectorDestroy(&user->qceed); 1835 CeedVectorDestroy(&user->qdotceed); 1836 CeedVectorDestroy(&user->gceed); 1837 CeedVectorDestroy(&xcorners); 1838 CeedBasisDestroy(&basisq); 1839 CeedBasisDestroy(&basisx); 1840 CeedBasisDestroy(&basisxc); 1841 CeedElemRestrictionDestroy(&restrictq); 1842 CeedElemRestrictionDestroy(&restrictx); 1843 CeedElemRestrictionDestroy(&restrictqdi); 1844 CeedQFunctionDestroy(&qf_setupVol); 1845 CeedQFunctionDestroy(&qf_ics); 1846 CeedQFunctionDestroy(&qf_rhsVol); 1847 CeedQFunctionDestroy(&qf_ifunctionVol); 1848 CeedQFunctionContextDestroy(&ctxSetup); 1849 CeedQFunctionContextDestroy(&ctxNS); 1850 CeedQFunctionContextDestroy(&ctxAdvection2d); 1851 CeedQFunctionContextDestroy(&ctxSurface); 1852 CeedQFunctionContextDestroy(&ctxEuler); 1853 CeedOperatorDestroy(&op_setupVol); 1854 CeedOperatorDestroy(&op_ics); 1855 CeedOperatorDestroy(&user->op_rhs_vol); 1856 CeedOperatorDestroy(&user->op_ifunction_vol); 1857 CeedDestroy(&ceed); 1858 CeedBasisDestroy(&basisqSur); 1859 CeedBasisDestroy(&basisxSur); 1860 CeedBasisDestroy(&basisxcSur); 1861 CeedQFunctionDestroy(&qf_setupSur); 1862 CeedQFunctionDestroy(&qf_applySur); 1863 CeedOperatorDestroy(&user->op_rhs); 1864 CeedOperatorDestroy(&user->op_ifunction); 1865 1866 // Clean up PETSc 1867 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1868 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1869 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1870 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1871 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1872 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1873 ierr = PetscFree(units); CHKERRQ(ierr); 1874 ierr = PetscFree(user); CHKERRQ(ierr); 1875 ierr = PetscFree(ctxEulerData); CHKERRQ(ierr); 1876 return PetscFinalize(); 1877 } 1878