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