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/core/utility/Units.hpp>
43 #include <opm/core/utility/SparseTable.hpp>
44 #include <cmath>
45 #include <fstream>
46 #include <iostream>
47 
48 
49 namespace 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::parameter::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  const Opm::SparseTable<double>::row_type 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:137
SteadyStateUpscalerImplicit.hpp
Opm
Definition: applier.hpp:18
UpscalingTraits.hpp
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
Opm::SteadyStateUpscalerManagerImplicit::upscale
void upscale(const Opm::parameter::ParameterGroup &param)
Definition: SteadyStateUpscalerManagerImplicit.hpp:140
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::writeRelPerm
void writeRelPerm(Ostream &os, const Tensor &K, double sat, double pdrop)
Definition: SteadyStateUpscalerManager.hpp:100
Opm::Elasticity::j
int j
Definition: elasticity_upscale_impl.hpp:658