diff --git a/test/decoupled/2p/buckleyleverettanalyticsolution.hh b/test/decoupled/2p/buckleyleverettanalyticsolution.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8d8dfb7b0ee869f74c8e71e393030b9397ead5da
--- /dev/null
+++ b/test/decoupled/2p/buckleyleverettanalyticsolution.hh
@@ -0,0 +1,353 @@
+// -*- 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
diff --git a/test/decoupled/2p/grids/test_mpfa2p.dgf b/test/decoupled/2p/grids/test_mpfa2p.dgf
new file mode 100644
index 0000000000000000000000000000000000000000..3ec5a022cdd2bc1193484c6922ec6320d0459f93
--- /dev/null
+++ b/test/decoupled/2p/grids/test_mpfa2p.dgf
@@ -0,0 +1,461 @@
+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
+84 83 92 93
+85 84 93 94
+86 85 94 95
+87 86 95 96
+88 87 96 97
+16 88 97 15
+89 27 28 98
+90 89 98 99
+91 90 99 100
+92 91 100 101
+93 92 101 102
+94 93 102 103
+95 94 103 104
+96 95 104 105
+97 96 105 106
+15 97 106 14
+98 28 29 107
+99 98 107 108
+100 99 108 109
+101 100 109 110
+102 101 110 111
+103 102 111 112
+104 103 112 113
+105 104 113 114
+106 105 114 115
+14 106 115 13
+107 29 30 116
+108 107 116 117
+109 108 117 118
+110 109 118 119
+111 110 119 120
+112 111 120 121
+113 112 121 122
+114 113 122 123
+115 114 123 124
+13 115 124 12
+116 30 31 125
+117 116 125 126
+118 117 126 127
+119 118 127 128
+120 119 128 129
+121 120 129 130
+122 121 130 131
+123 122 131 132
+124 123 132 133
+12 124 133 11
+125 31 32 134
+126 125 134 135
+127 126 135 136
+128 127 136 137
+129 128 137 138
+130 129 138 139
+131 130 139 140
+132 131 140 141
+133 132 141 142
+11 133 142 10
+134 32 33 143
+135 134 143 144
+136 135 144 145
+137 136 145 146
+138 137 146 147
+139 138 147 148
+140 139 148 149
+141 140 149 150
+142 141 150 151
+10 142 151 9
+143 33 34 152
+144 143 152 153
+145 144 153 154
+146 145 154 155
+147 146 155 156
+148 147 156 157
+149 148 157 158
+150 149 158 159
+151 150 159 160
+9 151 160 8
+152 34 2 35
+153 152 35 36
+154 153 36 37
+155 154 37 38
+156 155 38 39
+157 156 39 40
+158 157 40 41
+159 158 41 42
+160 159 42 43
+8 160 43 0
+3 5 44 26
+26 44 45 27
+27 45 46 28
+28 46 47 29
+29 47 48 30
+30 48 49 31
+31 49 50 32
+32 50 51 33
+33 51 52 34
+34 52 4 2
+61 7 62 161
+60 61 161 162
+59 60 162 163
+58 59 163 164
+57 58 164 165
+56 57 165 166
+55 56 166 167
+54 55 167 168
+53 54 168 169
+5 53 169 44
+161 62 63 170
+162 161 170 171
+163 162 171 172
+164 163 172 173
+165 164 173 174
+166 165 174 175
+167 166 175 176
+168 167 176 177
+169 168 177 178
+44 169 178 45
+170 63 64 179
+171 170 179 180
+172 171 180 181
+173 172 181 182
+174 173 182 183
+175 174 183 184
+176 175 184 185
+177 176 185 186
+178 177 186 187
+45 178 187 46
+179 64 65 188
+180 179 188 189
+181 180 189 190
+182 181 190 191
+183 182 191 192
+184 183 192 193
+185 184 193 194
+186 185 194 195
+187 186 195 196
+46 187 196 47
+188 65 66 197
+189 188 197 198
+190 189 198 199
+191 190 199 200
+192 191 200 201
+193 192 201 202
+194 193 202 203
+195 194 203 204
+196 195 204 205
+47 196 205 48
+197 66 67 206
+198 197 206 207
+199 198 207 208
+200 199 208 209
+201 200 209 210
+202 201 210 211
+203 202 211 212
+204 203 212 213
+205 204 213 214
+48 205 214 49
+206 67 68 215
+207 206 215 216
+208 207 216 217
+209 208 217 218
+210 209 218 219
+211 210 219 220
+212 211 220 221
+213 212 221 222
+214 213 222 223
+49 214 223 50
+215 68 69 224
+216 215 224 225
+217 216 225 226
+218 217 226 227
+219 218 227 228
+220 219 228 229
+221 220 229 230
+222 221 230 231
+223 222 231 232
+50 223 232 51
+224 69 70 233
+225 224 233 234
+226 225 234 235
+227 226 235 236
+228 227 236 237
+229 228 237 238
+230 229 238 239
+231 230 239 240
+232 231 240 241
+51 232 241 52
+233 70 6 71
+234 233 71 72
+235 234 72 73
+236 235 73 74
+237 236 74 75
+238 237 75 76
+239 238 76 77
+240 239 77 78
+241 240 78 79
+52 241 79 4
+#
+Boundarydomain
+default 1
+#
diff --git a/test/decoupled/2p/mcwhorteranalyticsolution.hh b/test/decoupled/2p/mcwhorteranalyticsolution.hh
new file mode 100644
index 0000000000000000000000000000000000000000..91278161fcf3e9163cac4f1b3dc30a649086099a
--- /dev/null
+++ b/test/decoupled/2p/mcwhorteranalyticsolution.hh
@@ -0,0 +1,427 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ *
+ * 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_MCWHORTER_ANALYTIC_HH
+#define DUMUX_MCWHORTER_ANALYTIC_HH
+
+/**
+ * @file
+ * @brief Analytic solution of
+ * the McWhorter problem
+ * @author Markus Wolff, Anneli Schöniger
+ */
+
+namespace Dumux
+{
+/**
+ * @brief Analytic solution of
+ * the McWhorter problem
+ *
+ * for naming of variables see "An Improved Semi-Analytical Solution for Verification
+ * of Numerical Models of Two-Phase Flow in Porous Media"
+ * (R. Fucik, J. Mikyska, M. Benes, T. H. Illangasekare; 2007)
+ */
+
+template<typename TypeTag>
+class McWhorterAnalytic
+{
+ 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(TypeTag, SolutionTypes)::PrimaryVariables PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum
+ {
+ dimworld = GridView::dimensionworld //dim = GridView::dimension,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+// pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx//,
+// pressEqIdx = Indices::pressEqIdx,
+// 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;
+ }
+
+ 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 prepareAnalytic()
+ {
+ ElementIterator dummyElement = problem_.gridView().template begin<0>();
+ const MaterialLawParams& materialLawParams(problem_.spatialParams().materialLawParams(*dummyElement));
+
+ swr_ = materialLawParams.Swr();
+ snr_ = materialLawParams.Snr();
+ porosity_ = problem_.spatialParams().porosity(*dummyElement);
+ permeability_ = problem_.spatialParams().intrinsicPermeability(*dummyElement)[0][0];
+ PrimaryVariables initVec;
+ problem_.initial(initVec, *dummyElement);
+ sInit_ = initVec[saturationIdx];
+ Scalar s0 =(1 - snr_ - swr_);
+// std::cerr<<"\n Swr, Snr: "<< swr_ << snr_ << "\n";
+// std::cerr<<"\n sInit "<< sInit_<< "\n";
+
+ // h_= (s0 - sInit_)/intervalNum_;
+ h_= (s0)/intervalNum_;
+// std::cout<<"h_= "<<h_<<std::endl;
+
+ // define saturation range for analytic solution
+ satVec_= swr_;
+ for (int i=1; i<pointNum_; i++)
+ {
+ satVec_[i]=satVec_[i-1]+h_;
+ }
+ 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);
+
+ // get fractional flow function vector
+ for (int i=0; i<pointNum_; i++)
+ {
+ fractionalW_[i] = MaterialLaw::krw(materialLawParams, satVec_[i])/viscosityW;
+ fractionalW_[i] /= (fractionalW_[i] + MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW);
+ }
+
+ // get capillary pressure derivatives
+ dpcdsw_=0;
+
+ for (int i=0; i<pointNum_; i++)
+ {
+ dpcdsw_[i] = MaterialLaw::dpC_dSw(materialLawParams, satVec_[i]);
+ }
+// std::cout<<"dpcdsw = "<<dpcdsw_<<std::endl;
+
+ // set initial fW
+ if (sInit_ == 0)
+ fInit_=0;
+ else
+ fInit_=fractionalW_[0];
+
+ fractionalW_[0]=0;
+
+ // normalize fW
+ // with r_ = qt/q0
+ // qt: total volume flux, q0: displacing phase flux at boundary
+ // --> r_ = 1 for unidirectional displacement; r_ = 0 for impermeable boundary
+ for (int i=0; i<pointNum_; i++)
+ {
+ fn_[i]= r_ * (fractionalW_[i] - fInit_)/ (1 - r_ * fInit_);
+ }
+
+ // std::cout<<"r_ = "<<r_<<std::endl;
+ // std::cout<<"fn_ = "<<fn_<<std::endl;
+
+ // diffusivity function
+ for (int i=0; i<pointNum_; i++)
+ {d_[i] = fractionalW_[i]*(-dpcdsw_[i])*(MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW)*permeability_;
+ }
+
+ // std::cout<<"fractionalW_ = "<<fractionalW_<<std::endl;
+ // std::cout<<"permeability_ = "<<permeability_<<std::endl;
+ // std::cout<<"d_ = "<<d_<<std::endl;
+
+
+ // gk_: part of fractional flow function
+ // initial guess for gk_
+ for (int i=0; i<pointNum_; i++)
+ {
+ gk_[i] = d_[i]/(1-fn_[i]);
+ }
+
+ gk_[0] = 0;
+
+// std::cout<<"gk_ = "<<gk_<<std::endl;
+
+ return;
+ }
+
+ void updateExSol()
+ {
+ Scalar time = problem_.timeManager().time() + problem_.timeManager().timeStepSize();
+
+ // with displacing phase flux at boundary q0 = A * 1/sqrt(t)
+ // Akm1, Ak: successive approximations of A
+ double Ak = 0;
+ double Akm1 = 0;
+ double diff = 1e100;
+
+ // partial numerical integrals a_, b_
+ a_=0, b_=0;
+ fp_=0;
+
+ // approximation of integral I
+ double I0 = 0;
+ double Ii = 0;
+
+ int k = 0;
+
+ while (diff> tolAnalytic_)
+ {
+ k++;
+// std::cout<<" k = "<<k<<std::endl;
+ if (k> 50000)
+ {
+ std::cout<<"Analytic solution: Too many iterations!"<<std::endl;
+ break;
+ }
+
+ Akm1=Ak;
+ I0=0;
+ for (int i=0; i<intervalNum_; i++)
+ {
+ a_[i] = 0.5 * h_ * sInit_ *(gk_[i] + gk_[i+1])+ pow(h_, 2) / 6* ((3* i + 1) * gk_[i]
+ + (3 * i + 2) * gk_[i+1]);
+ b_[i] = 0.5 * h_ * (gk_[i] + gk_[i+1]);
+ I0 += (a_[i] - sInit_ * b_[i]);
+ }
+ // std::cout<<" I0 = "<<I0<<std::endl;
+
+ gk_[0]=0;
+ for (int i=1; i<pointNum_; i++)
+ {
+ Ii=0;
+ for (int j = i; j<intervalNum_; j++)
+ Ii += (a_[j] - satVec_[i] * b_[j]);
+ //gk_[i] = d_[i] + gk_[i]*(fn_[i] + Ii/I0); // method A
+ gk_[i] = (d_[i] + gk_[i]*fn_[i])/(1 - Ii/I0); // method B
+ }
+
+ // with f(sInit) = 0: relationship between A and sInit
+ Ak = pow((0.5*porosity_/pow((1 - fInit_), 2)*I0), 0.5);
+ diff=fabs(Ak - Akm1);
+ // std::cout<<"diff = "<<diff<<std::endl;
+ }
+
+ // std::cout<<" b_ = "<<b_<<std::endl;
+ // std::cout<<" Ak = "<<Ak<<std::endl;
+
+
+ // fp_: first derivative of f
+ for (int i = 0; i<pointNum_; i++)
+ {
+ for (int j = i; j<intervalNum_; j++)
+ fp_[i] += b_[j]/I0;
+ }
+
+ // std::cout<<" fp_ = "<<fp_<<std::endl;
+ // std::cout<<fInit_<<std::endl;
+
+ for (int i = 0; i<pointNum_; i++)
+ {
+ xf_[i]= 2 * Ak * (1 - fInit_ * r_)/ porosity_ * fp_[i]* pow(time, 0.5);
+ }
+
+// std::cout<<" xf_ = "<<xf_<<std::endl;
+
+ // iterate over vertices and get analytical 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
+ const GlobalPosition& globalPos = eIt->geometry().center();
+
+// std::cout<<"globalPos = "<<globalPos[0]<<", x0 = "<<xf_[0]<<"\n";
+
+ // find index of current vertex
+ int index = problem_.variables().index(*eIt);
+
+ // find x_f next to global coordinate of the vertex
+ int xnext = 0;
+ for (int i=intervalNum_; i>0; i--)
+ {
+ if (globalPos[0] <= xf_[i])
+ {
+ xnext = i;
+ break;
+ }
+ }
+
+ // account for the area not yet reached by the front
+ if (globalPos[0]> xf_[0])
+ {
+ analyticSolution_[index] = sInit_;
+ continue;
+ }
+
+ if (globalPos[0] <= xf_[0])
+ {
+ analyticSolution_[index] = satVec_[xnext];
+ continue;
+ }
+
+// std::cout<<"Analytical = "<<satVec_[xnext]<<std::endl;
+ }
+
+ // 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;
+ }
+
+
+ void initialize()
+ {
+ initializeAnalytic();
+ prepareAnalytic();
+ }
+
+ McWhorterAnalytic(Problem& problem, Scalar tolAnalytic = 1e-14) :
+ problem_(problem), analyticSolution_(0), errorGlobal_(0), errorLocal_(0), tolAnalytic_(tolAnalytic)
+ {}
+
+private:
+ Problem& problem_;
+
+ BlockVector analyticSolution_;
+ BlockVector errorGlobal_;
+ BlockVector errorLocal_;
+
+ Scalar tolAnalytic_;
+ Scalar r_;
+
+ Scalar swr_;
+ Scalar snr_;
+ Scalar porosity_;
+ Scalar sInit_;
+ Scalar permeability_;
+ enum
+ { intervalNum_ = 1000, pointNum_ = intervalNum_+1};
+ Dune::FieldVector<Scalar, pointNum_> satVec_;
+ Dune::FieldVector<Scalar,pointNum_> fractionalW_;
+ Dune::FieldVector<Scalar, pointNum_> dpcdsw_;
+ Dune::FieldVector<Scalar, pointNum_> fn_;
+ Dune::FieldVector<Scalar, pointNum_> d_;
+ Dune::FieldVector<Scalar, pointNum_> gk_;
+ Dune::FieldVector<Scalar, pointNum_> xf_;
+ Scalar fInit_;
+ Scalar h_;
+ Dune::FieldVector<Scalar,intervalNum_> a_;
+ Dune::FieldVector<Scalar,intervalNum_> b_;
+ Dune::FieldVector<Scalar,pointNum_> fp_;
+
+};
+}
+#endif
diff --git a/test/decoupled/2p/test_mpfa2p.cc b/test/decoupled/2p/test_mpfa2p.cc
index 3dbf94340a48fb95aa9464b87caffffe2c003712..d18ca5f77e49f779aacd6b0bd2fb2580af4210b8 100644
--- a/test/decoupled/2p/test_mpfa2p.cc
+++ b/test/decoupled/2p/test_mpfa2p.cc
@@ -19,7 +19,7 @@
*
* \brief Test for the sequential 2p models
*/
-#define STRUCTUREDGRID 1
+#define PROBLEM 2 // 0 = Buckley-Leverett, 1 = McWhorter, 2 = 2D Lense problem
#include "config.h"
diff --git a/test/decoupled/2p/test_mpfa2p.input b/test/decoupled/2p/test_mpfa2p.input
index 1cfca93af00637adc000feee1a93f776d21a853c..72e6fc708308ab1dcbd9690fb832a3a3b9966ec8 100644
--- a/test/decoupled/2p/test_mpfa2p.input
+++ b/test/decoupled/2p/test_mpfa2p.input
@@ -9,10 +9,13 @@
############################################################
[TimeManager]
TEnd = 5e4# [s]
-DtInitial = 0# [s]
+DtInitial = 5e5# [s]
+SubTimestepVerbosity = 0
[Grid]
-#For nine-spot problem
+#File = grids/test_mpfa2p.dgf#grid for buckley-leverett or mcwhorter problem
+
+#For 2d-lense-problem
NumberOfCellsX = 10# [-] level 0 resolution in x-direction
NumberOfCellsY = 10# [-] level 0 resolution in y-direction
@@ -23,18 +26,21 @@ UpperRightY = 10# [m] dimension of the grid # 2D
# Optional arguments
############################################################
[Problem]
-EnableGravity = true
OutputInterval =0
OutputTimeInterval = 1e4
+#For 2d-lense-problem
InletWidth = 2
InjectionFlux = 0.1
-#OutputfileName =
+[Impet]
+#CFLFactor = 10.0 # Use larger time steps and
+#SubCFLFactor = 0.95 # enable local time-stepping
[Vtk]
-#OutputLevel = 1 # 0 -> only primary variables (default), 1 -> also secondary variables
+OutputLevel = 1 # 0 -> only primary variables (default), 1 -> also secondary variables
+#For 2d-lense-problem
[SpatialParams]
BackgroundPermeabilityXX = 1e-10
BackgroundPermeabilityXY = 0
@@ -72,4 +78,4 @@ CoarsenTolerance = 0.01 # threshold for coarsening criterion
EnableInitializationIndicator = true#
RefineAtDirichletBC = false
RefineAtFluxBC = true
-RefineAtSource = false
\ No newline at end of file
+RefineAtSource = false
diff --git a/test/decoupled/2p/test_mpfa2pproblem.hh b/test/decoupled/2p/test_mpfa2pproblem.hh
index 8b69b31499a693cf2528072d0cd95f67f44800ed..19b645b28876ecf324da356c4d684a8e72727a9d 100644
--- a/test/decoupled/2p/test_mpfa2pproblem.hh
+++ b/test/decoupled/2p/test_mpfa2pproblem.hh
@@ -47,6 +47,8 @@
#include<dumux/decoupled/2p/impes/gridadaptionindicator2plocal.hh>
#include "test_mpfa2pspatialparams.hh"
+#include "buckleyleverettanalyticsolution.hh"
+#include "mcwhorteranalyticsolution.hh"
namespace Dumux
{
@@ -68,8 +70,10 @@ SET_PROP(MPFATwoPTestProblem, Grid)
typedef Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming> type;
};
-// set the GridCreator property
+#if PROBLEM == 2
+//set the GridCreator property
SET_TYPE_PROP(MPFATwoPTestProblem, GridCreator, CubeGridCreator<TypeTag>);
+#endif
// Set the problem property
@@ -88,6 +92,7 @@ public:
typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
};
+#if PROBLEM == 2
// Set the non-wetting phase
SET_PROP(MPFATwoPTestProblem, NonwettingPhase)
{
@@ -96,14 +101,34 @@ private:
public:
typedef Dumux::LiquidPhase<Scalar, Dumux::DNAPL<Scalar> > type;
};
+#else
+// Set the non-wetting phase
+SET_PROP(MPFATwoPTestProblem, NonwettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
+};
+#endif
+
+#if PROBLEM == 1
+SET_INT_PROP(MPFATwoPTestProblem, Formulation, DecoupledTwoPCommonIndices::pnSw);
+#endif
+#if PROBLEM == 2
// Enable gravity
SET_BOOL_PROP(MPFATwoPTestProblem, ProblemEnableGravity, true);
+#else
+SET_BOOL_PROP(MPFATwoPTestProblem, ProblemEnableGravity, false);
+#endif
SET_TYPE_PROP(MPFATwoPTestProblem, EvalCflFluxFunction, Dumux::EvalCflFluxCoats<TypeTag>);
SET_SCALAR_PROP(MPFATwoPTestProblem, ImpetCFLFactor, 1.0);
SET_TYPE_PROP(MPFATwoPTestProblem, AdaptionIndicator, Dumux::GridAdaptionIndicator2PLocal<TypeTag>);
+SET_TYPE_PROP(MPFATwoPTestProblem, LinearSolver, Dumux::SuperLUBackend<TypeTag>);
+
NEW_TYPE_TAG(FVTwoPTestProblem, INHERITS_FROM(FVPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTestProblem));
NEW_TYPE_TAG(FVAdaptiveTwoPTestProblem, INHERITS_FROM(FVPressureTwoPAdaptive, FVTransportTwoP, IMPESTwoPAdaptive, MPFATwoPTestProblem));
NEW_TYPE_TAG(MPFAOTwoPTestProblem, INHERITS_FROM(FVMPFAOPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTestProblem));
@@ -167,7 +192,11 @@ enum
{
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
+#if PROBLEM == 1
+ pNIdx = Indices::pnIdx,
+#else
pWIdx = Indices::pwIdx,
+#endif
SwIdx = Indices::SwIdx,
eqIdxPress = Indices::pressEqIdx,
eqIdxSat = Indices::satEqIdx
@@ -184,9 +213,20 @@ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
public:
MPFATwoPTestProblem(TimeManager &timeManager,const GridView &gridView) :
ParentType(timeManager, gridView)
+#if PROBLEM != 2
+, analyticSolution_(*this)
+#endif
{
this->setGrid(GridCreator::grid());
+ int refinementFactor = 0;
+ if (ParameterTree::tree().hasKey("Grid.RefinementFactor"))
+ {
+ refinementFactor = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Grid, RefinementFactor);
+ }
+
+ this->grid().globalRefine(refinementFactor);
+
Scalar inletWidth = 1.0;
if (ParameterTree::tree().hasKey("Problem.InletWidth"))
{
@@ -221,6 +261,29 @@ ParentType(timeManager, gridView)
this->setOutputTimeInterval(outputTimeInterval);
}
+#if PROBLEM != 2
+void init()
+{
+ ParentType::init();
+
+#if PROBLEM == 0
+ Scalar vTot = 3e-6;
+ analyticSolution_.initialize(vTot);
+#endif
+#if PROBLEM == 1
+ analyticSolution_.initialize();
+#endif
+}
+
+void addOutputVtkFields()
+{
+ analyticSolution_.calculateAnalyticSolution();
+
+ ParentType::addOutputVtkFields();
+ analyticSolution_.addOutputVtkFields(this->resultWriter());
+}
+#endif
+
/*!
* \name Problem parameters
*/
@@ -282,6 +345,7 @@ void source(PrimaryVariables &values,const Element& element) const
*/
void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const
{
+#if PROBLEM == 2
if (isInlet(globalPos))
{
bcTypes.setNeumann(eqIdxPress);
@@ -296,11 +360,37 @@ void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos)
bcTypes.setDirichlet(eqIdxPress);
bcTypes.setOutflow(eqIdxSat);
}
+#elif PROBLEM == 0
+ if (globalPos[0] < eps_)
+ {
+ bcTypes.setAllDirichlet();
+ }
+ else if (globalPos[0] > this->bboxMax()[0] - eps_)
+ {
+ bcTypes.setNeumann(eqIdxPress);
+ bcTypes.setOutflow(eqIdxSat);
+ }
+ else
+ {
+ bcTypes.setAllNeumann();
+ }
+#else
+ if (globalPos[0] < eps_)
+ {
+ bcTypes.setAllDirichlet();
+ }
+ else
+ {
+ bcTypes.setAllNeumann();
+ }
+#endif
}
//! set dirichlet condition (pressure [Pa], saturation [-])
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
+ values = 0;
+#if PROBLEM == 2
Scalar pRef = referencePressureAtPos(globalPos);
Scalar temp = temperatureAtPos(globalPos);
@@ -311,25 +401,52 @@ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) c
{
values[SwIdx] = 0.0;
}
+#elif PROBLEM == 0
+ if (globalPos[0] < eps_)
+ {
+ values[SwIdx] = 0.8;
+ values[pWIdx] = 1;
+ }
+#else
+ if (globalPos[0] < eps_)
+ {
+ values[SwIdx] = 1.0;
+ values[pNIdx] = 1e5;
+ }
+#endif
}
//! set neumann condition for phases (flux, [kg/(m^2 s)])
void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
values = 0;
+#if PROBLEM == 2
if (isInlet(globalPos))
{
values[nPhaseIdx] = -inFlux_;
}
-
+#elif PROBLEM == 0
+ if (globalPos[0] > this->bboxMax()[0] - eps_)
+ {
+ values[nPhaseIdx] = 3e-3;
+ }
+#endif
}
//! return initial solution -> only saturation values have to be given!
void initialAtPos(PrimaryVariables &values,
const GlobalPosition& globalPos) const
{
+#if PROBLEM == 2
values[pWIdx] = 0;
values[SwIdx] = 1.0;
+#elif PROBLEM == 0
+ values[pWIdx] = 0;
+ values[SwIdx] = 0.2;
+#else
+ values[pNIdx] = 0;
+ values[SwIdx] = 0.0;
+#endif
}
private:
@@ -383,6 +500,12 @@ Scalar inFlux_;
GlobalPosition inletLeftCoord_;
GlobalPosition inletRightCoord_;
static constexpr Scalar eps_ = 1e-6;
+#if PROBLEM == 0
+BuckleyLeverettAnalytic<TypeTag> analyticSolution_;
+#endif
+#if PROBLEM == 1
+McWhorterAnalytic<TypeTag> analyticSolution_;
+#endif
};
} //end namespace
diff --git a/test/decoupled/2p/test_mpfa2pspatialparams.hh b/test/decoupled/2p/test_mpfa2pspatialparams.hh
index 5cd2f69c6507d9f9b3fbd60d71f79f373f8c3c43..068ef74b453b8a6d59083932982197a2d92a814b 100644
--- a/test/decoupled/2p/test_mpfa2pspatialparams.hh
+++ b/test/decoupled/2p/test_mpfa2pspatialparams.hh
@@ -94,7 +94,13 @@ public:
Scalar porosity(const Element& element) const
{
+#if PROBLEM == 0
+ return 0.2;
+#elif PROBLEM == 1
+ return 0.3;
+#else
return 0.4;
+#endif
}
// return the brooks-corey context depending on the position
@@ -114,6 +120,20 @@ public:
ParentType(gridView), permBackground_(0), permLenses_(0),
lensOneLowerLeft_(0), lensOneUpperRight_(0), lensTwoLowerLeft_(0), lensTwoUpperRight_(0), lensThreeLowerLeft_(0), lensThreeUpperRight_(0)
{
+#if PROBLEM == 0
+ // residual saturations
+ materialLawParamsBackground_.setSwr(0.2);
+ materialLawParamsBackground_.setSnr(0.2);
+
+ materialLawParamsLenses_.setSwr(0.2);
+ materialLawParamsLenses_.setSnr(0.2);
+
+ //parameters for Brooks-Corey law
+
+ // entry pressures function
+ materialLawParamsBackground_.setPe(0.);
+ materialLawParamsLenses_.setPe(0.);
+#elif PROBLEM == 1
// residual saturations
materialLawParamsBackground_.setSwr(0.);
materialLawParamsBackground_.setSnr(0.);
@@ -122,6 +142,20 @@ public:
materialLawParamsLenses_.setSnr(0.);
//parameters for Brooks-Corey law
+
+ // entry pressures function
+ materialLawParamsBackground_.setPe(5000.);
+ materialLawParamsLenses_.setPe(5000.);
+#else
+ // residual saturations
+ materialLawParamsBackground_.setSwr(0.);
+ materialLawParamsBackground_.setSnr(0.);
+
+ materialLawParamsLenses_.setSwr(0.);
+ materialLawParamsLenses_.setSnr(0.);
+
+ //parameters for Brooks-Corey law
+
// entry pressures function
materialLawParamsBackground_.setPe(0.);
if (ParameterTree::tree().hasKey("SpatialParams.BackgroundEntryPressure"))
@@ -134,16 +168,40 @@ public:
{
materialLawParamsLenses_.setPe(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LenseEntryPressure));
}
+#endif
// Brooks-Corey shape parameters
+
+#if PROBLEM == 2
materialLawParamsBackground_.setLambda(3);
+ if (ParameterTree::tree().hasKey("SpatialParams.BackgroundLambda"))
+ {
+ materialLawParamsBackground_.setLambda(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, BackgroundLambda));
+ }
materialLawParamsLenses_.setLambda(2);
+ if (ParameterTree::tree().hasKey("SpatialParams.LenseLambda"))
+ {
+ materialLawParamsLenses_.setLambda(GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LenseLambda));
+ }
+#else
+ materialLawParamsBackground_.setLambda(2);
+ materialLawParamsLenses_.setLambda(2);
+#endif
+#if PROBLEM == 0
+ permBackground_[0][0] = 1e-7;
+ permBackground_[1][1] = 1e-7;
+ permLenses_[0][0] = 1e-7;
+ permLenses_[1][1] = 1e-7;
+#else
permBackground_[0][0] = 1e-10;
permBackground_[1][1] = 1e-10;
permLenses_[0][0] = 1e-10;
permLenses_[1][1] = 1e-10;
+#endif
+
+#if PROBLEM == 2
if (ParameterTree::tree().hasKey("SpatialParams.BackgroundPermeabilityXX"))
{
permBackground_[0][0] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, BackgroundPermeabilityXX);
@@ -232,6 +290,7 @@ public:
{
lensThreeUpperRight_[1] = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, SpatialParams, LensThreeUpperRightY);
}
+#endif
}
private: