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