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