3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
1p/sequential/diffusion/cellcentered/pressure.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_FVPRESSURE1P_HH
26#define DUMUX_FVPRESSURE1P_HH
27
28
29// dumux environment
32
33namespace Dumux {
34
56template<class TypeTag> class FVPressure1P: public FVPressure<TypeTag>
57{
59
63
65
67
70 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
72 using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
73
74 enum
75 {
76 dim = GridView::dimension, dimWorld = GridView::dimensionworld
77 };
78
79 enum
80 {
81 pressEqIdx = Indices::pressureEqIdx // only one equation!
82 };
83
84 enum
85 {
87 };
88
89 using Element = typename GridView::Traits::template Codim<0>::Entity;
90 using Intersection = typename GridView::Intersection;
91
92 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
93 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
94 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
95
96
97public:
98 // Function which calculates the source entry
99 void getSource(Dune::FieldVector<Scalar, 2>&, const Element&, const CellData&, const bool);
100 // Function which calculates the storage entry
102 void getStorage(Dune::FieldVector<Scalar, 2>& entry, const Element& element,
103 const CellData& cellData, const bool first)
104 {
105 entry = 0;
106 }
108 // Function which calculates the flux entry
109 void getFlux(Dune::FieldVector<Scalar, 2>&, const Intersection&, const CellData&, const bool);
110 // Function which calculates the boundary flux entry
111 void getFluxOnBoundary(Dune::FieldVector<Scalar, 2>&,
112 const Intersection&, const CellData&, const bool);
113
121 void initialize(bool solveTwice = true)
122 {
124 this->assemble(true);
125 this->solve();
126 if (solveTwice)
127 {
128 this->assemble(false);
129 this-> solve();
130 }
132 return;
133 }
134
140 void update()
141 {
144 }
145
150 {
151 int size = problem_.gridView().size(0);
152 for (int i = 0; i < size; i++)
153 {
154 CellData& cellData = problem_.variables().cellData(i);
155 storePressureSolution(i, cellData);
156 cellData.fluxData().resetVelocity();
157 }
158 }
159
166 void storePressureSolution(int eIdxGlobal, CellData& cellData)
167 {
168 Scalar press = this->pressure()[eIdxGlobal];
169
170 cellData.setPressure(press);
171 }
172
181 template<class MultiWriter>
182 void addOutputVtkFields(MultiWriter &writer)
183 {
184 ScalarSolutionType *pressure = writer.allocateManagedBuffer (
185 problem_.gridView().size(0));
186
187 *pressure = this->pressure();
188
189 writer.attachCellData(*pressure, "pressure");
190
191 return;
192 }
193
199 FVPressure1P(Problem& problem) :
200 ParentType(problem), problem_(problem),
201 gravity_(
202 problem.gravity())
203 {
204 auto element = *problem_.gridView().template begin<0> ();
205 Scalar temperature = problem_.temperature(element);
206 Scalar referencePress = problem_.referencePressure(element);
207
208 density_ = FluidSystem::density(temperature, referencePress);
209 viscosity_ = FluidSystem::viscosity(temperature, referencePress);
210 }
211
212private:
213 Problem& problem_;
214 const GravityVector& gravity_;
215 Scalar density_;
216 Scalar viscosity_;
217};
218
226template<class TypeTag>
227void FVPressure1P<TypeTag>::getSource(Dune::FieldVector<Scalar, 2>& entry, const Element& element
228 , const CellData& cellData, const bool first)
229{
230 // cell volume, assume linear map here
231 Scalar volume = element.geometry().volume();
232
233 // get sources from problem
234 PrimaryVariables sourcePhase(0.0);
235 problem_.source(sourcePhase, element);
236 sourcePhase /= density_;
237
238 entry[rhs] = volume * sourcePhase;
239
240 return;
241}
242
248template<class TypeTag>
249void FVPressure1P<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entry, const Intersection& intersection
250 , const CellData& cellData, const bool first)
251{
252 auto elementI = intersection.inside();
253 auto elementJ = intersection.outside();
254
255 // get global coordinates of cell centers
256 const GlobalPosition& globalPosI = elementI.geometry().center();
257 const GlobalPosition& globalPosJ = elementJ.geometry().center();
258
259 //get face normal
260 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
261
262 // get face area
263 Scalar faceArea = intersection.geometry().volume();
264
265 // distance vector between barycenters
266 GlobalPosition distVec = globalPosJ - globalPosI;
267
268 // compute distance between cell centers
269 Scalar dist = distVec.two_norm();
270
271 // compute vectorized permeabilities
272 DimMatrix meanPermeability(0);
273
274 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
275 problem_.spatialParams().intrinsicPermeability(elementJ));
276
277 Dune::FieldVector<Scalar, dim> permeability(0);
278 meanPermeability.mv(unitOuterNormal, permeability);
279
280 permeability/=viscosity_;
281
282 //calculate current matrix entry
283 entry[matrix] = ((permeability * unitOuterNormal) / dist) * faceArea;
284
285 //calculate right hand side
286 entry[rhs] = density_ * (permeability * gravity_) * faceArea;
287
288 return;
289}
290
299template<class TypeTag>
300void FVPressure1P<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& entry,
301const Intersection& intersection, const CellData& cellData, const bool first)
302{
303 auto element = intersection.inside();
304
305 // get global coordinates of cell centers
306 const GlobalPosition& globalPosI = element.geometry().center();
307
308 // center of face in global coordinates
309 const GlobalPosition& globalPosJ = intersection.geometry().center();
310
311 //get face normal
312 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
313
314 // get face area
315 Scalar faceArea = intersection.geometry().volume();
316
317 // distance vector between barycenters
318 GlobalPosition distVec = globalPosJ - globalPosI;
319
320 // compute distance between cell centers
321 Scalar dist = distVec.two_norm();
322
323 BoundaryTypes bcType;
324 problem_.boundaryTypes(bcType, intersection);
325 PrimaryVariables boundValues(0.0);
326
327 if (bcType.isDirichlet(pressEqIdx))
328 {
329 problem_.dirichlet(boundValues, intersection);
330
331 // permeability vector at boundary
332 // compute vectorized permeabilities
333 DimMatrix meanPermeability(0);
334
335 problem_.spatialParams().meanK(meanPermeability,
336 problem_.spatialParams().intrinsicPermeability(element));
337
338 Dune::FieldVector<Scalar, dim> permeability(0);
339 meanPermeability.mv(unitOuterNormal, permeability);
340
341 permeability/= viscosity_;
342
343 // get Dirichlet pressure boundary condition
344 Scalar pressBound = boundValues;
345
346 // calculate current matrix entry
347 entry[matrix] = ((permeability * unitOuterNormal) / dist) * faceArea;
348 entry[rhs] = entry[matrix] * pressBound;
349
350 // calculate right hand side
351 entry[rhs] -= density_ * (permeability * gravity_) * faceArea;
352
353 }
354 // set Neumann boundary condition
355 else if (bcType.isNeumann(pressEqIdx))
356 {
357 problem_.neumann(boundValues, intersection);
358 Scalar J = boundValues /= density_;
359
360 entry[rhs] = -(J * faceArea);
361 }
362 else
363 {
364 DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
365 }
366
367 return;
368}
369
370} // end namespace Dumux
371#endif
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
std::string temperature() noexcept
I/O name of temperature for equilibrium models.
Definition: name.hh:51
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Sequential OneP Model solving the equations for pressure and velocity separately.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:57
void initialize(bool solveTwice=true)
Initializes the pressure model.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:121
FVPressure1P(Problem &problem)
Constructs a FVPressure1P object.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:199
void storePressureSolution(int eIdxGlobal, CellData &cellData)
Stores the pressure solution of a cell.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:166
void getSource(Dune::FieldVector< Scalar, 2 > &, const Element &, const CellData &, const bool)
Function which calculates the source entry.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:227
void storePressureSolution()
Globally stores the pressure solution.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:149
void update()
Pressure update.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:140
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:182
void getFluxOnBoundary(Dune::FieldVector< Scalar, 2 > &, const Intersection &, const CellData &, const bool)
Function which calculates the flux entry at a boundary.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:300
void getFlux(Dune::FieldVector< Scalar, 2 > &, const Intersection &, const CellData &, const bool)
Function which calculates the flux entry.
Definition: 1p/sequential/diffusion/cellcentered/pressure.hh:249
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:49
void assemble(bool first)
Function which assembles the system of equations to be solved.
Definition: sequential/cellcentered/pressure.hh:401
void getStorage(EntryType &entry, const Element &element, const CellData &cellData, const bool first)
Function which calculates the storage entry.
void initialize()
Initialize pressure model.
Definition: sequential/cellcentered/pressure.hh:213
PressureSolution & pressure()
Returns the vector containing the pressure solution.
Definition: sequential/cellcentered/pressure.hh:120
@ rhs
index for the right hand side entry
Definition: sequential/cellcentered/pressure.hh:88
@ matrix
index for the global matrix entry
Definition: sequential/cellcentered/pressure.hh:89
void solve()
Solves the global system of equations to get the spatial distribution of the pressure.
Definition: sequential/cellcentered/pressure.hh:527
void update()
Pressure update.
Definition: sequential/cellcentered/pressure.hh:228
Defines the properties required for the single phase sequential model.
Finite Volume Diffusion Model.