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 eulercontext; 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->eulercontext->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) { 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, ctxSetup); 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 eulercontext; 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; 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, &eulercontext); 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); CHKERRQ(ierr); 1323 1324 // Refine DM for high-order viz 1325 dmviz = NULL; 1326 interpviz = NULL; 1327 if (viz_refine) { 1328 DM dmhierarchy[viz_refine+1]; 1329 1330 ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); CHKERRQ(ierr); 1331 dmhierarchy[0] = dm; 1332 for (PetscInt i = 0, d = degree; i < viz_refine; i++) { 1333 Mat interp_next; 1334 1335 ierr = DMRefine(dmhierarchy[i], MPI_COMM_NULL, &dmhierarchy[i+1]); 1336 CHKERRQ(ierr); 1337 ierr = DMClearDS(dmhierarchy[i+1]); CHKERRQ(ierr); 1338 ierr = DMClearFields(dmhierarchy[i+1]); CHKERRQ(ierr); 1339 ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); 1340 d = (d + 1) / 2; 1341 if (i + 1 == viz_refine) d = 1; 1342 ierr = SetUpDM(dmhierarchy[i+1], problem, d, &bc, &ctxSetupData); 1343 CHKERRQ(ierr); 1344 ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], 1345 &interp_next, NULL); CHKERRQ(ierr); 1346 if (!i) interpviz = interp_next; 1347 else { 1348 Mat C; 1349 ierr = MatMatMult(interp_next, interpviz, MAT_INITIAL_MATRIX, 1350 PETSC_DECIDE, &C); CHKERRQ(ierr); 1351 ierr = MatDestroy(&interp_next); CHKERRQ(ierr); 1352 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1353 interpviz = C; 1354 } 1355 } 1356 for (PetscInt i=1; i<viz_refine; i++) { 1357 ierr = DMDestroy(&dmhierarchy[i]); CHKERRQ(ierr); 1358 } 1359 dmviz = dmhierarchy[viz_refine]; 1360 } 1361 ierr = DMCreateGlobalVector(dm, &Q); CHKERRQ(ierr); 1362 ierr = DMGetLocalVector(dm, &Qloc); CHKERRQ(ierr); 1363 ierr = VecGetSize(Qloc, &lnodes); CHKERRQ(ierr); 1364 lnodes /= ncompq; 1365 1366 // Initialize CEED 1367 CeedInit(ceedresource, &ceed); 1368 // Set memtype 1369 CeedMemType memtypebackend; 1370 CeedGetPreferredMemType(ceed, &memtypebackend); 1371 // Check memtype compatibility 1372 if (!setmemtyperequest) 1373 memtyperequested = memtypebackend; 1374 else if (!petschavecuda && memtyperequested == CEED_MEM_DEVICE) 1375 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, 1376 "PETSc was not built with CUDA. " 1377 "Requested MemType CEED_MEM_DEVICE is not supported.", NULL); 1378 1379 // Set number of 1D nodes and quadrature points 1380 numP = degree + 1; 1381 numQ = numP + qextra; 1382 1383 // Print summary 1384 if (!test) { 1385 CeedInt gdofs, odofs; 1386 int comm_size; 1387 char box_faces_str[PETSC_MAX_PATH_LEN] = "NONE"; 1388 ierr = VecGetSize(Q, &gdofs); CHKERRQ(ierr); 1389 ierr = VecGetLocalSize(Q, &odofs); CHKERRQ(ierr); 1390 gnodes = gdofs/ncompq; 1391 ierr = MPI_Comm_size(comm, &comm_size); CHKERRQ(ierr); 1392 ierr = PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, 1393 sizeof(box_faces_str), NULL); CHKERRQ(ierr); 1394 const char *usedresource; 1395 CeedGetResource(ceed, &usedresource); 1396 1397 ierr = PetscPrintf(comm, 1398 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 1399 " rank(s) : %d\n" 1400 " Problem:\n" 1401 " Problem Name : %s\n" 1402 " Stabilization : %s\n" 1403 " PETSc:\n" 1404 " Box Faces : %s\n" 1405 " libCEED:\n" 1406 " libCEED Backend : %s\n" 1407 " libCEED Backend MemType : %s\n" 1408 " libCEED User Requested MemType : %s\n" 1409 " Mesh:\n" 1410 " Number of 1D Basis Nodes (P) : %d\n" 1411 " Number of 1D Quadrature Points (Q) : %d\n" 1412 " Global DoFs : %D\n" 1413 " Owned DoFs : %D\n" 1414 " DoFs per node : %D\n" 1415 " Global nodes : %D\n" 1416 " Owned nodes : %D\n", 1417 comm_size, problemTypes[problemChoice], 1418 StabilizationTypes[stab], box_faces_str, usedresource, 1419 CeedMemTypes[memtypebackend], 1420 (setmemtyperequest) ? 1421 CeedMemTypes[memtyperequested] : "none", 1422 numP, numQ, gdofs, odofs, ncompq, gnodes, lnodes); 1423 CHKERRQ(ierr); 1424 } 1425 1426 // Set up global mass vector 1427 ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr); 1428 1429 // Set up libCEED 1430 // CEED Bases 1431 CeedInit(ceedresource, &ceed); 1432 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompq, numP, numQ, CEED_GAUSS, 1433 &basisq); 1434 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numQ, CEED_GAUSS, 1435 &basisx); 1436 CeedBasisCreateTensorH1Lagrange(ceed, dim, ncompx, 2, numP, 1437 CEED_GAUSS_LOBATTO, &basisxc); 1438 ierr = DMGetCoordinateDM(dm, &dmcoord); CHKERRQ(ierr); 1439 ierr = DMPlexSetClosurePermutationTensor(dmcoord, PETSC_DETERMINE, NULL); 1440 CHKERRQ(ierr); 1441 1442 // CEED Restrictions 1443 ierr = GetRestrictionForDomain(ceed, dm, 0, 0, 0, numP, numQ, 1444 qdatasizeVol, &restrictq, &restrictx, 1445 &restrictqdi); CHKERRQ(ierr); 1446 1447 ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); 1448 ierr = CreateVectorFromPetscVec(ceed, Xloc, &xcorners); CHKERRQ(ierr); 1449 1450 // Create the CEED vectors that will be needed in setup 1451 CeedInt NqptsVol; 1452 CeedBasisGetNumQuadraturePoints(basisq, &NqptsVol); 1453 CeedElemRestrictionGetNumElements(restrictq, &localNelemVol); 1454 CeedVectorCreate(ceed, qdatasizeVol*localNelemVol*NqptsVol, &qdata); 1455 CeedElemRestrictionCreateVector(restrictq, &q0ceed, NULL); 1456 1457 // Create the Q-function that builds the quadrature data for the NS operator 1458 CeedQFunctionCreateInterior(ceed, 1, problem->setupVol, problem->setupVol_loc, 1459 &qf_setupVol); 1460 CeedQFunctionAddInput(qf_setupVol, "dx", ncompx*dim, CEED_EVAL_GRAD); 1461 CeedQFunctionAddInput(qf_setupVol, "weight", 1, CEED_EVAL_WEIGHT); 1462 CeedQFunctionAddOutput(qf_setupVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1463 1464 // Create the Q-function that sets the ICs of the operator 1465 CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); 1466 CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); 1467 CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); 1468 1469 qf_rhsVol = NULL; 1470 if (problem->applyVol_rhs) { // Create the Q-function that defines the action of the RHS operator 1471 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_rhs, 1472 problem->applyVol_rhs_loc, &qf_rhsVol); 1473 CeedQFunctionAddInput(qf_rhsVol, "q", ncompq, CEED_EVAL_INTERP); 1474 CeedQFunctionAddInput(qf_rhsVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1475 CeedQFunctionAddInput(qf_rhsVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1476 CeedQFunctionAddInput(qf_rhsVol, "x", ncompx, CEED_EVAL_INTERP); 1477 CeedQFunctionAddOutput(qf_rhsVol, "v", ncompq, CEED_EVAL_INTERP); 1478 CeedQFunctionAddOutput(qf_rhsVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1479 } 1480 1481 qf_ifunctionVol = NULL; 1482 if (problem->applyVol_ifunction) { // Create the Q-function that defines the action of the IFunction 1483 CeedQFunctionCreateInterior(ceed, 1, problem->applyVol_ifunction, 1484 problem->applyVol_ifunction_loc, &qf_ifunctionVol); 1485 CeedQFunctionAddInput(qf_ifunctionVol, "q", ncompq, CEED_EVAL_INTERP); 1486 CeedQFunctionAddInput(qf_ifunctionVol, "dq", ncompq*dim, CEED_EVAL_GRAD); 1487 CeedQFunctionAddInput(qf_ifunctionVol, "qdot", ncompq, CEED_EVAL_INTERP); 1488 CeedQFunctionAddInput(qf_ifunctionVol, "qdata", qdatasizeVol, CEED_EVAL_NONE); 1489 CeedQFunctionAddInput(qf_ifunctionVol, "x", ncompx, CEED_EVAL_INTERP); 1490 CeedQFunctionAddOutput(qf_ifunctionVol, "v", ncompq, CEED_EVAL_INTERP); 1491 CeedQFunctionAddOutput(qf_ifunctionVol, "dv", ncompq*dim, CEED_EVAL_GRAD); 1492 } 1493 1494 // Create the operator that builds the quadrature data for the NS operator 1495 CeedOperatorCreate(ceed, qf_setupVol, NULL, NULL, &op_setupVol); 1496 CeedOperatorSetField(op_setupVol, "dx", restrictx, basisx, CEED_VECTOR_ACTIVE); 1497 CeedOperatorSetField(op_setupVol, "weight", CEED_ELEMRESTRICTION_NONE, 1498 basisx, CEED_VECTOR_NONE); 1499 CeedOperatorSetField(op_setupVol, "qdata", restrictqdi, 1500 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1501 1502 // Create the operator that sets the ICs 1503 CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics); 1504 CeedOperatorSetField(op_ics, "x", restrictx, basisxc, CEED_VECTOR_ACTIVE); 1505 CeedOperatorSetField(op_ics, "q0", restrictq, 1506 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 1507 1508 CeedElemRestrictionCreateVector(restrictq, &user->qceed, NULL); 1509 CeedElemRestrictionCreateVector(restrictq, &user->qdotceed, NULL); 1510 CeedElemRestrictionCreateVector(restrictq, &user->gceed, NULL); 1511 1512 if (qf_rhsVol) { // Create the RHS physics operator 1513 CeedOperator op; 1514 CeedOperatorCreate(ceed, qf_rhsVol, NULL, NULL, &op); 1515 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1516 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1517 CeedOperatorSetField(op, "qdata", restrictqdi, 1518 CEED_BASIS_COLLOCATED, qdata); 1519 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1520 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1521 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1522 user->op_rhs_vol = op; 1523 } 1524 1525 if (qf_ifunctionVol) { // Create the IFunction operator 1526 CeedOperator op; 1527 CeedOperatorCreate(ceed, qf_ifunctionVol, NULL, NULL, &op); 1528 CeedOperatorSetField(op, "q", restrictq, basisq, CEED_VECTOR_ACTIVE); 1529 CeedOperatorSetField(op, "dq", restrictq, basisq, CEED_VECTOR_ACTIVE); 1530 CeedOperatorSetField(op, "qdot", restrictq, basisq, user->qdotceed); 1531 CeedOperatorSetField(op, "qdata", restrictqdi, 1532 CEED_BASIS_COLLOCATED, qdata); 1533 CeedOperatorSetField(op, "x", restrictx, basisx, xcorners); 1534 CeedOperatorSetField(op, "v", restrictq, basisq, CEED_VECTOR_ACTIVE); 1535 CeedOperatorSetField(op, "dv", restrictq, basisq, CEED_VECTOR_ACTIVE); 1536 user->op_ifunction_vol = op; 1537 } 1538 1539 // Set up CEED for the boundaries 1540 CeedInt height = 1; 1541 CeedInt dimSur = dim - height; 1542 CeedInt numP_Sur = degree + 1; 1543 CeedInt numQ_Sur = numP_Sur + qextraSur; 1544 const CeedInt qdatasizeSur = problem->qdatasizeSur; 1545 CeedBasis basisxSur, basisxcSur, basisqSur; 1546 CeedInt NqptsSur; 1547 CeedQFunction qf_setupSur, qf_applySur; 1548 1549 // CEED bases for the boundaries 1550 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompq, numP_Sur, numQ_Sur, 1551 CEED_GAUSS, 1552 &basisqSur); 1553 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numQ_Sur, CEED_GAUSS, 1554 &basisxSur); 1555 CeedBasisCreateTensorH1Lagrange(ceed, dimSur, ncompx, 2, numP_Sur, 1556 CEED_GAUSS_LOBATTO, &basisxcSur); 1557 CeedBasisGetNumQuadraturePoints(basisqSur, &NqptsSur); 1558 1559 // Create the Q-function that builds the quadrature data for the Surface operator 1560 CeedQFunctionCreateInterior(ceed, 1, problem->setupSur, problem->setupSur_loc, 1561 &qf_setupSur); 1562 CeedQFunctionAddInput(qf_setupSur, "dx", ncompx*dimSur, CEED_EVAL_GRAD); 1563 CeedQFunctionAddInput(qf_setupSur, "weight", 1, CEED_EVAL_WEIGHT); 1564 CeedQFunctionAddOutput(qf_setupSur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1565 1566 // Creat Q-Function for Boundaries 1567 // -- Defined for Advection(2d) test cases 1568 qf_applySur = NULL; 1569 if (problem->applySur) { 1570 CeedQFunctionCreateInterior(ceed, 1, problem->applySur, 1571 problem->applySur_loc, &qf_applySur); 1572 CeedQFunctionAddInput(qf_applySur, "q", ncompq, CEED_EVAL_INTERP); 1573 CeedQFunctionAddInput(qf_applySur, "qdataSur", qdatasizeSur, CEED_EVAL_NONE); 1574 CeedQFunctionAddInput(qf_applySur, "x", ncompx, CEED_EVAL_INTERP); 1575 CeedQFunctionAddOutput(qf_applySur, "v", ncompq, CEED_EVAL_INTERP); 1576 } 1577 1578 // Create CEED Operator for the whole domain 1579 if (!implicit) 1580 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1581 user->op_rhs_vol, qf_applySur, 1582 qf_setupSur, height, numP_Sur, numQ_Sur, 1583 qdatasizeSur, NqptsSur, basisxSur, 1584 basisqSur, &user->op_rhs); 1585 CHKERRQ(ierr); 1586 if (implicit) 1587 ierr = CreateOperatorForDomain(ceed, dm, &bc, problemChoice, wind_type, 1588 user->op_ifunction_vol, qf_applySur, 1589 qf_setupSur, height, numP_Sur, numQ_Sur, 1590 qdatasizeSur, NqptsSur, basisxSur, 1591 basisqSur, &user->op_ifunction); 1592 CHKERRQ(ierr); 1593 // Set up contex for QFunctions 1594 CeedQFunctionContextCreate(ceed, &ctxSetup); 1595 CeedQFunctionContextSetData(ctxSetup, CEED_MEM_HOST, CEED_USE_POINTER, 1596 sizeof ctxSetupData, &ctxSetupData); 1597 CeedQFunctionSetContext(qf_ics, ctxSetup); 1598 1599 CeedScalar ctxNSData[8] = {lambda, mu, k, cv, cp, g, Rd, , user->currentTime}; 1600 CeedQFunctionContextCreate(ceed, &ctxNS); 1601 CeedQFunctionContextSetData(ctxNS, CEED_MEM_HOST, CEED_USE_POINTER, 1602 sizeof ctxNSData, &ctxNSData); 1603 1604 struct Advection2dContext_ ctxAdvection2d = { 1605 .CtauS = CtauS, 1606 .strong_form = strong_form, 1607 .stabilization = stab, 1608 }; 1609 CeedQFunctionContextCreate(ceed, &ctxAdvection2d); 1610 CeedQFunctionContextSetData(ctxAdvection2d, CEED_MEM_HOST, CEED_USE_POINTER, 1611 sizeof ctxAdvection2dData, &ctxAdvection2dData); 1612 1613 struct SurfaceContext_ ctxSurfaceData = { 1614 .E_wind = E_wind, 1615 .strong_form = strong_form, 1616 .implicit = implicit, 1617 }; 1618 struct EulerContext_ eulerctx = { 1619 .currentTime = 0., 1620 }; 1621 CeedQFunctionContextCreate(ceed, &ctxSurface); 1622 CeedQFunctionContextSetData(ctxSurface, CEED_MEM_HOST, CEED_USE_POINTER, 1623 sizeof ctxSurfaceData, &ctxSurfaceData); 1624 1625 switch (problemChoice) { 1626 case NS_DENSITY_CURRENT: 1627 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxNS, sizeof ctxNS); 1628 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxNS, 1629 sizeof ctxNS); 1630 break; 1631 case NS_ADVECTION: 1632 case NS_ADVECTION2D: 1633 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &ctxAdvection2d, 1634 sizeof ctxAdvection2d); 1635 if (qf_ifunctionVol) CeedQFunctionSetContext(qf_ifunctionVol, &ctxAdvection2d, 1636 sizeof ctxAdvection2d); 1637 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &ctxSurface, 1638 sizeof ctxSurface); 1639 case NS_EULER_VORTEX: 1640 if (qf_rhsVol) CeedQFunctionSetContext(qf_rhsVol, &eulercontext, 1641 sizeof eulercontext); 1642 if (qf_applySur) CeedQFunctionSetContext(qf_applySur, &eulercontext, 1643 sizeof eulercontext); 1644 } 1645 1646 // Set up PETSc context 1647 // Set up units structure 1648 units->meter = meter; 1649 units->kilogram = kilogram; 1650 units->second = second; 1651 units->Kelvin = Kelvin; 1652 units->Pascal = Pascal; 1653 units->JperkgK = JperkgK; 1654 units->mpersquareds = mpersquareds; 1655 units->WpermK = WpermK; 1656 units->kgpercubicm = kgpercubicm; 1657 units->kgpersquaredms = kgpersquaredms; 1658 units->Joulepercubicm = Joulepercubicm; 1659 units->Joule = Joule; 1660 1661 // Set up user structure 1662 user->comm = comm; 1663 user->outputfreq = outputfreq; 1664 user->contsteps = contsteps; 1665 user->units = units; 1666 user->dm = dm; 1667 user->dmviz = dmviz; 1668 user->interpviz = interpviz; 1669 user->ceed = ceed; 1670 1671 // Calculate qdata and ICs 1672 // Set up state global and local vectors 1673 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 1674 1675 ierr = VectorPlacePetscVec(q0ceed, Qloc); CHKERRQ(ierr); 1676 1677 // Apply Setup Ceed Operators 1678 ierr = VectorPlacePetscVec(xcorners, Xloc); CHKERRQ(ierr); 1679 CeedOperatorApply(op_setupVol, xcorners, qdata, CEED_REQUEST_IMMEDIATE); 1680 ierr = ComputeLumpedMassMatrix(ceed, dm, restrictq, basisq, restrictqdi, qdata, 1681 user->M); CHKERRQ(ierr); 1682 1683 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qloc, Q, restrictq, 1684 ctxSetup, 0.0); CHKERRQ(ierr); 1685 if (1) { // Record boundary values from initial condition and override DMPlexInsertBoundaryValues() 1686 // We use this for the main simulation DM because the reference DMPlexInsertBoundaryValues() is very slow. If we 1687 // disable this, we should still get the same results due to the problem->bc function, but with potentially much 1688 // slower execution. 1689 Vec Qbc; 1690 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1691 ierr = VecCopy(Qloc, Qbc); CHKERRQ(ierr); 1692 ierr = VecZeroEntries(Qloc); CHKERRQ(ierr); 1693 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Qloc); CHKERRQ(ierr); 1694 ierr = VecAXPY(Qbc, -1., Qloc); CHKERRQ(ierr); 1695 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 1696 ierr = PetscObjectComposeFunction((PetscObject)dm, 1697 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 1698 CHKERRQ(ierr); 1699 } 1700 1701 MPI_Comm_rank(comm, &rank); 1702 if (!rank) {ierr = PetscMkdir(user->outputdir); CHKERRQ(ierr);} 1703 // Gather initial Q values 1704 // In case of continuation of simulation, set up initial values from binary file 1705 if (contsteps) { // continue from existent solution 1706 PetscViewer viewer; 1707 char filepath[PETSC_MAX_PATH_LEN]; 1708 // Read input 1709 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-solution.bin", 1710 user->outputdir); 1711 CHKERRQ(ierr); 1712 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1713 CHKERRQ(ierr); 1714 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 1715 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1716 } 1717 ierr = DMRestoreLocalVector(dm, &Qloc); CHKERRQ(ierr); 1718 1719 // Create and setup TS 1720 ierr = TSCreate(comm, &ts); CHKERRQ(ierr); 1721 ierr = TSSetDM(ts, dm); CHKERRQ(ierr); 1722 if (implicit) { 1723 ierr = TSSetType(ts, TSBDF); CHKERRQ(ierr); 1724 if (user->op_ifunction) { 1725 ierr = TSSetIFunction(ts, NULL, IFunction_NS, &user); CHKERRQ(ierr); 1726 } else { // Implicit integrators can fall back to using an RHSFunction 1727 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1728 } 1729 } else { 1730 if (!user->op_rhs) SETERRQ(comm, PETSC_ERR_ARG_NULL, 1731 "Problem does not provide RHSFunction"); 1732 ierr = TSSetType(ts, TSRK); CHKERRQ(ierr); 1733 ierr = TSRKSetType(ts, TSRK5F); CHKERRQ(ierr); 1734 ierr = TSSetRHSFunction(ts, NULL, RHS_NS, &user); CHKERRQ(ierr); 1735 } 1736 ierr = TSSetMaxTime(ts, 500. * units->second); CHKERRQ(ierr); 1737 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr); 1738 ierr = TSSetTimeStep(ts, 1.e-2 * units->second); CHKERRQ(ierr); 1739 if (test) {ierr = TSSetMaxSteps(ts, 10); CHKERRQ(ierr);} 1740 ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); 1741 ierr = TSAdaptSetStepLimits(adapt, 1742 1.e-12 * units->second, 1743 1.e2 * units->second); CHKERRQ(ierr); 1744 ierr = TSSetFromOptions(ts); CHKERRQ(ierr); 1745 if (!contsteps) { // print initial condition 1746 if (!test) { 1747 ierr = TSMonitor_NS(ts, 0, 0., Q, user); CHKERRQ(ierr); 1748 } 1749 } else { // continue from time of last output 1750 PetscReal time; 1751 PetscInt count; 1752 PetscViewer viewer; 1753 char filepath[PETSC_MAX_PATH_LEN]; 1754 ierr = PetscSNPrintf(filepath, sizeof filepath, "%s/ns-time.bin", 1755 user->outputdir); CHKERRQ(ierr); 1756 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1757 CHKERRQ(ierr); 1758 ierr = PetscViewerBinaryRead(viewer, &time, 1, &count, PETSC_REAL); 1759 CHKERRQ(ierr); 1760 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1761 ierr = TSSetTime(ts, time * user->units->second); CHKERRQ(ierr); 1762 } 1763 if (!test) { 1764 ierr = TSMonitorSet(ts, TSMonitor_NS, user, NULL); CHKERRQ(ierr); 1765 } 1766 1767 // Solve 1768 start = MPI_Wtime(); 1769 ierr = PetscBarrier((PetscObject)ts); CHKERRQ(ierr); 1770 ierr = TSSolve(ts, Q); CHKERRQ(ierr); 1771 cpu_time_used = MPI_Wtime() - start; 1772 ierr = TSGetSolveTime(ts, &ftime); CHKERRQ(ierr); 1773 ierr = MPI_Allreduce(MPI_IN_PLACE, &cpu_time_used, 1, MPI_DOUBLE, MPI_MIN, 1774 comm); CHKERRQ(ierr); 1775 if (!test) { 1776 ierr = PetscPrintf(PETSC_COMM_WORLD, 1777 "Time taken for solution (sec): %g\n", 1778 (double)cpu_time_used); CHKERRQ(ierr); 1779 } 1780 1781 // Get error 1782 if (problem->non_zero_time && !test) { 1783 Vec Qexact, Qexactloc; 1784 PetscReal norm; 1785 ierr = DMCreateGlobalVector(dm, &Qexact); CHKERRQ(ierr); 1786 ierr = DMGetLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1787 ierr = VecGetSize(Qexactloc, &lnodes); CHKERRQ(ierr); 1788 1789 ierr = ICs_FixMultiplicity(op_ics, xcorners, q0ceed, dm, Qexactloc, Qexact, 1790 restrictq, ctxSetup, ftime); CHKERRQ(ierr); 1791 1792 ierr = VecAXPY(Q, -1.0, Qexact); CHKERRQ(ierr); 1793 ierr = VecNorm(Q, NORM_MAX, &norm); CHKERRQ(ierr); 1794 CeedVectorDestroy(&q0ceed); 1795 ierr = PetscPrintf(PETSC_COMM_WORLD, 1796 "Max Error: %g\n", 1797 (double)norm); CHKERRQ(ierr); 1798 // Clean up vectors 1799 ierr = DMRestoreLocalVector(dm, &Qexactloc); CHKERRQ(ierr); 1800 ierr = VecDestroy(&Qexact); CHKERRQ(ierr); 1801 } 1802 1803 // Output Statistics 1804 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 1805 if (!test) { 1806 ierr = PetscPrintf(PETSC_COMM_WORLD, 1807 "Time integrator took %D time steps to reach final time %g\n", 1808 steps, (double)ftime); CHKERRQ(ierr); 1809 } 1810 // Output numerical values from command line 1811 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 1812 1813 // Compare reference solution values with current test run for CI 1814 if (test) { 1815 PetscViewer viewer; 1816 // Read reference file 1817 Vec Qref; 1818 PetscReal error, Qrefnorm; 1819 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 1820 ierr = PetscViewerBinaryOpen(comm, filepath, FILE_MODE_READ, &viewer); 1821 CHKERRQ(ierr); 1822 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 1823 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1824 1825 // Compute error with respect to reference solution 1826 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 1827 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 1828 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 1829 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 1830 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 1831 // Check error 1832 if (error > testtol) { 1833 ierr = PetscPrintf(PETSC_COMM_WORLD, 1834 "Test failed with error norm %g\n", 1835 (double)error); CHKERRQ(ierr); 1836 } 1837 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 1838 } 1839 1840 // Clean up libCEED 1841 CeedVectorDestroy(&qdata); 1842 CeedVectorDestroy(&user->qceed); 1843 CeedVectorDestroy(&user->qdotceed); 1844 CeedVectorDestroy(&user->gceed); 1845 CeedVectorDestroy(&xcorners); 1846 CeedBasisDestroy(&basisq); 1847 CeedBasisDestroy(&basisx); 1848 CeedBasisDestroy(&basisxc); 1849 CeedElemRestrictionDestroy(&restrictq); 1850 CeedElemRestrictionDestroy(&restrictx); 1851 CeedElemRestrictionDestroy(&restrictqdi); 1852 CeedQFunctionDestroy(&qf_setupVol); 1853 CeedQFunctionDestroy(&qf_ics); 1854 CeedQFunctionDestroy(&qf_rhsVol); 1855 CeedQFunctionDestroy(&qf_ifunctionVol); 1856 CeedQFunctionContextDestroy(&ctxSetup); 1857 CeedQFunctionContextDestroy(&ctxNS); 1858 CeedQFunctionContextDestroy(&ctxAdvection2d); 1859 CeedQFunctionContextDestroy(&ctxSurface); 1860 CeedOperatorDestroy(&op_setupVol); 1861 CeedOperatorDestroy(&op_ics); 1862 CeedOperatorDestroy(&user->op_rhs_vol); 1863 CeedOperatorDestroy(&user->op_ifunction_vol); 1864 CeedDestroy(&ceed); 1865 CeedBasisDestroy(&basisqSur); 1866 CeedBasisDestroy(&basisxSur); 1867 CeedBasisDestroy(&basisxcSur); 1868 CeedQFunctionDestroy(&qf_setupSur); 1869 CeedQFunctionDestroy(&qf_applySur); 1870 CeedOperatorDestroy(&user->op_rhs); 1871 CeedOperatorDestroy(&user->op_ifunction); 1872 1873 // Clean up PETSc 1874 ierr = VecDestroy(&Q); CHKERRQ(ierr); 1875 ierr = VecDestroy(&user->M); CHKERRQ(ierr); 1876 ierr = MatDestroy(&interpviz); CHKERRQ(ierr); 1877 ierr = DMDestroy(&dmviz); CHKERRQ(ierr); 1878 ierr = TSDestroy(&ts); CHKERRQ(ierr); 1879 ierr = DMDestroy(&dm); CHKERRQ(ierr); 1880 ierr = PetscFree(units); CHKERRQ(ierr); 1881 ierr = PetscFree(user); CHKERRQ(ierr); 1882 ierr = PetscFree(eulercontext); CHKERRQ(ierr); 1883 return PetscFinalize(); 1884 } 1885