All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
UpscalerBase_impl.hpp
1 //===========================================================================
2 //
3 // File: UpscalerBase_impl.hpp
4 //
5 // Created: Thu Apr 29 10:22:06 2010
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 //
9 // $Date$
10 //
11 // $Revision$
12 //
13 //===========================================================================
14 
15 /*
16  Copyright 2010 SINTEF ICT, Applied Mathematics.
17  Copyright 2010 Statoil ASA.
18 
19  This file is part of The Open Reservoir Simulator Project (OpenRS).
20 
21  OpenRS is free software: you can redistribute it and/or modify
22  it under the terms of the GNU General Public License as published by
23  the Free Software Foundation, either version 3 of the License, or
24  (at your option) any later version.
25 
26  OpenRS is distributed in the hope that it will be useful,
27  but WITHOUT ANY WARRANTY; without even the implied warranty of
28  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29  GNU General Public License for more details.
30 
31  You should have received a copy of the GNU General Public License
32  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
33 */
34 
35 #ifndef OPM_UPSCALERBASE_IMPL_HEADER
36 #define OPM_UPSCALERBASE_IMPL_HEADER
37 
38 #include <opm/porsol/common/setupGridAndProps.hpp>
39 #include <opm/porsol/common/setupBoundaryConditions.hpp>
40 #include <opm/porsol/common/ReservoirPropertyTracerFluid.hpp>
41 
42 #include <iostream>
43 
44 namespace Opm
45 {
46 
47  template <class Traits>
49  : bctype_(Fixed),
50  twodim_hack_(false),
51  residual_tolerance_(1e-8),
52  linsolver_maxit_(0),
53  linsolver_prolongate_factor_(1.0),
54  linsolver_verbosity_(0),
55  linsolver_type_(3),
56  linsolver_smooth_steps_(1)
57  {
58  }
59 
60 
61 
62 
63  template <class Traits>
64  inline void UpscalerBase<Traits>::init(const Opm::ParameterGroup& param)
65  {
66  initImpl(param);
67  initFinal(param);
68  }
69 
70 
71 
72 
73  template <class Traits>
74  inline void UpscalerBase<Traits>::initImpl(const Opm::ParameterGroup& param)
75  {
76  // Request the boundary condition type parameter early since,
77  // depending on the actual type, we may have to manufacture
78  // and insert other parameters into the ParameterGroup prior
79  // to constructing the grid and associated properties.
80  //
81  int bct = param.get<int>("boundary_condition_type");
82  bctype_ = static_cast<BoundaryConditionType>(bct);
83 
84  // Import control parameters pertaining to reduced physical
85  // dimensionality ("2d_hack = true" precludes periodic
86  // boundary conditions in the Z direction), and linear solves.
87  //
88  twodim_hack_ = param.getDefault("2d_hack", twodim_hack_);
89  residual_tolerance_ = param.getDefault("residual_tolerance", residual_tolerance_);
90  linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_);
91  linsolver_type_ = param.getDefault("linsolver_type", linsolver_type_);
92  linsolver_maxit_ = param.getDefault("linsolver_max_iterations", linsolver_maxit_);
93  linsolver_prolongate_factor_ = param.getDefault("linsolver_prolongate_factor", linsolver_prolongate_factor_);
94  linsolver_smooth_steps_ = param.getDefault("linsolver_smooth_steps", linsolver_smooth_steps_);
95 
96  // Ensure sufficient grid support for requested boundary
97  // condition type.
98  //
99  Opm::ParameterGroup temp_param = param;
100  if (bctype_ == Linear || bctype_ == Periodic) {
101  if (!temp_param.has("use_unique_boundary_ids")) {
102  temp_param.insertParameter("use_unique_boundary_ids", "true");
103  }
104  if (!temp_param.has("clip_z")) {
105  temp_param.insertParameter("clip_z", "true");
106  }
107  }
108  if (bctype_ == Periodic) {
109  if (!temp_param.has("periodic_extension")) {
110  temp_param.insertParameter("periodic_extension", "true");
111  }
112  }
113 
114  setupGridAndProps(temp_param, grid_, res_prop_);
115  ginterf_.init(grid_);
116  }
117 
118 
119 
120 
121  template <class Traits>
122  inline void UpscalerBase<Traits>::initFinal(const Opm::ParameterGroup& param)
123  {
124  // Report any unused parameters.
125  std::cout << "==================== Unused parameters: ====================\n";
126  param.displayUsage();
127  std::cout << "================================================================\n";
128  }
129 
130 
131 
132 
133  template <class Traits>
134  inline void UpscalerBase<Traits>::init(const Opm::Deck& deck,
135  BoundaryConditionType bctype,
136  double perm_threshold,
137  double residual_tolerance,
138  int linsolver_verbosity,
139  int linsolver_type,
140  bool twodim_hack,
141  int linsolver_maxit,
142  double linsolver_prolongate_factor,
143  int linsolver_smooth_steps,
144  const double gravity)
145  {
146  bctype_ = bctype;
147  residual_tolerance_ = residual_tolerance;
148  linsolver_verbosity_ = linsolver_verbosity;
149  linsolver_type_ = linsolver_type;
150  linsolver_maxit_ = linsolver_maxit;
151  linsolver_prolongate_factor_ = linsolver_prolongate_factor;
152  linsolver_smooth_steps_ = linsolver_smooth_steps;
153  twodim_hack_ = twodim_hack;
154  gravity_ = gravity;
155 
156 
157  // Faking some parameters depending on bc type.
158  bool periodic_ext = (bctype_ == Periodic);
159  bool turn_normals = false;
160  bool clip_z = (bctype_ == Linear || bctype_ == Periodic);
161  bool unique_bids = (bctype_ == Linear || bctype_ == Periodic);
162  std::string rock_list("no_list");
164  periodic_ext, turn_normals, clip_z, unique_bids,
165  perm_threshold, rock_list,
166  useJ<ResProp>(), 1.0, 0.0,
167  grid_, res_prop_);
168  ginterf_.init(grid_);
169  }
170 
171 
172 
173 
174  template <class Traits>
175  inline const typename UpscalerBase<Traits>::GridType&
177  {
178  return grid_;
179  }
180 
181 
182 
183 
184  template <class Traits>
185  inline void
187  {
188  if ((type == Periodic && bctype_ != Periodic)
189  || (type != Periodic && bctype_ == Periodic)) {
190  OPM_THROW(std::runtime_error, "Cannot switch to or from Periodic boundary condition, "
191  "periodic must be set in init() params.");
192  } else {
193  bctype_ = type;
194  if (type == Periodic || type == Linear) {
195  grid_.setUniqueBoundaryIds(true);
196  } else {
197  grid_.setUniqueBoundaryIds(false);
198  }
199  }
200  }
201 
202 
203 
204 
205  template <class Traits>
206  inline void
207  UpscalerBase<Traits>::setPermeability(const int cell_index, const permtensor_t& k)
208  {
209  res_prop_.permeabilityModifiable(cell_index) = k;
210  }
211 
212 
213 
214 
215  template <class Traits>
216  inline typename UpscalerBase<Traits>::permtensor_t
218  {
220  return upscaleEffectivePerm(fluid);
221  }
222 
223 
224 
225 
226  template <class Traits>
227  template <class FluidInterface>
228  inline typename UpscalerBase<Traits>::permtensor_t
229  UpscalerBase<Traits>::upscaleEffectivePerm(const FluidInterface& fluid)
230  {
231  int num_cells = ginterf_.numberOfCells();
232  // No source or sink.
233  std::vector<double> src(num_cells, 0.0);
234  // Just water.
235  std::vector<double> sat(num_cells, 1.0);
236  // Gravity.
237  Dune::FieldVector<double, 3> gravity(0.0);
238  gravity[2] = gravity_;
239 
240  permtensor_t upscaled_K(3, 3, (double*)0);
241  for (int pdd = 0; pdd < Dimension; ++pdd) {
242  setupUpscalingConditions(ginterf_, bctype_, pdd, 1.0, 1.0, twodim_hack_, bcond_);
243  if (pdd == 0) {
244  // Only on first iteration, since we do not change the
245  // structure of the system, the way the flow solver is
246  // implemented.
247  flow_solver_.init(ginterf_, res_prop_, gravity, bcond_);
248  }
249 
250  // Run pressure solver.
251  bool same_matrix = (bctype_ != Fixed) && (pdd != 0);
252  flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_,
253  linsolver_verbosity_,
254  linsolver_type_, same_matrix,
255  linsolver_maxit_, linsolver_prolongate_factor_,
256  linsolver_smooth_steps_);
257  double max_mod = flow_solver_.postProcessFluxes();
258  std::cout << "Max mod = " << max_mod << std::endl;
259 
260  // Compute upscaled K.
261  double Q[Dimension] = { 0 };
262  switch (bctype_) {
263  case Fixed:
264  Q[pdd] = computeAverageVelocity(flow_solver_.getSolution(), pdd, pdd);
265  break;
266  case Linear:
267  case Periodic:
268  for (int i = 0; i < Dimension; ++i) {
269  Q[i] = computeAverageVelocity(flow_solver_.getSolution(), i, pdd);
270  }
271  break;
272  default:
273  OPM_THROW(std::runtime_error, "Unknown boundary type: " << bctype_);
274  }
275  double delta = computeDelta(pdd);
276  for (int i = 0; i < Dimension; ++i) {
277  upscaled_K(i, pdd) = Q[i] * delta;
278  }
279  }
280  return upscaled_K;
281  }
282 
283 
284 
285 
286  template <class Traits>
287  template <class FlowSol>
288  double UpscalerBase<Traits>::computeAverageVelocity(const FlowSol& flow_solution,
289  const int flow_dir,
290  const int pdrop_dir) const
291  {
292  double side1_flux = 0.0;
293  double side2_flux = 0.0;
294  double side1_area = 0.0;
295  double side2_area = 0.0;
296 
297  int num_faces = 0;
298  int num_bdyfaces = 0;
299  int num_side1 = 0;
300  int num_side2 = 0;
301 
302  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
303  for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
304  ++num_faces;
305  if (f->boundary()) {
306  ++num_bdyfaces;
307  int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
308  if ((canon_bid - 1)/2 == flow_dir) {
309  double flux = flow_solution.outflux(f);
310  double area = f->area();
311  double norm_comp = f->normal()[flow_dir];
312  // std::cout << "bid " << f->boundaryId() << " area " << area << " n " << norm_comp << std::endl;
313  if (canon_bid - 1 == 2*flow_dir) {
314  ++num_side1;
315  if (flow_dir == pdrop_dir && flux > 0.0) {
316 #ifdef VERBOSE
317  std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()<<" (canonical: "<<canon_bid
318  << ") Magnitude: " << std::fabs(flux) << std::endl;
319 #endif
320  // OPM_THROW(std::runtime_error, "Detected outflow at entry face: " << face);
321  }
322  side1_flux += flux*norm_comp;
323  side1_area += area;
324  } else {
325  assert(canon_bid - 1 == 2*flow_dir + 1);
326  ++num_side2;
327  if (flow_dir == pdrop_dir && flux < 0.0) {
328 #ifdef VERBOSE
329  std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()
330  << " Magnitude: " << std::fabs(flux) << std::endl;
331 #endif
332  // OPM_THROW(std::runtime_error, "Detected inflow at exit face: " << face);
333  }
334  side2_flux += flux*norm_comp;
335  side2_area += area;
336  }
337  }
338  }
339  }
340  }
341 // std::cout << "Faces: " << num_faces << " Boundary faces: " << num_bdyfaces
342 // << " Side 1 faces: " << num_side1 << " Side 2 faces: " << num_side2 << std::endl;
343  // q is the average velocity.
344  return 0.5*(side1_flux/side1_area + side2_flux/side2_area);
345  }
346 
347 
348 
349 
350  template <class Traits>
351  inline double UpscalerBase<Traits>::computeDelta(const int flow_dir) const
352  {
353  double side1_pos = 0.0;
354  double side2_pos = 0.0;
355  double side1_area = 0.0;
356  double side2_area = 0.0;
357  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
358  for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
359  if (f->boundary()) {
360  int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
361  if ((canon_bid - 1)/2 == flow_dir) {
362  double area = f->area();
363  double pos_comp = f->centroid()[flow_dir];
364  if (canon_bid - 1 == 2*flow_dir) {
365  side1_pos += area*pos_comp;
366  side1_area += area;
367  } else {
368  side2_pos += area*pos_comp;
369  side2_area += area;
370  }
371  }
372  }
373  }
374  }
375  // delta is the average length.
376  return side2_pos/side2_area - side1_pos/side1_area;
377  }
378 
379 
380 
381 
382  template <class Traits>
384  {
385  double total_vol = 0.0;
386  double total_pore_vol = 0.0;
387  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
388  total_vol += c->volume();
389  total_pore_vol += c->volume()*res_prop_.porosity(c->index());
390  }
391  return total_pore_vol/total_vol;
392  }
393 
394 
395  template <class Traits>
397  {
398  double total_net_vol = 0.0;
399  double total_pore_vol = 0.0;
400  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
401  total_net_vol += c->volume()*res_prop_.ntg(c->index());
402  total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
403  }
404  if (total_net_vol>0.0) return total_pore_vol/total_net_vol;
405  else return 0.0;
406  }
407 
408  template <class Traits>
410  {
411  double total_vol = 0.0;
412  double total_net_vol = 0.0;
413  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
414  total_vol += c->volume();
415  total_net_vol += c->volume()*res_prop_.ntg(c->index());
416  }
417  return total_net_vol/total_vol;
418  }
419 
420  template <class Traits>
421  double UpscalerBase<Traits>::upscaleSWCR(const bool NTG) const
422  {
423  double total_swcr = 0.0;
424  double total_pore_vol = 0.0;
425  if (NTG) {
426  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
427  total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.swcr(c->index());
428  total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
429  }
430  }
431  else {
432  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
433  total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.swcr(c->index());
434  total_pore_vol += c->volume()*res_prop_.porosity(c->index());
435  }
436  }
437  return total_swcr/total_pore_vol;
438  }
439 
440  template <class Traits>
441  double UpscalerBase<Traits>::upscaleSOWCR(const bool NTG) const
442  {
443  double total_sowcr = 0.0;
444  double total_pore_vol = 0.0;
445  if (NTG) {
446  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
447  total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.sowcr(c->index());
448  total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
449  }
450  }
451  else {
452  for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
453  total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.sowcr(c->index());
454  total_pore_vol += c->volume()*res_prop_.porosity(c->index());
455  }
456  }
457  return total_sowcr/total_pore_vol;
458  }
459 
460 } // namespace Opm
461 
462 
463 
464 #endif // OPM_UPSCALERBASE_IMPL_HEADER
void setupGridAndPropsEclipse(const Opm::Deck &deck, bool periodic_extension, bool turn_normals, bool clip_z, bool unique_bids, double perm_threshold, const std::string &rock_list, bool use_jfunction_scaling, double sigma, double theta, Dune::CpGrid &grid, ResProp< 3 > &res_prop)
Definition: setupGridAndProps.hpp:151
void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition: setupBoundaryConditions.hpp:99
Definition: GridInterfaceEuler.hpp:376
A base class for upscaling.
Definition: UpscalerBase.hpp:55
void setPermeability(const int cell_index, const permtensor_t &k)
Set the permeability of a cell directly.
Definition: UpscalerBase_impl.hpp:207
void setBoundaryConditionType(BoundaryConditionType type)
Set boundary condition type.
Definition: UpscalerBase_impl.hpp:186
const GridType & grid() const
Access the grid.
Definition: UpscalerBase_impl.hpp:176
void setupGridAndProps(const Opm::ParameterGroup &param, Dune::CpGrid &grid, ResProp< 3 > &res_prop)
Definition: setupGridAndProps.hpp:70
double upscaleNetPorosity() const
Compute upscaled net porosity.
Definition: UpscalerBase_impl.hpp:396
Definition: ReservoirPropertyTracerFluid.hpp:40
double upscalePorosity() const
Compute upscaled porosity.
Definition: UpscalerBase_impl.hpp:383
UpscalerBase()
Default constructor.
Definition: UpscalerBase_impl.hpp:48
double upscaleSOWCR(const bool NTG) const
Compute upscaled SOWCR.
Definition: UpscalerBase_impl.hpp:441
ResProp::MutablePermTensor permtensor_t
A type for the upscaled permeability.
Definition: UpscalerBase.hpp:66
double upscaleNTG() const
Compute upscaled NTG.
Definition: UpscalerBase_impl.hpp:409
void init(const Opm::ParameterGroup &param)
Initializes the upscaler from parameters.
Definition: UpscalerBase_impl.hpp:64
permtensor_t upscaleSinglePhase()
Does a single-phase upscaling.
Definition: UpscalerBase_impl.hpp:217
double upscaleSWCR(const bool NTG) const
Compute upscaled SWCR.
Definition: UpscalerBase_impl.hpp:421