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 if (problem->bc == Exact_Euler) 855 bc->nwall = bc->nslip[0] = bc->nslip[1] = bc->nslip[2] = 0; 856 { 857 PetscInt comps[1] = {1}; 858 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipx", "Face Sets", 0, 859 1, comps, (void(*)(void))NULL, NULL, bc->nslip[0], 860 bc->slips[0], ctxSetupData); CHKERRQ(ierr); 861 comps[0] = 2; 862 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipy", "Face Sets", 0, 863 1, comps, (void(*)(void))NULL, NULL, bc->nslip[1], 864 bc->slips[1], ctxSetupData); CHKERRQ(ierr); 865 comps[0] = 3; 866 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", "Face Sets", 0, 867 1, comps, (void(*)(void))NULL, NULL, bc->nslip[2], 868 bc->slips[2], ctxSetupData); CHKERRQ(ierr); 869 } 870 if (bc->userbc == PETSC_TRUE) { 871 for (PetscInt c = 0; c < 3; c++) { 872 for (PetscInt s = 0; s < bc->nslip[c]; s++) { 873 for (PetscInt w = 0; w < bc->nwall; w++) { 874 if (bc->slips[c][s] == bc->walls[w]) 875 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, 876 "Boundary condition already set on face %D!\n", 877 bc->walls[w]); 878 } 879 } 880 } 881 } 882 // Wall boundary conditions are zero energy density and zero flux for 883 // velocity in advection/advection2d, and zero velocity and zero flux 884 // for mass density and energy density in density_current 885 { 886 if (problem->bc == Exact_Advection || problem->bc == Exact_Advection2d) { 887 PetscInt comps[1] = {4}; 888 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 889 1, comps, (void(*)(void))problem->bc, NULL, 890 bc->nwall, bc->walls, ctxSetup); CHKERRQ(ierr); 891 } else if (problem->bc == Exact_DC) { 892 PetscInt comps[3] = {1, 2, 3}; 893 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "Face Sets", 0, 894 3, comps, (void(*)(void))problem->bc, NULL, 895 bc->nwall, bc->walls, ctxSetupData); CHKERRQ(ierr); 896 } else if (problem->bc == Exact_Euler) { 897 PetscInt bcMMS[4] = {1, 2, 3, 4}; 898 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "mms", "Face Sets", 0, 899 0, NULL, (void(*)(void))problem->bc, 900 4, bcMMS, ctxSetup); CHKERRQ(ierr); 901 } else 902 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, 903 "Undefined boundary conditions for this problem"); 904 } 905 ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); 906 CHKERRQ(ierr); 907 ierr = PetscFEDestroy(&fe); CHKERRQ(ierr); 908 } 909 { 910 // Empty name for conserved field (because there is only one field) 911 PetscSection section; 912 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 913 ierr = PetscSectionSetFieldName(section, 0, ""); CHKERRQ(ierr); 914 ierr = PetscSectionSetComponentName(section, 0, 0, "Density"); 915 CHKERRQ(ierr); 916 ierr = PetscSectionSetComponentName(section, 0, 1, "MomentumX"); 917 CHKERRQ(ierr); 918 ierr = PetscSectionSetComponentName(section, 0, 2, "MomentumY"); 919 CHKERRQ(ierr); 920 ierr = PetscSectionSetComponentName(section, 0, 3, "MomentumZ"); 921 CHKERRQ(ierr); 922 ierr = PetscSectionSetComponentName(section, 0, 4, "EnergyDensity"); 923 CHKERRQ(ierr); 924 } 925 PetscFunctionReturn(0); 926 } 927 928 int main(int argc, char **argv) { 929 PetscInt ierr; 930 MPI_Comm comm; 931 DM dm, dmcoord, dmviz; 932 Mat interpviz; 933 TS ts; 934 TSAdapt adapt; 935 User user; 936 Units units; 937 char ceedresource[4096] = "/cpu/self"; 938 PetscInt localNelemVol, lnodes, gnodes, steps; 939 const PetscInt ncompq = 5; 940 PetscMPIInt rank; 941 PetscScalar ftime; 942 Vec Q, Qloc, Xloc; 943 Ceed ceed; 944 CeedInt numP, numQ; 945 CeedVector xcorners, qdata, q0ceed; 946 CeedBasis basisx, basisxc, basisq; 947 CeedElemRestriction restrictx, restrictq, restrictqdi; 948 CeedQFunction qf_setupVol, qf_ics, qf_rhsVol, qf_ifunctionVol; 949 CeedQFunctionContext ctxSetup, ctxNS, ctxAdvection2d, ctxSurface; 950 CeedOperator op_setupVol, op_ics; 951 CeedScalar Rd; 952 CeedMemType memtyperequested; 953 PetscScalar WpermK, Pascal, JperkgK, mpersquareds, kgpercubicm, 954 kgpersquaredms, Joulepercubicm, Joule; 955 problemType problemChoice; 956 problemData *problem = NULL; 957 WindType wind_type; 958 StabilizationType stab; 959 PetscBool implicit; 960 PetscInt viz_refine = 0; 961 struct SimpleBC_ bc = { 962 .nslip = {2, 2, 2}, 963 .slips = {{5, 6}, {3, 4}, {1, 2}} 964 }; 965 double start, cpu_time_used; 966 // Test variables 967 PetscBool test; 968 PetscScalar testtol = 0.; 969 char filepath[PETSC_MAX_PATH_LEN]; 970 // Check PETSc CUDA support 971 PetscBool petschavecuda, setmemtyperequest = PETSC_FALSE; 972 // *INDENT-OFF* 973 #ifdef PETSC_HAVE_CUDA 974 petschavecuda = PETSC_TRUE; 975 #else 976 petschavecuda = PETSC_FALSE; 977 #endif 978 // *INDENT-ON* 979 980 // Create the libCEED contexts 981 PetscScalar meter = 1e-2; // 1 meter in scaled length units 982 PetscScalar second = 1e-2; // 1 second in scaled time units 983 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 984 PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units 985 CeedScalar theta0 = 300.; // K 986 CeedScalar thetaC = -15.; // K 987 CeedScalar P0 = 1.e5; // Pa 988 CeedScalar E_wind = 1.e6; // J 989 CeedScalar N = 0.01; // 1/s 990 CeedScalar cv = 717.; // J/(kg K) 991 CeedScalar cp = 1004.; // J/(kg K) 992 CeedScalar vortex_strength = 5.; // - 993 CeedScalar rho_enter = 1.2; // Kg/m^3 994 CeedScalar g = 9.81; // m/s^2 995 CeedScalar lambda = -2./3.; // - 996 CeedScalar mu = 75.; // Pa s, dynamic viscosity 997 // mu = 75 is not physical for air, but is good for numerical stability 998 CeedScalar k = 0.02638; // W/(m K) 999 CeedScalar CtauS = 0.; // dimensionless 1000 CeedScalar strong_form = 0.; // [0,1] 1001 PetscScalar lx = 8000.; // m 1002 PetscScalar ly = 8000.; // m 1003 PetscScalar lz = 4000.; // m 1004 CeedScalar rc = 1000.; // m (Radius of bubble) 1005 PetscScalar resx = 1000.; // m (resolution in x) 1006 PetscScalar resy = 1000.; // m (resolution in y) 1007 PetscScalar resz = 1000.; // m (resolution in z) 1008 PetscInt outputfreq = 10; // - 1009 PetscInt contsteps = 0; // - 1010 PetscInt degree = 1; // - 1011 PetscInt qextra = 2; // - 1012 PetscInt qextraSur = 2; // - 1013 PetscReal center[3], dc_axis[3] = {0, 0, 0}, wind[3] = {1., 0, 0}, 1014 u_enter[3] = {-1.2, 0, 0}; 1015 1016 ierr = PetscInitialize(&argc, &argv, NULL, help); 1017 if (ierr) return ierr; 1018 1019 // Allocate PETSc context 1020 ierr = PetscCalloc1(1, &user); CHKERRQ(ierr); 1021 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 1022 1023 // Parse command line options 1024 comm = PETSC_COMM_WORLD; 1025 ierr = PetscOptionsBegin(comm, NULL, "Navier-Stokes in PETSc with libCEED", 1026 NULL); CHKERRQ(ierr); 1027 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 1028 NULL, ceedresource, ceedresource, 1029 sizeof(ceedresource), NULL); CHKERRQ(ierr); 1030 ierr = PetscOptionsBool("-test", "Run in test mode", 1031 NULL, test=PETSC_FALSE, &test, NULL); CHKERRQ(ierr); 1032 ierr = PetscOptionsScalar("-compare_final_state_atol", 1033 "Test absolute tolerance", 1034 NULL, testtol, &testtol, NULL); CHKERRQ(ierr); 1035 ierr = PetscOptionsString("-compare_final_state_filename", "Test filename", 1036 NULL, filepath, filepath, 1037 sizeof(filepath), NULL); CHKERRQ(ierr); 1038 problemChoice = NS_DENSITY_CURRENT; 1039 ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, 1040 problemTypes, (PetscEnum)problemChoice, 1041 (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); 1042 problem = &problemOptions[problemChoice]; 1043 ierr = PetscOptionsEnum("-problem_advection_wind", "Wind type in Advection", 1044 NULL, WindTypes, 1045 (PetscEnum)(wind_type = ADVECTION_WIND_ROTATION), 1046 (PetscEnum *)&wind_type, NULL); CHKERRQ(ierr); 1047 PetscInt n = problem->dim; 1048 PetscBool userWind; 1049 ierr = PetscOptionsRealArray("-problem_advection_wind_translation", 1050 "Constant wind vector", 1051 NULL, wind, &n, &userWind); CHKERRQ(ierr); 1052 if (wind_type == ADVECTION_WIND_ROTATION && userWind) { 1053 ierr = PetscPrintf(comm, 1054 "Warning! Use -problem_advection_wind_translation only with -problem_advection_wind translation\n"); 1055 CHKERRQ(ierr); 1056 } 1057 if (wind_type == ADVECTION_WIND_TRANSLATION 1058 && (problemChoice == NS_DENSITY_CURRENT || 1059 problemChoice == NS_EULER_VORTEX)) { 1060 SETERRQ(comm, PETSC_ERR_ARG_INCOMP, 1061 "-problem_advection_wind translation is not defined for -problem density_current or -problem euler_vortex"); 1062 } 1063 ierr = PetscOptionsRealArray("-problem_euler_vortex_velocity", 1064 "Incoming velocity vector", 1065 NULL, u_enter, &n, NULL); CHKERRQ(ierr); 1066 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 1067 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 1068 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 1069 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 1070 NULL, implicit=PETSC_FALSE, &implicit, NULL); 1071 CHKERRQ(ierr); 1072 if (!implicit && stab != STAB_NONE) { 1073 ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); 1074 CHKERRQ(ierr); 1075 } 1076 { 1077 PetscInt len; 1078 PetscBool flg; 1079 ierr = PetscOptionsIntArray("-bc_wall", 1080 "Use wall boundary conditions on this list of faces", 1081 NULL, bc.walls, 1082 (len = sizeof(bc.walls) / sizeof(bc.walls[0]), 1083 &len), &flg); CHKERRQ(ierr); 1084 if (flg) { 1085 bc.nwall = len; 1086 // Using a no-slip wall disables automatic slip walls (they must be set explicitly) 1087 bc.nslip[0] = bc.nslip[1] = bc.nslip[2] = 0; 1088 } 1089 for (PetscInt j=0; j<3; j++) { 1090 const char *flags[3] = {"-bc_slip_x", "-bc_slip_y", "-bc_slip_z"}; 1091 ierr = PetscOptionsIntArray(flags[j], 1092 "Use slip boundary conditions on this list of faces", 1093 NULL, bc.slips[j], 1094 (len = sizeof(bc.slips[j]) / sizeof(bc.slips[j][0]), 1095 &len), &flg); 1096 CHKERRQ(ierr); 1097 if (flg) { 1098 bc.nslip[j] = len; 1099 bc.userbc = PETSC_TRUE; 1100 } 1101 } 1102 } 1103 ierr = PetscOptionsInt("-viz_refine", 1104 "Regular refinement levels for visualization", 1105 NULL, viz_refine, &viz_refine, NULL); 1106 CHKERRQ(ierr); 1107 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 1108 NULL, meter, &meter, NULL); CHKERRQ(ierr); 1109 meter = fabs(meter); 1110 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 1111 NULL, second, &second, NULL); CHKERRQ(ierr); 1112 second = fabs(second); 1113 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 1114 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 1115 kilogram = fabs(kilogram); 1116 ierr = PetscOptionsScalar("-units_Kelvin", 1117 "1 Kelvin in scaled temperature units", 1118 NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr); 1119 Kelvin = fabs(Kelvin); 1120 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", 1121 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 1122 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 1123 NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr); 1124 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", 1125 NULL, P0, &P0, NULL); CHKERRQ(ierr); 1126 ierr = PetscOptionsScalar("-E_wind", "Total energy of inflow wind", 1127 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 1128 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", 1129 NULL, N, &N, NULL); CHKERRQ(ierr); 1130 ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume", 1131 NULL, cv, &cv, NULL); CHKERRQ(ierr); 1132 ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure", 1133 NULL, cp, &cp, NULL); CHKERRQ(ierr); 1134 PetscBool userVortex; 1135 ierr = PetscOptionsScalar("-vortex_strength", "Strength of Vortex", 1136 NULL, vortex_strength, &vortex_strength, &userVortex); 1137 CHKERRQ(ierr); 1138 if (problemChoice != NS_EULER_VORTEX && userVortex) { 1139 ierr = PetscPrintf(comm, 1140 "Warning! Use -vortex_strength only with -problem euler_vortex\n"); 1141 CHKERRQ(ierr); 1142 } 1143 ierr = PetscOptionsScalar("-problem_euler_vortex_rho", "Incoming density", 1144 NULL, rho_enter, &rho_enter, NULL); 1145 CHKERRQ(ierr); 1146 ierr = PetscOptionsScalar("-g", "Gravitational acceleration", 1147 NULL, g, &g, NULL); CHKERRQ(ierr); 1148 ierr = PetscOptionsScalar("-lambda", 1149 "Stokes hypothesis second viscosity coefficient", 1150 NULL, lambda, &lambda, NULL); CHKERRQ(ierr); 1151 ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient", 1152 NULL, mu, &mu, NULL); CHKERRQ(ierr); 1153 ierr = PetscOptionsScalar("-k", "Thermal conductivity", 1154 NULL, k, &k, NULL); CHKERRQ(ierr); 1155 ierr = PetscOptionsScalar("-CtauS", 1156 "Scale coefficient for tau (nondimensional)", 1157 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 1158 if (stab == STAB_NONE && CtauS != 0) { 1159 ierr = PetscPrintf(comm, 1160 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 1161 CHKERRQ(ierr); 1162 } 1163 ierr = PetscOptionsScalar("-strong_form", 1164 "Strong (1) or weak/integrated by parts (0) advection residual", 1165 NULL, strong_form, &strong_form, NULL); 1166 CHKERRQ(ierr); 1167 if (problemChoice == NS_DENSITY_CURRENT && (CtauS != 0 || strong_form != 0)) { 1168 ierr = PetscPrintf(comm, 1169 "Warning! Problem density_current does not support -CtauS or -strong_form\n"); 1170 CHKERRQ(ierr); 1171 } 1172 ierr = PetscOptionsScalar("-lx", "Length scale in x direction", 1173 NULL, lx, &lx, NULL); CHKERRQ(ierr); 1174 ierr = PetscOptionsScalar("-ly", "Length scale in y direction", 1175 NULL, ly, &ly, NULL); CHKERRQ(ierr); 1176 ierr = PetscOptionsScalar("-lz", "Length scale in z direction", 1177 NULL, lz, &lz, NULL); CHKERRQ(ierr); 1178 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 1179 NULL, rc, &rc, NULL); CHKERRQ(ierr); 1180 ierr = PetscOptionsScalar("-resx","Target resolution in x", 1181 NULL, resx, &resx, NULL); CHKERRQ(ierr); 1182 ierr = PetscOptionsScalar("-resy","Target resolution in y", 1183 NULL, resy, &resy, NULL); CHKERRQ(ierr); 1184 ierr = PetscOptionsScalar("-resz","Target resolution in z", 1185 NULL, resz, &resz, NULL); CHKERRQ(ierr); 1186 n = problem->dim; 1187 center[0] = 0.5 * lx; 1188 center[1] = 0.5 * ly; 1189 center[2] = 0.5 * lz; 1190 ierr = PetscOptionsRealArray("-center", "Location of bubble center", 1191 NULL, center, &n, NULL); CHKERRQ(ierr); 1192 n = problem->dim; 1193 ierr = PetscOptionsRealArray("-dc_axis", 1194 "Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric", 1195 NULL, dc_axis, &n, NULL); CHKERRQ(ierr); 1196 { 1197 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + 1198 PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 1199 if (norm > 0) { 1200 for (int i=0; i<3; i++) dc_axis[i] /= norm; 1201 } 1202 } 1203 ierr = PetscOptionsInt("-output_freq", 1204 "Frequency of output, in number of steps", 1205 NULL, outputfreq, &outputfreq, NULL); CHKERRQ(ierr); 1206 ierr = PetscOptionsInt("-continue", "Continue from previous solution", 1207 NULL, contsteps, &contsteps, NULL); CHKERRQ(ierr); 1208 ierr = PetscOptionsInt("-degree", "Polynomial degree of finite elements", 1209 NULL, degree, °ree, NULL); CHKERRQ(ierr); 1210 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 1211 NULL, qextra, &qextra, NULL); CHKERRQ(ierr); 1212 PetscBool userQextraSur; 1213 ierr = PetscOptionsInt("-qextra_boundary", 1214 "Number of extra quadrature points on in/outflow faces", 1215 NULL, qextraSur, &qextraSur, &userQextraSur); 1216 CHKERRQ(ierr); 1217 if ((wind_type == ADVECTION_WIND_ROTATION 1218 || problemChoice == NS_DENSITY_CURRENT) && userQextraSur) { 1219 ierr = PetscPrintf(comm, 1220 "Warning! Use -qextra_boundary only with -problem_advection_wind translation\n"); 1221 CHKERRQ(ierr); 1222 } 1223 ierr = PetscStrncpy(user->outputdir, ".", 2); CHKERRQ(ierr); 1224 ierr = PetscOptionsString("-output_dir", "Output directory", 1225 NULL, user->outputdir, user->outputdir, 1226 sizeof(user->outputdir), NULL); CHKERRQ(ierr); 1227 memtyperequested = petschavecuda ? CEED_MEM_DEVICE : CEED_MEM_HOST; 1228 ierr = PetscOptionsEnum("-memtype", 1229 "CEED MemType requested", NULL, 1230 memTypes, (PetscEnum)memtyperequested, 1231 (PetscEnum *)&memtyperequested, &setmemtyperequest); 1232 CHKERRQ(ierr); 1233 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 1234 1235 // Define derived units 1236 Pascal = kilogram / (meter * PetscSqr(second)); 1237 JperkgK = PetscSqr(meter) / (PetscSqr(second) * Kelvin); 1238 mpersquareds = meter / PetscSqr(second); 1239 WpermK = kilogram * meter / (pow(second,3) * Kelvin); 1240 kgpercubicm = kilogram / pow(meter,3); 1241 kgpersquaredms = kilogram / (PetscSqr(meter) * second); 1242 Joulepercubicm = kilogram / (meter * PetscSqr(second)); 1243 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 1244 1245 // Scale variables to desired units 1246 theta0 *= Kelvin; 1247 thetaC *= Kelvin; 1248 P0 *= Pascal; 1249 E_wind *= Joule; 1250 rho_enter *= kgpercubicm; 1251 N *= (1./second); 1252 cv *= JperkgK; 1253 cp *= JperkgK; 1254 Rd = cp - cv; 1255 g *= mpersquareds; 1256 mu *= Pascal * second; 1257 k *= WpermK; 1258 lx = fabs(lx) * meter; 1259 ly = fabs(ly) * meter; 1260 lz = fabs(lz) * meter; 1261 rc = fabs(rc) * meter; 1262 resx = fabs(resx) * meter; 1263 resy = fabs(resy) * meter; 1264 resz = fabs(resz) * meter; 1265 for (int i=0; i<3; i++) center[i] *= meter; 1266 1267 const CeedInt dim = problem->dim, ncompx = problem->dim, 1268 qdatasizeVol = problem->qdatasizeVol; 1269 // Set up the libCEED context 1270 struct SetupContext_ ctxSetupData = { 1271 .theta0 = theta0, 1272 .thetaC = thetaC, 1273 .P0 = P0, 1274 .N = N, 1275 .cv = cv, 1276 .cp = cp, 1277 .Rd = Rd, 1278 .g = g, 1279 .rc = rc, 1280 .lx = lx, 1281 .ly = ly, 1282 .lz = lz, 1283 .center[0] = center[0], 1284 .center[1] = center[1], 1285 .center[2] = center[2], 1286 .dc_axis[0] = dc_axis[0], 1287 .dc_axis[1] = dc_axis[1], 1288 .dc_axis[2] = dc_axis[2], 1289 .wind[0] = wind[0], 1290 .wind[1] = wind[1], 1291 .wind[2] = wind[2], 1292 .time = 0, 1293 .vortex_strength = vortex_strength, 1294 .wind_type = wind_type, 1295 }; 1296 1297 // Create the mesh 1298 { 1299 const PetscReal scale[3] = {lx, ly, lz}; 1300 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, NULL, NULL, scale, 1301 NULL, PETSC_TRUE, &dm); 1302 CHKERRQ(ierr); 1303 } 1304 1305 // Distribute the mesh over processes 1306 { 1307 DM dmDist = NULL; 1308 PetscPartitioner part; 1309 1310 ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr); 1311 ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr); 1312 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist); CHKERRQ(ierr); 1313 if (dmDist) { 1314 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1315 dm = dmDist; 1316 } 1317 } 1318 ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); 1319 1320 // Setup DM 1321 ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); 1322 ierr = DMSetFromOptions(dm); CHKERRQ(ierr); 1323 ierr = SetUpDM(dm, problem, degree, &bc, &ctxSetupData); 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 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 CeedQFunctionSetContext(qf_ics, ctxSetup); 1599 1600 CeedScalar ctxNSData[8] = {lambda, mu, k, cv, cp, g, Rd}; 1601 CeedQFunctionContextCreate(ceed, &ctxNS); 1602 CeedQFunctionContextSetData(ctxNS, CEED_MEM_HOST, CEED_USE_POINTER, 1603 sizeof ctxNSData, &ctxNSData); 1604 1605 struct Advection2dContext_ ctxAdvection2dData = { 1606 .CtauS = CtauS, 1607 .strong_form = strong_form, 1608 .stabilization = stab, 1609 }; 1610 CeedQFunctionContextCreate(ceed, &ctxAdvection2d); 1611 CeedQFunctionContextSetData(ctxAdvection2d, CEED_MEM_HOST, CEED_USE_POINTER, 1612 sizeof ctxAdvection2dData, &ctxAdvection2dData); 1613 1614 struct SurfaceContext_ ctxSurfaceData = { 1615 .E_wind = E_wind, 1616 .strong_form = strong_form, 1617 .implicit = implicit, 1618 }; 1619 struct EulerContext_ ctxEuler = { 1620 .rho_enter = rho_enter, 1621 .u_enter[0] = u_enter[0], 1622 .u_enter[1] = u_enter[1], 1623 .u_enter[2] = u_enter[2], 1624 .implicit = implicit, 1625 }; 1626 CeedQFunctionContextCreate(ceed, &ctxSurface); 1627 CeedQFunctionContextSetData(ctxSurface, CEED_MEM_HOST, CEED_USE_POINTER, 1628 sizeof ctxSurfaceData, &ctxSurfaceData); 1629 1630 switch (problemChoice) { 1631 case NS_DENSITY_CURRENT: 1632 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1633 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1634 sizeof ctxNS); 1635 break; 1636 case NS_ADVECTION: 1637 case NS_ADVECTION2D: 1638 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1639 sizeof ctxAdvection2d); 1640 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1641 sizeof ctxAdvection2d); 1642 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &ctxSurface, 1643 sizeof ctxSurface); 1644 case NS_EULER_VORTEX: 1645 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &ctxSetup, 1646 sizeof ctxSetup); 1647 } 1648 1649 // Set up PETSc context 1650 // Set up units structure 1651 units->meter = meter; 1652 units->kilogram = kilogram; 1653 units->second = second; 1654 units->Kelvin = Kelvin; 1655 units->Pascal = Pascal; 1656 units->JperkgK = JperkgK; 1657 units->mpersquareds = mpersquareds; 1658 units->WpermK = WpermK; 1659 units->kgpercubicm = kgpercubicm; 1660 units->kgpersquaredms = kgpersquaredms; 1661 units->Joulepercubicm = Joulepercubicm; 1662 units->Joule = Joule; 1663 1664 // Set up user structure 1665 user->comm = comm; 1666 user->outputfreq = outputfreq; 1667 user->contsteps = contsteps; 1668 user->units = units; 1669 user->dm = dm; 1670 user->dmviz = dmviz; 1671 user->interpviz = interpviz; 1672 user->ceed = ceed; 1673 1674 // Calculate qdata and ICs 1675 // Set up state global and local vectors 1676 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1677 1678 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1679 1680 // Apply Setup Ceed Operators 1681 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1682 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1683 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1684 user->M); CHKERRQ(ierr); 1685 1686 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1687 ctxSetup, 0.0); CHKERRQ(ierr); 1688 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1689 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1690 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1691 // slower execution. 1692 Vec Qbc; 1693 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1694 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1695 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1696 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1697 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1698 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1699 ierr = PetscObjectComposeFunction((PetscObject)dm, 1700 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1701 CHKERRQ(ierr); 1702 } 1703 1704 MPI_Comm_rank(comm, &rank); 1705 if (!rank) {ierr = PetscMkdir(user->outputdir); CHKERRQ(ierr);} 1706 // Gather initial Q values 1707 // In case of continuation of simulation, set up initial values from binary file 1708 if (contsteps) { // continue from existent solution 1709 PetscViewer viewer; 1710 char filepath[PETSC_MAX_PATH_LEN]; 1711 // Read input 1712 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1713 user->outputdir); 1714 CHKERRQ(ierr); 1715 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1716 CHKERRQ(ierr); 1717 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1718 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1719 } 1720 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1721 1722 // Create and setup TS 1723 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1724 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1725 if (implicit) { 1726 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1727 if (user->op_ifunction) { 1728 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1729 } else { // Implicit integrators can fall back to using an RHSFunction 1730 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1731 } 1732 } else { 1733 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1734 "Problem does not provide RHSFunction"); 1735 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1736 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1737 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1738 } 1739 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1740 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1741 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1742 if (test) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1743 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1744 ierr = TSAdaptSetStepLimits(adapt, 1745 1.e-12 * units->second, 1746 1.e2 * units->second); CHKERRQ(ierr); 1747 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1748 if (!contsteps) { // print initial condition 1749 if (!test) { 1750 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1751 } 1752 } else { // continue from time of last output 1753 PetscReal time; 1754 PetscInt count; 1755 PetscViewer viewer; 1756 char filepath[PETSC_MAX_PATH_LEN]; 1757 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1758 user->outputdir); CHKERRQ(ierr); 1759 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1760 CHKERRQ(ierr); 1761 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1762 CHKERRQ(ierr); 1763 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1764 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1765 } 1766 if (!test) { 1767 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1768 } 1769 // Get the current time 1770 PetscReal timeNow; 1771 ierr = TSGetTime(ts,&timeNow);CHKERRQ(ierr); 1772 ctxSetup.time = timeNow; 1773 1774 // Solve 1775 start = MPI_Wtime(); 1776 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1777 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1778 cpu_time_used = MPI_Wtime() - start; 1779 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1780 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1781 comm); CHKERRQ(ierr); 1782 if (!test) { 1783 ierr = PetscPrintf(PETSC_COMM_WORLD, 1784 "Time taken for solution (sec): %g\n", 1785 (double)cpu_time_used); CHKERRQ(ierr); 1786 } 1787 1788 // Get error 1789 if (problem->non_zero_time && !test) { 1790 Vec Qexact, Qexactloc; 1791 PetscReal norm; 1792 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1793 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1794 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1795 1796 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1797 restrictq, ctxSetup, ftime); CHKERRQ(ierr); 1798 1799 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1800 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1801 CeedVectorDestroy(&q0ceed); 1802 ierr = PetscPrintf(PETSC_COMM_WORLD, 1803 "Max Error: %g\n", 1804 (double)norm); CHKERRQ(ierr); 1805 // Clean up vectors 1806 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1807 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1808 } 1809 1810 // Output Statistics 1811 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1812 if (!test) { 1813 ierr = PetscPrintf(PETSC_COMM_WORLD, 1814 "Time integrator took %D time steps to reach final time %g\n", 1815 steps, (double)ftime); CHKERRQ(ierr); 1816 } 1817 // Output numerical values from command line 1818 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1819 1820 // Compare reference solution values with current test run for CI 1821 if (test) { 1822 PetscViewer viewer; 1823 // Read reference file 1824 Vec Qref; 1825 PetscReal error, Qrefnorm; 1826 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1827 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1828 CHKERRQ(ierr); 1829 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1830 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1831 1832 // Compute error with respect to reference solution 1833 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1834 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1835 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1836 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1837 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1838 // Check error 1839 if (error > testtol) { 1840 ierr = PetscPrintf(PETSC_COMM_WORLD, 1841 "Test failed with error norm %g\n", 1842 (double)error); CHKERRQ(ierr); 1843 } 1844 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1845 } 1846 1847 // Clean up libCEED 1848 CeedVectorDestroy(&qdata); 1849 CeedVectorDestroy(&user->qceed); 1850 CeedVectorDestroy(&user->qdotceed); 1851 CeedVectorDestroy(&user->gceed); 1852 CeedVectorDestroy(&xcorners); 1853 CeedBasisDestroy(&basisq); 1854 CeedBasisDestroy(&basisx); 1855 CeedBasisDestroy(&basisxc); 1856 CeedElemRestrictionDestroy(&restrictq); 1857 CeedElemRestrictionDestroy(&restrictx); 1858 CeedElemRestrictionDestroy(&restrictqdi); 1859 CeedQFunctionDestroy(&qf_setupVol); 1860 CeedQFunctionDestroy(&qf_ics); 1861 CeedQFunctionDestroy(&qf_rhsVol); 1862 CeedQFunctionDestroy(&qf_ifunctionVol); 1863 CeedQFunctionContextDestroy(&ctxSetup); 1864 CeedQFunctionContextDestroy(&ctxNS); 1865 CeedQFunctionContextDestroy(&ctxAdvection2d); 1866 CeedQFunctionContextDestroy(&ctxSurface); 1867 CeedOperatorDestroy(&op_setupVol); 1868 CeedOperatorDestroy(&op_ics); 1869 CeedOperatorDestroy(&user->op_rhs_vol); 1870 CeedOperatorDestroy(&user->op_ifunction_vol); 1871 CeedDestroy(&ceed); 1872 CeedBasisDestroy(&basisqSur); 1873 CeedBasisDestroy(&basisxSur); 1874 CeedBasisDestroy(&basisxcSur); 1875 CeedQFunctionDestroy(&qf_setupSur); 1876 CeedQFunctionDestroy(&qf_applySur); 1877 CeedOperatorDestroy(&user->op_rhs); 1878 CeedOperatorDestroy(&user->op_ifunction); 1879 1880 // Clean up PETSc 1881 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1882 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1883 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1884 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1885 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1886 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1887 ierr = PetscFree(units); CHKERRQ(ierr); 1888 ierr = PetscFree(user); CHKERRQ(ierr); 1889 return PetscFinalize(); 1890 } 1891