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