3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
sequential/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 *****************************************************************************/
19#ifndef DUMUX_FVPRESSURE_HH
20#define DUMUX_FVPRESSURE_HH
21
22// dumux environment
23#include <type_traits>
24#include <dumux/common/math.hh>
27#include <map>
33namespace Dumux
34{
36
48template<class TypeTag> class FVPressure
49{
50 //the model implementation
52
57
58 // using declarations to abbreviate several dune classes...
59 using Element = typename GridView::Traits::template Codim<0>::Entity;
60 using Intersection = typename GridView::Intersection;
61
62 // the typenames used for the stiffness matrix and solution vector
66
68 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
69
71
72protected:
73
78 using EntryType = Dune::FieldVector<Scalar, 2>;
79
81
86 enum
87 {
88 rhs = 1,
89 matrix = 0
90
91 };
92
93 enum
94 {
95 pressEqIdx = Indices::pressureEqIdx,
96 };
97
99 void initializeMatrix();
102
103
114 void assemble(bool first);
115
117 void solve();
118
120 PressureSolution& pressure()
121 { return pressure_;}
122
124 const PressureSolution& pressure() const
125 { return pressure_;}
126
128 //result in better convergence of the linear solver!
130 {
131 for (const auto& element : elements(problem_.gridView()))
132 {
133 PrimaryVariables initValues;
134 problem_.initial(initValues, element);
135
136 int eIdxGlobal = problem_.variables().index(element);
137
138 pressure_[eIdxGlobal] = initValues[pressEqIdx];
139 }
140 }
141
142public:
152 void getSource(EntryType& entry, const Element& element, const CellData& cellData, const bool first);
153
163 void getStorage(EntryType& entry, const Element& element, const CellData& cellData, const bool first);
164
174 void getFlux(EntryType& entry, const Intersection& intersection, const CellData& cellData, const bool first);
175
186 const Intersection& intersection, const CellData& cellData, const bool first);
187
194 const Scalar pressure(const int eIdxGlobal) const
195 { return pressure_[eIdxGlobal];}
196
198 const Matrix& globalMatrix()
199 {
200 return A_;
201 }
202
204 const RHSVector& rightHandSide()
205 {
206 return f_;
207 }
208
214 {
215 int size = problem_.gridView().size(0);//resize to make sure the final grid size (after the problem was completely built) is used!
216 A_.setSize(size, size);
217 A_.setBuildMode(Matrix::random);
218 f_.resize(size);
219 pressure_.resize(size);
221 asImp_().initializeMatrix();// initialize sparse matrix
222 }
223
228 void update()
229 {
230 asImp_().assemble(false); Dune::dinfo << "pressure calculation"<< std::endl;
231 solve();
232
233 return;
234 }
235
237 {
238 DUNE_THROW(Dune::NotImplemented,"Velocity calculation not implemented in pressure model!");
239 }
240
242 {} //empty function for the case the velocity is calculated in the transport model
243
251 void serializeEntity(std::ostream &outstream, const Element &element)
252 {
253 int eIdxGlobal = problem_.variables().index(element);
254 outstream << pressure_[eIdxGlobal][0];
255 }
256
264 void deserializeEntity(std::istream &instream, const Element &element)
265 {
266 int eIdxGlobal = problem_.variables().index(element);
267 instream >> pressure_[eIdxGlobal][0];
268 }
269
279 void setFixPressureAtIndex(Scalar pressure, int eIdxGlobal)
280 {
281 fixPressure_.insert(std::make_pair(eIdxGlobal, pressure));
282 }
283
290 void unsetFixPressureAtIndex(int eIdxGlobal)
291 {
292 fixPressure_.erase(eIdxGlobal);
293 }
294
296 {
297 fixPressure_.clear();
298 }
299
304 FVPressure(Problem& problem) :
305 problem_(problem)
306 {}
307
308private:
310 Implementation &asImp_()
311 { return *static_cast<Implementation *>(this);}
312
314 const Implementation &asImp_() const
315 { return *static_cast<const Implementation *>(this);}
316
317 Problem& problem_;
318
319 PressureSolution pressure_;
320
321 std::string solverName_;
322 std::string preconditionerName_;
323protected:
324 Matrix A_;
325 RHSVector f_;
326private:
327 std::map<int, Scalar> fixPressure_;
328};
329
331template<class TypeTag>
333{
334 initializeMatrixRowSize();
335 A_.endrowsizes();
336 initializeMatrixIndices();
337 A_.endindices();
338}
339
341template<class TypeTag>
343{
344 // determine matrix row sizes
345 for (const auto& element : elements(problem_.gridView()))
346 {
347 // cell index
348 int eIdxGlobalI = problem_.variables().index(element);
349
350 // initialize row size
351 int rowSize = 1;
352
353 // run through all intersections with neighbors
354 for (const auto& intersection : intersections(problem_.gridView(), element))
355 {
356 if (intersection.neighbor())
357 rowSize++;
358 }
359 A_.setrowsize(eIdxGlobalI, rowSize);
360 }
361}
362
364template<class TypeTag>
366{
367 // determine position of matrix entries
368 for (const auto& element : elements(problem_.gridView()))
369 {
370 // cell index
371 int eIdxGlobalI = problem_.variables().index(element);
372
373 // add diagonal index
374 A_.addindex(eIdxGlobalI, eIdxGlobalI);
375
376 // run through all intersections with neighbors
377 for (const auto& intersection : intersections(problem_.gridView(), element))
378 if (intersection.neighbor())
379 {
380 // access neighbor
381 int eIdxGlobalJ = problem_.variables().index(intersection.outside());
382
383 // add off diagonal index
384 A_.addindex(eIdxGlobalI, eIdxGlobalJ);
385 }
386 }
387}
388
389
400template<class TypeTag>
402{
403 // initialization: set matrix A_ to zero
404 A_ = 0;
405 f_ = 0;
406
407 for (const auto& element : elements(problem_.gridView()))
408 {
409 // get the global index of the cell
410 int eIdxGlobalI = problem_.variables().index(element);
411
412 // assemble interior element contributions
413 if (element.partitionType() == Dune::InteriorEntity)
414 {
415 // get the cell data
416 CellData& cellDataI = problem_.variables().cellData(eIdxGlobalI);
417
418 EntryType entries(0.);
419
420 /***** source term ***********/
421 asImp_().getSource(entries, element, cellDataI, first);
422 f_[eIdxGlobalI] += entries[rhs];
423
424 /***** flux term ***********/
425 // iterate over all faces of the cell
426 for (const auto& intersection : intersections(problem_.gridView(), element))
427 {
428 /************* handle interior face *****************/
429 if (intersection.neighbor())
430 {
431 auto elementNeighbor = intersection.outside();
432
433 int eIdxGlobalJ = problem_.variables().index(elementNeighbor);
434
435 // check for hanging nodes
436 bool haveSameLevel = (element.level() == elementNeighbor.level());
437 // calculate only from one side (except for hanging nodes), but add matrix entries for both sides
438 // the last condition is needed to properly assemble in the presence
439 // of ghost elements
440 if (getPropValue<TypeTag, Properties::VisitFacesOnlyOnce>()
441 && (eIdxGlobalI > eIdxGlobalJ) && haveSameLevel
442 && elementNeighbor.partitionType() == Dune::InteriorEntity)
443 continue;
444
445 entries = 0;
446 asImp_().getFlux(entries, intersection, cellDataI, first);
447
448 //set right hand side
449 f_[eIdxGlobalI] -= entries[rhs];
450
451 // set diagonal entry
452 A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
453
454 // set off-diagonal entry
455 A_[eIdxGlobalI][eIdxGlobalJ] -= entries[matrix];
456
457 // The second condition is needed to not spoil the ghost element entries
458 if (getPropValue<TypeTag, Properties::VisitFacesOnlyOnce>()
459 && elementNeighbor.partitionType() == Dune::InteriorEntity)
460 {
461 f_[eIdxGlobalJ] += entries[rhs];
462 A_[eIdxGlobalJ][eIdxGlobalJ] += entries[matrix];
463 A_[eIdxGlobalJ][eIdxGlobalI] -= entries[matrix];
464 }
465
466 } // end neighbor
467
468 /************* boundary face ************************/
469 else
470 {
471 entries = 0;
472 asImp_().getFluxOnBoundary(entries, intersection, cellDataI, first);
473
474 //set right hand side
475 f_[eIdxGlobalI] += entries[rhs];
476 // set diagonal entry
477 A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
478 }
479 } //end interfaces loop
480 // printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
481
482 /***** storage term ***********/
483 entries = 0;
484 asImp_().getStorage(entries, element, cellDataI, first);
485 f_[eIdxGlobalI] += entries[rhs];
486 // set diagonal entry
487 A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
488 }
489 // assemble overlap and ghost element contributions
490 else
491 {
492 A_[eIdxGlobalI] = 0.0;
493 A_[eIdxGlobalI][eIdxGlobalI] = 1.0;
494 f_[eIdxGlobalI] = pressure_[eIdxGlobalI];
495 }
496 } // end grid traversal
497// printmatrix(std::cout, A_, "global stiffness matrix after assempling", "row", 11,3);
498// printvector(std::cout, f_, "right hand side", "row", 10);
499}
500
501// forward declaration
502template<class T>
504
505namespace Detail {
506template<class T> struct isParallelAMGBackend : public std::false_type {};
507template<class T>
508struct isParallelAMGBackend<Dumux::AMGBiCGSTABBackend<T>> : public std::true_type {};
509} // end namespace Detail
510
511template<class Solver, class Problem>
512typename std::enable_if_t<!Detail::isParallelAMGBackend<Solver>::value, Solver>
513getSolver(const Problem& problem)
514{
515 return Solver();
516}
517
518template<class Solver, class Problem>
519typename std::enable_if_t<Detail::isParallelAMGBackend<Solver>::value, Solver>
520getSolver(const Problem& problem)
521{
522 return Solver(problem.gridView(), problem.model().dofMapper());
523}
524
526template<class TypeTag>
528{
530
531 int verboseLevelSolver = getParam<int>("LinearSolver.Verbosity", 0);
532
533 if (verboseLevelSolver)
534 std::cout << __FILE__ << ": solve for pressure" << std::endl;
535
536 //set a fixed pressure for a certain cell
537 if (fixPressure_.size() > 0)
538 {
539 for (auto it = fixPressure_.begin(); it != fixPressure_.end(); ++it)
540 {
541 A_[it->first] = 0;
542 A_[it->first][it->first] = 1;
543 f_[it->first] = it->second;
544 }
545 }
546
547// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
548// printvector(std::cout, f_, "right hand side", "row", 10, 1, 3);
549
550 auto solver = getSolver<Solver>(problem_);
551 bool converged = solver.solve(A_, pressure_, f_);
552
553 if (!converged)
554 DUNE_THROW(Dumux::NumericalProblem, "Pressure solver did not converge!");
555
556// printvector(std::cout, pressure_, "pressure", "row", 200, 1, 3);
557}
558
559} //end namespace Dumux
560#endif
Some exceptions thrown in DuMux
Define some often used mathematical functions.
Definition: adapt.hh:29
std::enable_if_t<!Detail::isParallelAMGBackend< Solver >::value, Solver > getSolver(const Problem &problem)
Definition: sequential/cellcentered/pressure.hh:513
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
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
A linear solver based on the ISTL AMG preconditioner and the ISTL BiCGSTAB solver.
Definition: amgbackend.hh:50
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:49
void getSource(EntryType &entry, const Element &element, const CellData &cellData, const bool first)
Function which calculates the source entry.
@ pressEqIdx
Definition: sequential/cellcentered/pressure.hh:95
void initializeMatrixIndices()
Initialize the global matrix of the system of equations to solve.
Definition: sequential/cellcentered/pressure.hh:365
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.
FVPressure(Problem &problem)
Constructs a FVPressure object.
Definition: sequential/cellcentered/pressure.hh:304
void initialize()
Initialize pressure model.
Definition: sequential/cellcentered/pressure.hh:213
void getFlux(EntryType &entry, const Intersection &intersection, const CellData &cellData, const bool first)
Function which calculates the flux entry.
void updateVelocity()
Definition: sequential/cellcentered/pressure.hh:241
void deserializeEntity(std::istream &instream, const Element &element)
Function for deserialization of the pressure field.
Definition: sequential/cellcentered/pressure.hh:264
void resetFixPressureAtIndex()
Definition: sequential/cellcentered/pressure.hh:295
void unsetFixPressureAtIndex(int eIdxGlobal)
Reset the fixed pressure state.
Definition: sequential/cellcentered/pressure.hh:290
void initializePressure()
Initialization of the pressure solution vector: Initialization with meaningful values may.
Definition: sequential/cellcentered/pressure.hh:129
PressureSolution & pressure()
Returns the vector containing the pressure solution.
Definition: sequential/cellcentered/pressure.hh:120
const Matrix & globalMatrix()
Returns the global matrix of the last pressure solution step.
Definition: sequential/cellcentered/pressure.hh:198
void serializeEntity(std::ostream &outstream, const Element &element)
Function for serialization of the pressure field.
Definition: sequential/cellcentered/pressure.hh:251
const PressureSolution & pressure() const
Returns the vector containing the pressure solution.
Definition: sequential/cellcentered/pressure.hh:124
@ 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
Matrix A_
Global stiffnes matrix (sparse matrix which is build by the initializeMatrix() function)
Definition: sequential/cellcentered/pressure.hh:324
void getFluxOnBoundary(EntryType &entry, const Intersection &intersection, const CellData &cellData, const bool first)
Function which calculates the boundary flux entry.
const RHSVector & rightHandSide()
Returns the right hand side vector of the last pressure solution step.
Definition: sequential/cellcentered/pressure.hh:204
void solve()
Solves the global system of equations to get the spatial distribution of the pressure.
Definition: sequential/cellcentered/pressure.hh:527
const Scalar pressure(const int eIdxGlobal) const
Public access function for the primary pressure variable.
Definition: sequential/cellcentered/pressure.hh:194
void initializeMatrixRowSize()
Initialize the global matrix of the system of equations to solve.
Definition: sequential/cellcentered/pressure.hh:342
void calculateVelocity()
Definition: sequential/cellcentered/pressure.hh:236
RHSVector f_
Right hand side vector.
Definition: sequential/cellcentered/pressure.hh:325
void setFixPressureAtIndex(Scalar pressure, int eIdxGlobal)
Set a pressure to be fixed at a certain cell.
Definition: sequential/cellcentered/pressure.hh:279
void update()
Pressure update.
Definition: sequential/cellcentered/pressure.hh:228
Dune::FieldVector< Scalar, 2 > EntryType
Definition: sequential/cellcentered/pressure.hh:78
void initializeMatrix()
Initialize the global matrix of the system of equations to solve.
Definition: sequential/cellcentered/pressure.hh:332
Definition: sequential/cellcentered/pressure.hh:506
Base file for properties related to sequential IMPET algorithms.