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