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