| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | |||
| 8 | #ifndef DUMUX_EXAMPLE_SHALLOWWATER_FRICTION_PROBLEM_HH | ||
| 9 | #define DUMUX_EXAMPLE_SHALLOWWATER_FRICTION_PROBLEM_HH | ||
| 10 | |||
| 11 | // ## The problem file (`problem.hh`) | ||
| 12 | // | ||
| 13 | // This file contains the __problem class__ which defines the initial and boundary | ||
| 14 | // conditions for the shallow water flow simulation with bottom friction. | ||
| 15 | // In addition, the analytical solution is defined here. | ||
| 16 | // | ||
| 17 | // [[content]] | ||
| 18 | // | ||
| 19 | // ### Include files | ||
| 20 | // [[details]] includes | ||
| 21 | // | ||
| 22 | // The first include we need here is the `ShallowWaterProblem` class, the base | ||
| 23 | // class from which we will derive. | ||
| 24 | #include <dumux/freeflow/shallowwater/problem.hh> | ||
| 25 | // In addition, we need the boundaryflux header, which handles the flux over | ||
| 26 | // the model boundaries. | ||
| 27 | #include <dumux/freeflow/shallowwater/boundaryfluxes.hh> | ||
| 28 | // Include the `BoundaryTypes` class which specifies the boundary types set in this problem. | ||
| 29 | #include <dumux/common/boundarytypes.hh> | ||
| 30 | // Include the `NumEqVector` class which specifies a field vector with size number of equations in this problem. | ||
| 31 | #include <dumux/common/numeqvector.hh> | ||
| 32 | // [[/details]] | ||
| 33 | |||
| 34 | // | ||
| 35 | // ### The problem class | ||
| 36 | // | ||
| 37 | // We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation. | ||
| 38 | // In addition the analytical solution of the problem is calculated. | ||
| 39 | // As this is a shallow water problem, we inherit from the basic ShallowWaterProblem. | ||
| 40 | // [[codeblock]] | ||
| 41 | namespace Dumux { | ||
| 42 | |||
| 43 | template <class TypeTag> | ||
| 44 | class RoughChannelProblem : public ShallowWaterProblem<TypeTag> | ||
| 45 | { | ||
| 46 | // [[/codeblock]] | ||
| 47 | // [[details]] convenience aliases | ||
| 48 | using ParentType = ShallowWaterProblem<TypeTag>; | ||
| 49 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 50 | using BoundaryTypes = Dumux::BoundaryTypes<PrimaryVariables::size()>; | ||
| 51 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 52 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 53 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 54 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 55 | using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView; | ||
| 56 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 57 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 58 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 59 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 60 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 61 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 62 | using NeumannFluxes = NumEqVector; | ||
| 63 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 64 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 65 | // [[/details]] | ||
| 66 | |||
| 67 | // In the constructor, we retrieve all required parameters from the input file | ||
| 68 | // [[codeblock]] | ||
| 69 | public: | ||
| 70 | 1 | RoughChannelProblem(std::shared_ptr<const GridGeometry> gridGeometry) | |
| 71 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
3 | : ParentType(gridGeometry) |
| 72 | { | ||
| 73 | // We read the parameters from the params.input file. | ||
| 74 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | constManningN_ = getParam<Scalar>("Problem.ManningN"); |
| 75 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | bedSlope_ = getParam<Scalar>("Problem.BedSlope"); |
| 76 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | discharge_ = getParam<Scalar>("Problem.Discharge"); |
| 77 | // We calculate the outflow boundary condition using the Gauckler-Manning-Strickler formula. | ||
| 78 | 1 | hBoundary_ = this->gaucklerManningStrickler(discharge_,constManningN_,bedSlope_); | |
| 79 | // We initialize the analytic solution to a vector of the appropriate size filled with zeros. | ||
| 80 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | exactWaterDepth_.resize(gridGeometry->numDofs(), 0.0); |
| 81 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | exactVelocityX_.resize(gridGeometry->numDofs(), 0.0); |
| 82 | 1 | } | |
| 83 | // [[/codeblock]] | ||
| 84 | |||
| 85 | // #### Analytical Solution | ||
| 86 | // | ||
| 87 | // The analytical solution is calculated using the equation of Gauckler, Manning and Strickler. | ||
| 88 | // [[codeblock]] | ||
| 89 | |||
| 90 | // Equation of Gauckler, Manning and Strickler | ||
| 91 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | Scalar gaucklerManningStrickler(Scalar discharge, Scalar manningN, Scalar bedSlope) |
| 92 | { | ||
| 93 | using std::pow; | ||
| 94 | using std::abs; | ||
| 95 | using std::sqrt; | ||
| 96 | |||
| 97 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return pow(abs(discharge)*manningN/sqrt(bedSlope), 0.6); |
| 98 | } | ||
| 99 | |||
| 100 | // Calculate the analytical solution | ||
| 101 | 1 | void analyticalSolution() | |
| 102 | { | ||
| 103 | using std::abs; | ||
| 104 | |||
| 105 |
1/2✓ Branch 1 taken 2501 times.
✗ Branch 2 not taken.
|
5001 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 106 | { | ||
| 107 | 2500 | const Scalar h = this->gaucklerManningStrickler(discharge_,constManningN_,bedSlope_); | |
| 108 | 2500 | const Scalar u = abs(discharge_)/h; | |
| 109 | |||
| 110 | 2500 | const auto eIdx = this->gridGeometry().elementMapper().index(element); | |
| 111 | 2500 | exactWaterDepth_[eIdx] = h; | |
| 112 | 2500 | exactVelocityX_[eIdx] = u; | |
| 113 | } | ||
| 114 | 1 | } | |
| 115 | |||
| 116 | // [[exclude]] | ||
| 117 | // Getter function for the analytical solution of the water depth | ||
| 118 | 1 | const std::vector<Scalar>& getExactWaterDepth() | |
| 119 | { | ||
| 120 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return exactWaterDepth_; |
| 121 | } | ||
| 122 | |||
| 123 | // Getter function for the analytical solution of the velocity in x-direction | ||
| 124 | 1 | const std::vector<Scalar>& getExactVelocityX() | |
| 125 | { | ||
| 126 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return exactVelocityX_; |
| 127 | } | ||
| 128 | // [[/exclude]] | ||
| 129 | // [[/codeblock]] | ||
| 130 | |||
| 131 | // #### Bottom friction | ||
| 132 | // | ||
| 133 | // The bottom friction is a source term and therefore handled by the `source` function. | ||
| 134 | // [[codeblock]] | ||
| 135 | 854447 | NumEqVector source(const Element& element, | |
| 136 | const FVElementGeometry& fvGeometry, | ||
| 137 | const ElementVolumeVariables& elemVolVars, | ||
| 138 | const SubControlVolume &scv) const | ||
| 139 | { | ||
| 140 | |||
| 141 | 854447 | NumEqVector source (0.0); | |
| 142 | |||
| 143 | // Since the bed slope source term is handles within the flux computation, | ||
| 144 | // in this model the bottom friction is the only source term. | ||
| 145 | 854447 | source += bottomFrictionSource(element, fvGeometry, elemVolVars, scv); | |
| 146 | |||
| 147 | 854447 | return source; | |
| 148 | } | ||
| 149 | // [[/codeblock]] | ||
| 150 | |||
| 151 | // The calculation of the source term due to bottom friction needs the bottom shear stess. | ||
| 152 | // This is the force per area, which works between the flow and the channel bed | ||
| 153 | // (1D vector with two entries) and is calculated within the `FrictionLaw` class. | ||
| 154 | // The bottom friction causes a loss of momentum. Thus, the first entry of the `bottomFrictionSource`, | ||
| 155 | // which is related to the mass balance equation is zero. | ||
| 156 | // The second entry of the `bottomFricitonSource` corresponds to the momentum equation in x-direction | ||
| 157 | // and is therefore equal to the first, the x-component, of the `bottomShearStress`. | ||
| 158 | // Accordingly, the third entry of the `bottomFrictionSource` is equal to the second component of the `bottomShearStress`. | ||
| 159 | // [[codeblock]] | ||
| 160 | 854447 | NumEqVector bottomFrictionSource(const Element& element, | |
| 161 | const FVElementGeometry& fvGeometry, | ||
| 162 | const ElementVolumeVariables& elemVolVars, | ||
| 163 | const SubControlVolume &scv) const | ||
| 164 | { | ||
| 165 | 854447 | NumEqVector bottomFrictionSource(0.0); | |
| 166 | 854447 | const auto& volVars = elemVolVars[scv]; | |
| 167 | |||
| 168 | // bottom shear stress vector | ||
| 169 | 854447 | Dune::FieldVector<Scalar, 2> bottomShearStress = this->spatialParams().frictionLaw(element, scv).bottomShearStress(volVars); | |
| 170 | |||
| 171 | // source term due to bottom friction | ||
| 172 | 854447 | bottomFrictionSource[Indices::massBalanceIdx] = 0.0; | |
| 173 | 854447 | bottomFrictionSource[Indices::momentumXBalanceIdx] = -bottomShearStress[0] / volVars.density(); | |
| 174 | 854447 | bottomFrictionSource[Indices::momentumYBalanceIdx] = -bottomShearStress[1] / volVars.density(); | |
| 175 | |||
| 176 | 854447 | return bottomFrictionSource; | |
| 177 | } | ||
| 178 | // [[/codeblock]] | ||
| 179 | |||
| 180 | // #### Boundary conditions | ||
| 181 | // | ||
| 182 | // We use Neumann boundary conditions on the entire boundary. | ||
| 183 | // [[codeblock]] | ||
| 184 | 435239 | BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const | |
| 185 | { | ||
| 186 |
2/2✓ Branch 0 taken 1305717 times.
✓ Branch 1 taken 435239 times.
|
1740956 | BoundaryTypes bcTypes; |
| 187 | 435239 | bcTypes.setAllNeumann(); | |
| 188 | 435239 | return bcTypes; | |
| 189 | } | ||
| 190 | // [[/codeblock]] | ||
| 191 | |||
| 192 | // In the following function we implement the __Neumann boundary conditions__. | ||
| 193 | // Due to the weak imposition we calculate the flux at the boundary with a Riemann solver. | ||
| 194 | // This needs the state of a virtual cell outside of the boundary (`boundaryStateVariables`), | ||
| 195 | // which is calculated with the Riemann invariants | ||
| 196 | // (see: Yoon and Kang, "Finite Volume Model for Two-Dimensional Shallow Water Flows on Unstructured Grids"). | ||
| 197 | // [[codeblock]] | ||
| 198 | 345349 | NeumannFluxes neumann(const Element& element, | |
| 199 | const FVElementGeometry& fvGeometry, | ||
| 200 | const ElementVolumeVariables& elemVolVars, | ||
| 201 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 202 | const SubControlVolumeFace& scvf) const | ||
| 203 | { | ||
| 204 | 345349 | NeumannFluxes values(0.0); | |
| 205 | |||
| 206 | 345349 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 207 |
2/2✓ Branch 0 taken 1780 times.
✓ Branch 1 taken 343569 times.
|
345349 | const auto& insideVolVars = elemVolVars[insideScv]; |
| 208 | 345349 | const auto& nxy = scvf.unitOuterNormal(); | |
| 209 |
2/2✓ Branch 0 taken 1780 times.
✓ Branch 1 taken 343569 times.
|
345349 | const auto gravity = this->spatialParams().gravity(scvf.center()); |
| 210 | std::array<Scalar, 3> boundaryStateVariables; | ||
| 211 | |||
| 212 | // Calculate the Riemann invariants for imposed discharge at the left side. | ||
| 213 |
2/2✓ Branch 0 taken 1780 times.
✓ Branch 1 taken 343569 times.
|
345349 | if (scvf.center()[0] < this->gridGeometry().bBoxMin()[0] + eps_) |
| 214 | { | ||
| 215 | 1780 | boundaryStateVariables = ShallowWater::fixedDischargeBoundary(discharge_, | |
| 216 | insideVolVars.waterDepth(), | ||
| 217 | insideVolVars.velocity(0), | ||
| 218 | insideVolVars.velocity(1), | ||
| 219 | gravity, | ||
| 220 | nxy); | ||
| 221 | } | ||
| 222 | // Calculate the Riemann invariants for imposed water depth at the right side. | ||
| 223 |
2/2✓ Branch 0 taken 1780 times.
✓ Branch 1 taken 341789 times.
|
343569 | else if (scvf.center()[0] > this->gridGeometry().bBoxMax()[0] - eps_) |
| 224 | { | ||
| 225 | 1780 | boundaryStateVariables = ShallowWater::fixedWaterDepthBoundary(hBoundary_, | |
| 226 | insideVolVars.waterDepth(), | ||
| 227 | insideVolVars.velocity(0), | ||
| 228 | insideVolVars.velocity(1), | ||
| 229 | gravity, | ||
| 230 | nxy); | ||
| 231 | } | ||
| 232 | // Calculate the Riemann invariants for the no-flow boundary. | ||
| 233 | else | ||
| 234 | { | ||
| 235 | 341789 | boundaryStateVariables[0] = insideVolVars.waterDepth(); | |
| 236 | 341789 | boundaryStateVariables[1] = -insideVolVars.velocity(0); | |
| 237 | 341789 | boundaryStateVariables[2] = -insideVolVars.velocity(1); | |
| 238 | } | ||
| 239 | // Calculate the boundary fluxes based on a Riemann problem. | ||
| 240 | 345349 | auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(), | |
| 241 | 345349 | boundaryStateVariables[0], | |
| 242 | insideVolVars.velocity(0), | ||
| 243 | 345349 | boundaryStateVariables[1], | |
| 244 | insideVolVars.velocity(1), | ||
| 245 | 345349 | boundaryStateVariables[2], | |
| 246 | insideVolVars.bedSurface(), | ||
| 247 | insideVolVars.bedSurface(), | ||
| 248 | gravity, | ||
| 249 | nxy); | ||
| 250 | |||
| 251 | 345349 | values[Indices::massBalanceIdx] = riemannFlux[0]; | |
| 252 | 345349 | values[Indices::velocityXIdx] = riemannFlux[1]; | |
| 253 | 345349 | values[Indices::velocityYIdx] = riemannFlux[2]; | |
| 254 | |||
| 255 | 345349 | return values; | |
| 256 | } | ||
| 257 | // [[/codeblock]] | ||
| 258 | |||
| 259 | // #### Initial conditions | ||
| 260 | // | ||
| 261 | // We specify the initial conditions for the primary variables (water depth, velocity in y-direction | ||
| 262 | // and velocity in x-direction). In this example constant initial conditions are used. Therefore the | ||
| 263 | //argument `globalPos` is not needed. If you want to impose spatial variable initial conditions, | ||
| 264 | // you have to use the `globalPos` argument. | ||
| 265 | // [[codeblock]] | ||
| 266 | 2500 | PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const | |
| 267 | { | ||
| 268 | 2500 | PrimaryVariables initialValues(0.0); | |
| 269 | initialValues[Indices::waterdepthIdx] = 1.0; | ||
| 270 | initialValues[Indices::velocityXIdx] = 0.0; | ||
| 271 | initialValues[Indices::velocityYIdx] = 0.0; | ||
| 272 | return initialValues; | ||
| 273 | } | ||
| 274 | // [[/codeblock]] | ||
| 275 | |||
| 276 | // [[details]] private variables | ||
| 277 | // [[codeblock]] | ||
| 278 | private: | ||
| 279 | // variables for the analytic solution. | ||
| 280 | std::vector<Scalar> exactWaterDepth_; | ||
| 281 | std::vector<Scalar> exactVelocityX_; | ||
| 282 | // constant friction value (an analytic solution is only available for const friction). | ||
| 283 | Scalar constManningN_; | ||
| 284 | // The constant channel bed slope. | ||
| 285 | Scalar bedSlope_; | ||
| 286 | // The water depth at the outflow boundary. | ||
| 287 | Scalar hBoundary_; | ||
| 288 | // The discharge at the inflow boundary. | ||
| 289 | Scalar discharge_; | ||
| 290 | // We assign a private global variable for the epsilon: | ||
| 291 | static constexpr Scalar eps_ = 1.0e-6; | ||
| 292 | |||
| 293 | }; // end class definition RoughChannelProblem | ||
| 294 | } // end namespace Dumux | ||
| 295 | // [[/codeblock]] | ||
| 296 | // [[/details]] | ||
| 297 | // [[/content]] | ||
| 298 | #endif | ||
| 299 |