My Project
ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp
1 //===========================================================================
2 //
3 // File: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp
4 //
5 // Created: Mon Oct 26 10:12:15 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // B�rd Skaflestad <bard.skaflestad@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 2009, 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 OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
37 #define OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
38 
39 namespace Opm
40 {
41 
42  template <int dim>
43  template <class MatrixType>
45  int cell_index,
46  double saturation,
47  MatrixType& phase_mob) const
48  {
49  const int region = Super::rock_.size() > 0 ? Super::cell_to_rock_[cell_index] : -1;
50  phaseMobilityByRock(phase_index, region, saturation, phase_mob);
51  }
52 
53 
54  template <int dim>
55  double ReservoirPropertyCapillaryAnisotropicRelperm<dim>::fractionalFlow(int cell_index, double saturation) const
56  {
57  // This method is a hack.
58  // Assumes that the relperm is diagonal.
59  Mobility m;
60  double ff_first = 0.0;
61  for (int direction = 0; direction < dim; ++direction) {
62  phaseMobility(0, cell_index, saturation, m.mob);
63  double l1 = m.mob(direction, direction);
64  phaseMobility(1, cell_index, saturation, m.mob);
65  double l2 = m.mob(direction, direction);
66  ff_first += l1/(l1 + l2);
67  }
68  ff_first /= double(dim);
69  return ff_first;
70  }
71 
72 
73  template <int dim>
74  template <class MatrixType>
76  int rock_index,
77  double saturation,
78  MatrixType& phase_mob) const
79  {
80  assert ((0 <= phase_index) && (phase_index < Super::NumberOfPhases));
81  static_assert(Super::NumberOfPhases == 2, "");
82 
83  double visc = phase_index == 0 ? Super::viscosity1_ : Super::viscosity2_;
84  if (rock_index != -1) {
85  assert (rock_index < int(Super::rock_.size()));
86  Super::rock_[rock_index].kr(phase_index, saturation, phase_mob);
87  std::transform(phase_mob.data(), phase_mob.data() + dim*dim,
88  phase_mob.data(), [visc](double param) { return param / visc; });
89  } else {
90  // HACK: With no rocks, we use quadratic relperm functions.
91  double kr = phase_index == 0 ? saturation*saturation : (1.0 - saturation)*(1.0 - saturation);
92  double mob = kr/visc;
93  zero(phase_mob);
94  for (int j = 0; j < dim; ++j) {
95  phase_mob(j,j) = mob;
96  }
97  }
98  }
99 
100 
101  template <int dim>
102  void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::cflFracFlows(int rock, int direction,
103  double s, double& ff_first, double& ff_gravity) const
104  {
105  // Assumes that the relperm is diagonal.
106  Mobility m;
107  phaseMobilityByRock(0, rock, s, m.mob);
108  double l1 = m.mob(direction, direction);
109  phaseMobilityByRock(1, rock, s, m.mob);
110  double l2 = m.mob(direction, direction);
111  ff_first = l1/(l1 + l2);
112  ff_gravity = l1*l2/(l1 + l2);
113  }
114 
115 
116 
117 
118  template <int dim>
119  std::array<double, 3> ReservoirPropertyCapillaryAnisotropicRelperm<dim>::computeSingleRockCflFactors(int rock) const
120  {
121  const int N = 257;
122  double delta = 1.0/double(N - 1);
123  double last_ff1, last_ffg;
124  double max_der1 = -1e100;
125  double max_derg = -1e100;
126  double min_ffg = -1e100;
127  for (int direction = 0; direction < dim; ++direction) {
128  cflFracFlows(rock, direction, 0.0, last_ff1, last_ffg);
129  min_ffg = std::min(min_ffg, last_ffg);
130  for (int i = 1; i < N; ++i) {
131  double s = double(i)*delta;
132  double ff1, ffg;
133  cflFracFlows(rock, direction, s, ff1, ffg);
134  double est_deriv_ff1 = std::fabs(ff1 - last_ff1)/delta;
135  double est_deriv_ffg = std::fabs(ffg - last_ffg)/delta;
136  max_der1 = std::max(max_der1, est_deriv_ff1);
137  max_derg = std::max(max_derg, est_deriv_ffg);
138  min_ffg = std::min(min_ffg, last_ffg);
139  last_ff1 = ff1;
140  last_ffg = ffg;
141  }
142  }
143  std::array<double, 3> retval = {{ 1.0/max_der1, 1.0/max_derg, min_ffg }};
144  return retval;
145  }
146 
147 
148 
149 
150  template <int dim>
151  void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::computeCflFactors()
152  {
153  if (Super::rock_.empty()) {
154  std::array<double, 3> fac = computeSingleRockCflFactors(-1);
155  Super::cfl_factor_ = fac[0];
156  Super::cfl_factor_gravity_ = fac[1];
157  Super::cfl_factor_capillary_ = fac[2];
158  } else {
159  Super::cfl_factor_ = 1e100;
160  Super::cfl_factor_gravity_ = 1e100;
161  Super::cfl_factor_capillary_ = 1e100;
162  for (int r = 0; r < int(Super::rock_.size()); ++r) {
163  std::array<double, 3> fac = computeSingleRockCflFactors(r);
164  Super::cfl_factor_ = std::min(Super::cfl_factor_, fac[0]);
165  Super::cfl_factor_gravity_ = std::min(Super::cfl_factor_gravity_, fac[1]);
166  Super::cfl_factor_capillary_ = std::max(Super::cfl_factor_capillary_, fac[2]);
167  }
168  }
169  }
170 
171 
172 
173 
174 
175 
176 } // namespace Opm
177 
178 
179 #endif // OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
A property class for incompressible two-phase flow.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm.hpp:106
double fractionalFlow(int cell_index, double saturation) const
Some approximation to a scalar fractional flow (of the first phase).
Definition: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:55
void phaseMobility(int phase_index, int cell_index, double saturation, MatrixType &phase_mob) const
Anisotropic phase mobility.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:44
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:602
A wrapper for a tensor.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm.hpp:50