3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
plotmateriallaw3p.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 *****************************************************************************/
24#ifndef DUMUX_PLOT_FLUID_MATRIX_LAW_HH
25#define DUMUX_PLOT_FLUID_MATRIX_LAW_HH
26
27#include <cmath>
28#include <string>
29#include <vector>
30
31namespace Dumux {
32
33// forward declaration
34template<class Scalar> class GnuplotInterface;
35
41template<class Scalar, class MaterialLaw>
43{
44 using MaterialLawParams = typename MaterialLaw::Params;
45
46public:
48 PlotMaterialLaw(bool interaction = true)
49 : numIntervals_(1000)
50 {}
51
62 const MaterialLawParams &params,
63 Scalar lowerSat = 0.0,
64 Scalar upperSat = 1.0,
65 std::string curveTitle = "")
66 {
67 addpcgw(gnuplot, params, lowerSat, upperSat, curveTitle);
68 addpcnw(gnuplot, params, lowerSat, upperSat, curveTitle);
69 addpcgn(gnuplot, params, lowerSat, upperSat, curveTitle);
70 }
71
82 const MaterialLawParams &params,
83 Scalar lowerSat = 0.0,
84 Scalar upperSat = 1.0,
85 std::string curveTitle = "")
86 {
87 std::vector<Scalar> sw;
88 std::vector<Scalar> pc;
89 Scalar satInterval = upperSat - lowerSat;
90 Scalar pcMin = 0.0;
91 Scalar pcMax = -1e100;
92
93 Scalar swTemp, pcTemp = 0.0;
94 for (int i = 0; i <= numIntervals_; i++)
95 {
96 swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_);
97 pcTemp = MaterialLaw::pcgw(params, swTemp);
98 using std::max;
99 using std::min;
100 if (checkValues_(swTemp, pcTemp))
101 {
102 sw.push_back(swTemp);
103 pc.push_back(pcTemp);
104 pcMin = min(pcMin, pcTemp);
105 pcMax = max(pcMax, pcTemp);
106 }
107 }
108
109 gnuplot.setXlabel("wetting phase saturation [-]");
110 gnuplot.setYlabel("capillary pressure [Pa]");
111 gnuplot.addDataSetToPlot(sw, pc, curveTitle + "_pcgw-Sw");
112 }
113
124 const MaterialLawParams &params,
125 Scalar lowerSat = 0.0,
126 Scalar upperSat = 1.0,
127 std::string curveTitle = "")
128 {
129 std::vector<Scalar> sw;
130 std::vector<Scalar> pc;
131 Scalar satInterval = upperSat - lowerSat;
132 Scalar pcMin = 0.0;
133 Scalar pcMax = -1e100;
134
135 Scalar swTemp, pcTemp = 0.0;
136 for (int i = 0; i <= numIntervals_; i++)
137 {
138 swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_);
139 pcTemp = MaterialLaw::pcnw(params, swTemp);
140 using std::max;
141 using std::min;
142 if (checkValues_(swTemp, pcTemp))
143 {
144 sw.push_back(swTemp);
145 pc.push_back(pcTemp);
146 pcMin = min(pcMin, pcTemp);
147 pcMax = max(pcMax, pcTemp);
148 }
149 }
150
151 gnuplot.setXlabel("wetting phase saturation [-]");
152 gnuplot.setYlabel("capillary pressure [Pa]");
153 gnuplot.addDataSetToPlot(sw, pc, curveTitle + "_pcnw-Sw");
154 }
155
166 const MaterialLawParams &params,
167 Scalar lowerSat = 0.0,
168 Scalar upperSat = 1.0,
169 std::string curveTitle = "")
170 {
171 std::vector<Scalar> st;
172 std::vector<Scalar> pc;
173 Scalar satInterval = upperSat - lowerSat;
174 Scalar pcMin = 0.0;
175 Scalar pcMax = -1e100;
176
177 Scalar stTemp, pcTemp = 0.0;
178 for (int i = 0; i <= numIntervals_; i++)
179 {
180 stTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_);
181 pcTemp = MaterialLaw::pcgn(params, stTemp);
182 using std::max;
183 using std::min;
184 if (checkValues_(stTemp, pcTemp))
185 {
186 st.push_back(stTemp);
187 pc.push_back(pcTemp);
188 pcMin = min(pcMin, pcTemp);
189 pcMax = max(pcMax, pcTemp);
190 }
191 }
192
193 gnuplot.setXlabel("wetting phase saturation [-]");
194 gnuplot.setYlabel("capillary pressure [Pa]");
195 gnuplot.addDataSetToPlot(st, pc, curveTitle + "_pcgn-St");
196 }
197
198
209 const MaterialLawParams &params,
210 Scalar lowerSat = 0.0,
211 Scalar upperSat = 1.0,
212 std::string curveTitle = "")
213 {
214 std::vector<Scalar> sw(numIntervals_ + 1);
215 std::vector<Scalar> krw(numIntervals_ + 1);
216 std::vector<Scalar> krn(numIntervals_ + 1);
217 std::vector<Scalar> krg(numIntervals_ + 1);
218 Scalar satInterval = upperSat - lowerSat;
219 Scalar krMin = 1e100;
220 Scalar krMax = -1e100;
221
222 Scalar swTemp, krwTemp, krnTemp, krgTemp = 0.0;
223 for (int i = 0; i <= numIntervals_; i++)
224 {
225 swTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_);
226 krwTemp = MaterialLaw::krw(params, swTemp, 0.0);
227 krnTemp = MaterialLaw::krn(params, swTemp, 1.0 - swTemp);
228 krgTemp = MaterialLaw::krg(params, swTemp, 0.0);
229 using std::max;
230 using std::min;
231 if (checkValues_(swTemp, krwTemp)
232 && checkValues_(swTemp, krnTemp)
233 && checkValues_(swTemp, krgTemp))
234 {
235 sw.push_back(swTemp);
236 krw.push_back(krwTemp);
237 krn.push_back(krnTemp);
238 krg.push_back(krgTemp);
239 krMin = min({krMin, krwTemp, krnTemp, krgTemp});
240 krMax = max({krMax, krwTemp, krnTemp, krgTemp});
241 }
242 }
243
244 gnuplot.setXlabel("wetting phase saturation [-]");
245 gnuplot.setYlabel("relative permeability [-]");
246 gnuplot.addDataSetToPlot(sw, krw, curveTitle + "_krw");
247 gnuplot.addDataSetToPlot(sw, krn, curveTitle + "_krn");
248 gnuplot.addDataSetToPlot(sw, krg, curveTitle + "_krg");
249 }
250
261 const MaterialLawParams &params,
262 Scalar lowerSat = 0.0,
263 Scalar upperSat = 1.0,
264 std::string curveTitle = "")
265 {
266 std::vector<Scalar> sn(numIntervals_ + 1);
267 std::vector<Scalar> alpha(numIntervals_ + 1);
268 Scalar satInterval = upperSat - lowerSat;
269 Scalar alphaMin = -2;
270 Scalar alphaMax = 2;
271
272 Scalar snTemp, alphaTemp = 0.0;
273 for (int i = 0; i <= numIntervals_; i++)
274 {
275 snTemp = lowerSat + satInterval * Scalar(i) / Scalar(numIntervals_);
276 alphaTemp = MaterialLaw::pcAlpha(params, snTemp);
277 using std::max;
278 using std::min;
279 if (checkValues_(snTemp, alphaTemp))
280 {
281 sn.push_back(snTemp);
282 alpha.push_back(alphaTemp);
283 alphaMin = min(alphaMin, alphaTemp);
284 alphaMax = max(alphaMax, alphaTemp);
285 }
286 }
287
288 gnuplot.setXlabel("nonwetting phase saturation [-]");
289 gnuplot.setYlabel("transition function [-]");
290 gnuplot.addDataSetToPlot(sn, alpha, curveTitle + "_alpha");
291 }
292
293private:
300 bool checkValues_(Scalar value1, Scalar value2)
301 {
302 using std::isnan;
303 using std::isinf;
304 return !isnan(value1) && !isinf(value1)
305 && !isnan(value2) && !isinf(value2);
306 }
307
308 int numIntervals_;
309};
310} // end of namespace Dumux
311
312#endif
Definition: adapt.hh:29
Interface for passing data sets to a file and plotting them, if gnuplot is installed.
Definition: gnuplotinterface.hh:57
void setXlabel(const std::string &label)
Sets the label for the x-axis.
Definition: gnuplotinterface.hh:279
void addDataSetToPlot(const std::vector< Scalar > &x, const std::vector< Scalar > &y, const std::string &fileName, const std::string &options="with lines")
Adds a data set and writes a data file.
Definition: gnuplotinterface.hh:242
void setYlabel(const std::string &label)
Sets the label for the y-axis.
Definition: gnuplotinterface.hh:289
Interface for plotting the three-phase fluid-matrix-interaction laws TODO: add theta head pressure pl...
Definition: plotmateriallaw3p.hh:43
void addPcAlpha(GnuplotInterface< Scalar > &gnuplot, const MaterialLawParams &params, Scalar lowerSat=0.0, Scalar upperSat=1.0, std::string curveTitle="")
Plot the transition (2P/3P) function.
Definition: plotmateriallaw3p.hh:260
void addpc(GnuplotInterface< Scalar > &gnuplot, const MaterialLawParams &params, Scalar lowerSat=0.0, Scalar upperSat=1.0, std::string curveTitle="")
Plot the capillary pressure-saturation curve for all phases.
Definition: plotmateriallaw3p.hh:61
void addpcgw(GnuplotInterface< Scalar > &gnuplot, const MaterialLawParams &params, Scalar lowerSat=0.0, Scalar upperSat=1.0, std::string curveTitle="")
Plot the capillary pressure-saturation curve for the water-gas interphase.
Definition: plotmateriallaw3p.hh:81
void addpcgn(GnuplotInterface< Scalar > &gnuplot, const MaterialLawParams &params, Scalar lowerSat=0.0, Scalar upperSat=1.0, std::string curveTitle="")
Plot the capillary pressure-saturation curve for the gas-NAPL interface.
Definition: plotmateriallaw3p.hh:165
PlotMaterialLaw(bool interaction=true)
Constructor.
Definition: plotmateriallaw3p.hh:48
void addkr(GnuplotInterface< Scalar > &gnuplot, const MaterialLawParams &params, Scalar lowerSat=0.0, Scalar upperSat=1.0, std::string curveTitle="")
Plot the relative permeabilities.
Definition: plotmateriallaw3p.hh:208
void addpcnw(GnuplotInterface< Scalar > &gnuplot, const MaterialLawParams &params, Scalar lowerSat=0.0, Scalar upperSat=1.0, std::string curveTitle="")
Plot the capillary pressure-saturation curve for the water-NAPL interface.
Definition: plotmateriallaw3p.hh:123