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