Commit 14b78ebe authored by Markus Wolff's avatar Markus Wolff
Browse files

Added additional standard test scenarios for testing the mpfa

implementation



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10706 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent db68877d
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#ifndef DUMUX_BUCKLEYLEVERETT_ANALYTICAL_HH
#define DUMUX_BUCKLEYLEVERETT_ANALYTICAL_HH
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
/**
* @file
* @brief Analytical solution of the buckley-leverett problem
* @author Markus Wolff
*/
namespace Dumux
{
/**
* \ingroup fracflow
* @brief IMplicit Pressure Explicit Saturation (IMPES) scheme for the solution of
* the Buckley-Leverett problem
*/
template<typename Scalar, typename Law>
struct CheckMaterialLaw
{
static bool isLinear()
{
return false;
}
};
template<typename Scalar>
struct CheckMaterialLaw<Scalar, LinearMaterial<Scalar> >
{
static bool isLinear()
{
return true;
}
};
template<typename Scalar>
struct CheckMaterialLaw<Scalar, EffToAbsLaw< LinearMaterial<Scalar> > >
{
static bool isLinear()
{
return true;
}
};
template<class TypeTag> class BuckleyLeverettAnalytic
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
typedef typename SpatialParams::MaterialLaw MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
enum
{
dim = GridView::dimension, dimworld = GridView::dimensionworld
};
enum
{
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
satEqIdx = Indices::satEqIdx,
};
typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > BlockVector;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef Dune::FieldVector<Scalar, dimworld> GlobalPosition;
private:
// functions needed for analytical solution
void initializeAnalytic()
{
int size = problem_.gridView().size(0);
analyticSolution_.resize(size);
analyticSolution_ = 0;
errorGlobal_.resize(size);
errorGlobal_ = 0;
errorLocal_.resize(size);
errorLocal_ = 0;
return;
}
/*!
* \brief DOC ME!
* \param materialLawParams DOC ME!
*/
void prepareAnalytic()
{
ElementIterator dummyElement = problem_.gridView().template begin<0>();
const MaterialLawParams& materialLawParams(problem_.spatialParams().materialLawParams(*dummyElement));
swr_ = materialLawParams.Swr();
snr_ = materialLawParams.Snr();
Scalar porosity = problem_.spatialParams().porosity(*dummyElement);
FluidState fluidState;
fluidState.setTemperature(problem_.temperature(*dummyElement));
fluidState.setPressure(wPhaseIdx, problem_.referencePressure(*dummyElement));
fluidState.setPressure(nPhaseIdx, problem_.referencePressure(*dummyElement));
Scalar viscosityW = FluidSystem::viscosity(fluidState, wPhaseIdx);
Scalar viscosityNW = FluidSystem::viscosity(fluidState, nPhaseIdx);
if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear() && viscosityW == viscosityNW)
{
std::pair<Scalar, Scalar> entry;
entry.first = 1 - snr_;
entry.second = vTot_ / porosity;
frontParams_.push_back(entry);
}
else
{
Scalar sw0 = swr_;
Scalar fw0 = MaterialLaw::krw(materialLawParams, sw0)/viscosityW;
fw0 /= (fw0 + MaterialLaw::krn(materialLawParams, sw0)/viscosityNW);
Scalar sw1 = sw0 + deltaS_;
Scalar fw1 = MaterialLaw::krw(materialLawParams, sw1)/viscosityW;
fw1 /= (fw1 + MaterialLaw::krn(materialLawParams, sw1)/viscosityNW);
Scalar tangentSlopeOld = (fw1 - fw0)/(sw1 - sw0);
sw1 += deltaS_;
fw1 = MaterialLaw::krw(materialLawParams, sw1)/viscosityW;
fw1 /= (fw1 + MaterialLaw::krn(materialLawParams, sw1)/viscosityNW);
Scalar tangentSlopeNew = (fw1 - fw0)/(sw1 - sw0);
while (tangentSlopeNew >= tangentSlopeOld && sw1 < (1.0 - snr_))
{
tangentSlopeOld = tangentSlopeNew;
sw1 += deltaS_;
fw1 = MaterialLaw::krw(materialLawParams, sw1)/viscosityW;
fw1 /= (fw1 + MaterialLaw::krn(materialLawParams, sw1)/viscosityNW);
tangentSlopeNew = (fw1 - fw0)/(sw1 - sw0);
}
sw0 = sw1 - deltaS_;
fw0 = MaterialLaw::krw(materialLawParams, sw0)/viscosityW;
fw0 /= (fw0 + MaterialLaw::krn(materialLawParams, sw0)/viscosityNW);
Scalar sw2 = sw1 + deltaS_;
Scalar fw2 = MaterialLaw::krw(materialLawParams, sw2)/viscosityW;
fw2 /= (fw2 + MaterialLaw::krn(materialLawParams, sw2)/viscosityNW);
while (sw1 <= (1.0 - snr_))
{
std::pair<Scalar, Scalar> entry;
entry.first = sw1;
Scalar dfwdsw = (fw2 - fw0)/(sw2 - sw0);
entry.second = vTot_ / porosity * dfwdsw;
frontParams_.push_back(entry);
sw0 = sw1;
sw1 = sw2;
fw0 = fw1;
fw1 = fw2;
sw2 += deltaS_;
fw2 = MaterialLaw::krw(materialLawParams, sw2)/viscosityW;
fw2 /= (fw2 + MaterialLaw::krn(materialLawParams, sw2)/viscosityNW);
}
}
return;
}
void calcSatError()
{
Scalar globalVolume = 0;
Scalar errorNorm = 0.0;
ElementIterator eItEnd = problem_.gridView().template end<0> ();
for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
{
// get entity
ElementPointer element = *eIt;
int index = problem_.variables().index(*element);
Scalar sat = problem_.variables().cellData(index).saturation(wPhaseIdx);
Scalar volume = element->geometry().volume();
Scalar error = analyticSolution_[index] - sat;
errorLocal_[index] = error;
if (sat > swr_ + 1e-6)
{
globalVolume += volume;
errorNorm += (volume*volume*error*error);
}
}
if (globalVolume > 0.0 && errorNorm > 0.0)
{
errorNorm = std::sqrt(errorNorm)/globalVolume;
errorGlobal_ = errorNorm;
}
else
{
errorGlobal_ = 0;
}
return;
}
void updateExSol()
{
Scalar time = problem_.timeManager().time() + problem_.timeManager().timeStepSize();
// position of the fluid front
Scalar xMax = frontParams_[0].second * time;
// iterate over vertices and get analytic saturation solution
ElementIterator eItEnd = problem_.gridView().template end<0> ();
for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
{
// get global coordinate of cell center
GlobalPosition globalPos = eIt->geometry().center();
int index = problem_.variables().index(*eIt);
if (globalPos[0] > xMax)
analyticSolution_[index] = swr_;
else
{
int size = frontParams_.size();
if (size == 1)
{
analyticSolution_[index] = frontParams_[0].first;
}
else
{
// find x_f next to global coordinate of the vertex
int xnext = 0;
for (int i = 0; i < size-1; i++)
{
Scalar x = frontParams_[i].second * time;
Scalar xMinus = frontParams_[i+1].second * time;
if (globalPos[0] <= x && globalPos[0] > xMinus)
{
analyticSolution_[index] = frontParams_[i].first - (frontParams_[i].first - frontParams_[i+1].first)/(x - xMinus) * (x - globalPos[0]);
break;
}
}
}
}
}
// call function to calculate the saturation error
calcSatError();
return;
}
public:
void calculateAnalyticSolution()
{
initializeAnalytic();
updateExSol();
}
BlockVector AnalyticSolution() const
{
return analyticSolution_;
}
//Write saturation and pressure into file
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
int size = problem_.gridView().size(0);
BlockVector *analyticSolution = writer.allocateManagedBuffer (size);
BlockVector *errorGlobal = writer.allocateManagedBuffer (size);
BlockVector *errorLocal = writer.allocateManagedBuffer (size);
*analyticSolution = analyticSolution_;
*errorGlobal = errorGlobal_;
*errorLocal = errorLocal_;
writer.attachCellData(*analyticSolution, "saturation (exact solution)");
writer.attachCellData(*errorGlobal, "global error");
writer.attachCellData(*errorLocal, "local error");
return;
}
//! Construct an IMPES object.
BuckleyLeverettAnalytic(Problem& problem) :
problem_(problem), analyticSolution_(0), errorGlobal_(0), errorLocal_(0), frontParams_(0), deltaS_(1e-3)
{}
void initialize(Scalar vTot)
{
vTot_ = vTot;
initializeAnalytic();
prepareAnalytic();
}
private:
Problem& problem_;
BlockVector analyticSolution_;
BlockVector errorGlobal_;
BlockVector errorLocal_;
std::vector<std::pair<Scalar, Scalar> > frontParams_;
const Scalar deltaS_;
Scalar swr_;
Scalar snr_;
Scalar vTot_;
int tangentIdx_;
};
}
#endif
DGF
Vertex
0 0
0 3
11.25 0
17.25 3
12.75 0
18.75 3
30 0
30 3
0 0.3
0 0.6
0 0.9
0 1.2
0 1.5
0 1.8
0 2.1
0 2.4
0 2.7
1.725 3
3.45 3
5.175 3
6.9 3
8.625 3
10.35 3
12.075 3
13.8 3
15.525 3
16.65 2.7
16.05 2.4
15.45 2.1
14.85 1.8
14.25 1.5
13.65 1.2
13.05 0.9
12.45 0.6
11.85 0.3
10.125 0
9 0
7.875 0
6.75 0
5.625 0
4.5 0
3.375 0
2.25 0
1.125 0
18.15 2.7
17.55 2.4
16.95 2.1
16.35 1.8
15.75 1.5
15.15 1.2
14.55 0.9
13.95 0.6
13.35 0.3
19.875 3
21 3
22.125 3
23.25 3
24.375 3
25.5 3
26.625 3
27.75 3
28.875 3
30 2.7
30 2.4
30 2.1
30 1.8
30 1.5
30 1.2
30 0.9
30 0.6
30 0.3
28.275 0
26.55 0
24.825 0
23.1 0
21.375 0
19.65 0
17.925 0
16.2 0
14.475 0
14.985 2.7
13.32 2.7
11.655 2.7
9.99 2.7
8.325 2.7
6.66 2.7
4.995 2.7
3.33 2.7
1.665 2.7
14.445 2.4
12.84 2.4
11.235 2.4
9.63 2.4
8.025 2.4
6.42 2.4
4.815 2.4
3.21 2.4
1.605 2.4
13.905 2.1
12.36 2.1
10.815 2.1
9.27 2.1
7.725 2.1
6.18 2.1
4.635 2.1
3.09 2.1
1.545 2.1
13.365 1.8
11.88 1.8
10.395 1.8
8.91 1.8
7.425 1.8
5.94 1.8
4.455 1.8
2.97 1.8
1.485 1.8
12.825 1.5
11.4 1.5
9.975 1.5
8.55 1.5
7.125 1.5
5.7 1.5
4.275 1.5
2.85 1.5
1.425 1.5
12.285 1.2
10.92 1.2
9.555 1.2
8.19 1.2
6.825 1.2
5.46 1.2
4.095 1.2
2.73 1.2
1.365 1.2
11.745 0.9
10.44 0.9
9.135 0.9
7.83 0.9
6.525 0.9
5.22 0.9
3.915 0.9
2.61 0.9
1.305 0.9
11.205 0.6
9.96 0.6
8.715 0.6
7.47 0.6
6.225 0.6
4.98 0.6
3.735 0.6
2.49 0.6
1.245 0.6
10.665 0.3
9.48 0.3
8.295 0.3
7.11 0.3
5.925 0.3
4.74 0.3
3.555 0.3
2.37 0.3
1.185 0.3
28.815 2.7
27.63 2.7
26.445 2.7
25.26 2.7
24.075 2.7
22.89 2.7
21.705 2.7
20.52 2.7
19.335 2.7
28.755 2.4
27.51 2.4
26.265 2.4
25.02 2.4
23.775 2.4
22.53 2.4
21.285 2.4
20.04 2.4
18.795 2.4
28.695 2.1
27.39 2.1
26.085 2.1
24.78 2.1
23.475 2.1
22.17 2.1
20.865 2.1
19.56 2.1
18.255 2.1
28.635 1.8
27.27 1.8
25.905 1.8
24.54 1.8
23.175 1.8
21.81 1.8
20.445 1.8
19.08 1.8
17.715 1.8
28.575 1.5
27.15 1.5
25.725 1.5
24.3 1.5
22.875 1.5
21.45 1.5
20.025 1.5
18.6 1.5
17.175 1.5
28.515 1.2
27.03 1.2
25.545 1.2
24.06 1.2
22.575 1.2
21.09 1.2
19.605 1.2
18.12 1.2
16.635 1.2
28.455 0.9
26.91 0.9
25.365 0.9
23.82 0.9
22.275 0.9
20.73 0.9
19.185 0.9
17.64 0.9
16.095 0.9
28.395 0.6
26.79 0.6
25.185 0.6
23.58 0.6
21.975 0.6
20.37 0.6
18.765 0.6
17.16 0.6
15.555 0.6
28.335 0.3
26.67 0.3
25.005 0.3
23.34 0.3
21.675 0.3
20.01 0.3
18.345 0.3
16.68 0.3
15.015 0.3
#
Cube
map 1 3 2 0
25 3 26 80
24 25 80 81
23 24 81 82
22 23 82 83
21 22 83 84
20 21 84 85
19 20 85 86
18 19 86 87
17 18 87 88
1 17 88 16
80 26 27 89
81 80 89 90
82 81 90 91
83 82 91 92