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