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