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