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