SteadyStateUpscalerManagerImplicit.hpp
Go to the documentation of this file.
1//===========================================================================
2//
3// File: SteadyStateUpscalerManager.hpp
4//
5// Created: Thu Apr 29 11:12:56 2010
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// Jostein R Natvig <jostein.r.natvig@sintef.no>
9//
10// $Date$
11//
12// $Revision$
13//
14//===========================================================================
15
16/*
17 Copyright 2010 SINTEF ICT, Applied Mathematics.
18 Copyright 2010 Statoil ASA.
19
20 This file is part of The Open Reservoir Simulator Project (OpenRS).
21
22 OpenRS is free software: you can redistribute it and/or modify
23 it under the terms of the GNU General Public License as published by
24 the Free Software Foundation, either version 3 of the License, or
25 (at your option) any later version.
26
27 OpenRS is distributed in the hope that it will be useful,
28 but WITHOUT ANY WARRANTY; without even the implied warranty of
29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 GNU General Public License for more details.
31
32 You should have received a copy of the GNU General Public License
33 along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34*/
35
36#ifndef OPM_STEADYSTATEUPSCALERMANAGERIMPLICIT_HEADER
37#define OPM_STEADYSTATEUPSCALERMANAGERIMPLICIT_HEADER
38
39
40#include <opm/upscaling/SteadyStateUpscalerImplicit.hpp>
41#include <opm/upscaling/UpscalingTraits.hpp>
42#include <opm/input/eclipse/Units/Units.hpp>
43#include <opm/grid/utility/SparseTable.hpp>
44#include <cmath>
45#include <fstream>
46#include <iostream>
47
48
49namespace Opm
50{
51
52
54 template <class Istream>
55 void readControl(Istream& is, std::vector<double>& saturations, Opm::SparseTable<double>& all_pdrops)
56 {
57 int num_sats;
58 is >> num_sats;
59 std::vector<double> sat(num_sats);
60 Opm::SparseTable<double> ap;
61 std::vector<double> tmp_pd;
62 for (int i = 0; i < num_sats; ++i) {
63 is >> sat[i];
64 int num_pdrops;
65 is >> num_pdrops;
66 tmp_pd.resize(num_pdrops);
67 for (int j = 0; j < num_pdrops; ++j) {
68 is >> tmp_pd[j];
69 }
70 ap.appendRow(tmp_pd.begin(), tmp_pd.end());
71 }
72 all_pdrops.swap(ap);
73 saturations.swap(sat);
74 }
75
76
77
78
80 template <class Ostream>
81 void writeControl(Ostream& os, const std::vector<double>& saturations, const Opm::SparseTable<double>& all_pdrops)
82 {
83 int num_sats = saturations.size();
84 os << num_sats << '\n';
85 for (int i = 0; i < num_sats; ++i) {
86 os << saturations[i] << '\n';
87 int num_pdrops = all_pdrops[i].size();
88 os << num_pdrops;
89 for (int j = 0; j < num_pdrops; ++j) {
90 os << ' ' << all_pdrops[i][j];
91 }
92 os << '\n';
93 }
94 }
95
96
97
98
99 template<class Ostream, class Tensor>
100 void writeRelPerm(Ostream& os, const Tensor& K, double sat, double pdrop,bool success)
101 {
102 /* const int num_rows = K.numRows(); */
103 /* const int num_cols = K.numCols(); */
104
105 /* We write tensor output in Voigt notation,
106 (but we also output the remainding three terms)
107 http://en.wikipedia.org/wiki/Voigt_notation
108 */
109
110 /*
111 TODO:
112 If fixed boundary condition, only output diagonal elements
113 If linear or periodic bc's, output all 9 elements (even though 6 is strictly necessary for periodic bc's)
114 Use the Tensor class to order output into Voigt notation, so that it works for more than 3x3 matrices
115 */
116
117 os << pdrop << '\t';
118 os << sat << '\t';
119 os << K(0,0) << '\t'; /* xx */
120 os << K(1,1) << '\t'; /* yy */
121 os << K(2,2) << '\t'; /* zz */
122 os << K(1,2) << '\t'; /* yz */
123 os << K(0,2) << '\t'; /* xz */
124 os << K(0,1) << '\t'; /* xy */
125 os << K(2,1) << '\t'; /* zy */
126 os << K(2,0) << '\t'; /* zx */
127 os << K(1,2) << '\t'; /* yz */
128 os << success;
129
130 os << std::endl;
131
132 }
133
134
135
136 template <class Upscaler>
137 class SteadyStateUpscalerManagerImplicit
138 {
139 public:
140 void upscale(const Opm::ParameterGroup& param)
141 {
142 // Control structure.
143 std::vector<double> saturations;
144 Opm::SparseTable<double> all_pdrops;
145 bool from_file = param.has("sat_pdrop_filename");
146 if (from_file) {
147 std::string filename = param.get<std::string>("sat_pdrop_filename");
148 std::ifstream file(filename.c_str());
149 if (!file) {
150 OPM_THROW(std::runtime_error, "Could not open file " + filename);
151 }
152 readControl(file, saturations, all_pdrops);
153 } else {
154 // Get a linear range of saturations.
155 int num_sats = param.getDefault("num_sats", 4);
156 double min_sat = param.getDefault("min_sat", 0.2);
157 double max_sat = param.getDefault("max_sat", 0.8);
158 saturations.resize(num_sats);
159 for (int i = 0; i < num_sats; ++i) {
160 double factor = num_sats == 1 ? 0 : double(i)/double(num_sats - 1);
161 saturations[i] = (1.0 - factor)*min_sat + factor*max_sat;
162 }
163 // Get a logarithmic range of pressure drops.
164 int num_pdrops = param.getDefault("num_pdrops", 5);
165 double log_min_pdrop = std::log(param.getDefault("min_pdrop", 1e2));
166 double log_max_pdrop = std::log(param.getDefault("max_pdrop", 1e6));
167 std::vector<double> pdrops;
168 pdrops.resize(num_pdrops);
169 for (int i = 0; i < num_pdrops; ++i) {
170 double factor = num_pdrops == 1 ? 0 : double(i)/double(num_pdrops - 1);
171 pdrops[i] = std::exp((1.0 - factor)*log_min_pdrop + factor*log_max_pdrop);
172 }
173 // Assign the same pressure drops to all saturations.
174 for (int i = 0; i < num_sats; ++i) {
175 all_pdrops.appendRow(pdrops.begin(), pdrops.end());
176 }
177 }
178 int flow_direction = param.getDefault("flow_direction", 0);
179 bool start_from_cl = param.getDefault("start_from_cl", true);
180
181 // Print the saturations and pressure drops.
182 // writeControl(std::cout, saturations, all_pdrops);
183
184 // Create output streams for upscaled relative permeabilities
185 std::string kr_filename = param.getDefault<std::string>("kr_filename", "upscaled_relperm");
186 std::string krwtmpname = kr_filename + "_water";
187 std::string krotmpname = kr_filename + "_oil";
188 std::string krw_filename = param.getDefault<std::string>("outputWater", krwtmpname);
189 std::string kro_filename = param.getDefault<std::string>("outputOil", krotmpname);
190 std::ofstream krw_out(krw_filename.c_str());
191 std::ofstream kro_out(kro_filename.c_str());
192 // std::stringstream krw_out;
193 // std::stringstream kro_out;
194 krw_out << "# Result from steady state upscaling" << std::endl;
195 krw_out << "# Pressuredrop Sw Krxx Kryy Krzz" << std::endl;
196 kro_out << "# Result from steady state upscaling" << std::endl;
197 kro_out << "# Pressuredrop Sw Krxx Kryy Krzz" << std::endl;
198 krw_out.precision(4); krw_out.setf(std::ios::scientific | std::ios::showpoint);
199 kro_out.precision(4); kro_out.setf(std::ios::scientific | std::ios::showpoint);
200
201 // Initialize upscaler.
202 Upscaler upscaler;
203 upscaler.init(param);
204
205 // First, compute an upscaled permeability.
206 typedef typename Upscaler::permtensor_t permtensor_t;
207 permtensor_t upscaled_K = upscaler.upscaleSinglePhase();
208 permtensor_t upscaled_K_copy = upscaled_K;
209 upscaled_K_copy *= (1.0/(Opm::prefix::milli*Opm::unit::darcy));
210 std::cout.precision(15);
211 std::cout << "Upscaled K in millidarcy:\n" << upscaled_K_copy << std::endl;
212 std::cout << "Upscaled porosity: " << upscaler.upscalePorosity() << std::endl;
213
214 // Then, compute some upscaled relative permeabilities.
215 int num_cells = upscaler.grid().size(0);
216 int num_sats = saturations.size();
217 //std::vector<bool> success_state;
218 for (int i = 0; i < num_sats; ++i) {
219 // Starting every computation with a trio of uniform profiles.
220 std::vector<double> init_sat;
221 if (start_from_cl) {
222 try {
223 upscaler.setToCapillaryLimit(saturations[i], init_sat);
224 } catch (...) {
225 init_sat.resize(num_cells, saturations[i]);
226 std::cout << "Failed to initialize with capillary limit for s = " << saturations[i]
227 << ". Init with uniform distribution." << std::endl;
228 }
229 } else {
230 init_sat.resize(num_cells, saturations[i]);
231 }
232
233 auto pdrops = all_pdrops[i];
234 int num_pdrops = pdrops.size();
235 for (int j = 0; j < num_pdrops; ++j) {
236 upscaler.initSatLimits(init_sat);
237 double pdrop = pdrops[j];
238 bool success = false;
239 std::pair<permtensor_t, permtensor_t> lambda
240 = upscaler.upscaleSteadyState(flow_direction, init_sat, saturations[i], pdrop, upscaled_K, success);
241 double usat = upscaler.lastSaturationUpscaled();
242 if(! success){
243 std::cout << "Upscaling failed for " << usat << std::endl;
244 }else{
245 init_sat = upscaler.lastSaturationState();
246 }
247 std::cout << "\n\nTensor of upscaled relperms for initial saturation " << saturations[i]
248 << ", real steady-state saturation " << usat
249 << " and pressure drop " << pdrop
250 << ":\n\n[water]\n" << lambda.first
251 << "\n[oil]\n" << lambda.second << std::endl;
252 // Changing initial saturations for next pressure drop to equal the steady state of the last
253 writeRelPerm(krw_out, lambda.first , usat, pdrop, success);
254 writeRelPerm(kro_out, lambda.second, usat, pdrop, success);
255
256 }
257 }
258 // std::cout << krw_out << std::endl;
259 // std::cout << kro_out << std::endl;
260 // if (!krw_filename.empty()) {
261 // // write water results to file
262 // std::ofstream krw_outfile;
263 // krw_outfile.open(krw_filename.c_str(), std::ios::out | std::ios::trunc);
264 // krw_outfile << krw_out.str();
265 // krw_outfile.close();
266 // }
267 // if (!kro_filename.empty()) {
268 // // write water results to file
269 // std::ofstream kro_outfile;
270 // kro_outfile.open(kro_filename.c_str());
271 // kro_outfile << kro_out.str();
272 // kro_outfile.close();
273 // }
274 }
275 };
276
277
278
279
280
281
282} // namespace Opm
283
284#endif // OPM_STEADYSTATEUPSCALERMANAGER_HEADER
Opm::SteadyStateUpscalerManagerImplicit
Definition: SteadyStateUpscalerManagerImplicit.hpp:138
Opm::SteadyStateUpscalerManagerImplicit::upscale
void upscale(const Opm::ParameterGroup &param)
Definition: SteadyStateUpscalerManagerImplicit.hpp:140
Opm::Elasticity::j
int j
Definition: elasticity_upscale_impl.hpp:658
Opm
Definition: ImplicitAssembly.hpp:43
Opm::writeRelPerm
void writeRelPerm(Ostream &os, const Tensor &K, double sat, double pdrop)
Definition: SteadyStateUpscalerManager.hpp:100
Opm::readControl
void readControl(Istream &is, std::vector< double > &saturations, Opm::SparseTable< double > &all_pdrops)
Reads saturation and pressure drop data from an input stream.
Definition: SteadyStateUpscalerManager.hpp:55
Opm::writeControl
void writeControl(Ostream &os, const std::vector< double > &saturations, const Opm::SparseTable< double > &all_pdrops)
Writes saturation and pressure drop data to an output stream.
Definition: SteadyStateUpscalerManager.hpp:81