ImplicitCapillarity_impl.hpp
Go to the documentation of this file.
1//===========================================================================
2//
3// File: ImplicitCapillarity_impl.hpp
4//
5// Created: Thu May 6 15:36:07 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 OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
37#define OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
38
39#include <opm/common/ErrorMacros.hpp>
40#include <opm/core/utility/Average.hpp>
41#include <opm/core/utility/Units.hpp>
42#include <opm/core/utility/RootFinders.hpp>
43#include <opm/core/utility/StopWatch.hpp>
44#include <dune/grid/common/Volumes.hpp>
45#include <opm/porsol/common/ReservoirPropertyFixedMobility.hpp>
46#include <opm/porsol/euler/MatchSaturatedVolumeFunctor.hpp>
47
48#include <cassert>
49#include <cmath>
50#include <algorithm>
51#include <iostream>
52
53namespace Opm
54{
55
56
57 template <class GI, class RP, class BC, template <class, class> class IP>
58 inline ImplicitCapillarity<GI, RP, BC, IP>::ImplicitCapillarity()
59 : method_viscous_(true),
60 method_gravity_(true),
61 check_sat_(true),
62 clamp_sat_(false),
63 residual_tolerance_(1e-8),
64 linsolver_verbosity_(1),
65 linsolver_type_(1),
66 update_relaxation_(1.0)
67 {
68 }
69
70 template <class GI, class RP, class BC, template <class, class> class IP>
71 inline ImplicitCapillarity<GI, RP, BC, IP>::ImplicitCapillarity(const GI& g, const RP& r, const BC& b)
72 : residual_(g, r, b),
73 method_viscous_(true),
74 method_gravity_(true),
75 check_sat_(true),
76 clamp_sat_(false),
77 residual_tolerance_(1e-8),
78 linsolver_verbosity_(1),
79 linsolver_type_(1),
80 update_relaxation_(1.0)
81 {
82 Dune::FieldVector<double, GI::Dimension> grav(0.0);
83 psolver_.init(g, r, grav, b);
84 }
85
86
87
88 template <class GI, class RP, class BC, template <class, class> class IP>
89 inline void ImplicitCapillarity<GI, RP, BC, IP>::init(const Opm::parameter::ParameterGroup& param)
90 {
91 method_viscous_ = param.getDefault("method_viscous", method_viscous_);
92 method_gravity_ = param.getDefault("method_gravity", method_gravity_);
93 check_sat_ = param.getDefault("check_sat", check_sat_);
94 clamp_sat_ = param.getDefault("clamp_sat", clamp_sat_);
95 residual_tolerance_ = param.getDefault("residual_tolerance", residual_tolerance_);
96 linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_);
97 linsolver_type_ = param.getDefault("linsolver_type", linsolver_type_);
98 update_relaxation_ = param.getDefault("update_relaxation", update_relaxation_);
99 }
100
101 template <class GI, class RP, class BC, template <class, class> class IP>
102 inline void ImplicitCapillarity<GI, RP, BC, IP>::init(const Opm::parameter::ParameterGroup& param,
103 const GI& g, const RP& r, const BC& b)
104 {
105 init(param);
106 initObj(g, r, b);
107 }
108
109
110 template <class GI, class RP, class BC, template <class, class> class IP>
111 inline void ImplicitCapillarity<GI, RP, BC, IP>::initObj(const GI& g, const RP& r, const BC& b)
112 {
113 residual_.initObj(g, r, b);
114 Dune::FieldVector<double, GI::Dimension> grav(0.0);
115 psolver_.init(g, r, grav, b);
116
117 }
118
119
120
121 namespace ImplicitCapillarityDetails {
122 void thresholdMobility(double& m, double threshold);
123
124 // The matrix variant expects diagonal mobilities.
125 template <class SomeMatrixType>
126 void thresholdMobility(SomeMatrixType& m, double threshold)
127 {
128 for (int i = 0; i < std::min(m.numRows(), m.numCols()); ++i) {
129 m(i, i) = std::max(m(i, i), threshold);
130 }
131 }
132 } // anon namespace
133
134
135
136 template <class GI, class RP, class BC, template <class, class> class IP>
137 template <class PressureSolution>
138 void ImplicitCapillarity<GI, RP, BC, IP>::transportSolve(std::vector<double>& saturation,
139 const double /*time*/,
140 const typename GI::Vector& gravity,
141 const PressureSolution& pressure_sol,
142 const Opm::SparseVector<double>& injection_rates) const
143 {
144 // Start a timer.
145 Opm::time::StopWatch clock;
146 clock.start();
147
148 // Compute capillary mobilities.
149 typedef typename RP::Mobility Mob;
150 int num_cells = saturation.size();
151 std::vector<Mob> cap_mob(num_cells);
152 for (int c = 0; c < num_cells; ++c) {
153 Mob& m = cap_mob[c];
154 residual_.reservoirProperties().phaseMobility(0, c, saturation[c], m.mob);
155 Mob mob2;
156 residual_.reservoirProperties().phaseMobility(1, c, saturation[c], mob2.mob);
157 Mob mob_tot;
158 mob_tot.setToSum(m, mob2);
159 Mob mob_totinv;
160 mob_totinv.setToInverse(mob_tot);
161 m *= mob_totinv;
162 m *= mob2;
163 ImplicitCapillarityDetails::thresholdMobility(m.mob, 1e-10); // @@TODO: User-set limit.
164 // std::cout << m.mob(0,0) << '\n';
165 }
166 ReservoirPropertyFixedMobility<Mob> capillary_mobilities(cap_mob);
167
168 // Set up boundary conditions.
169 BC cap_press_bcs(residual_.boundaryConditions());
170 for (int i = 0; i < cap_press_bcs.size(); ++i) {
171 if (cap_press_bcs.flowCond(i).isPeriodic()) {
172 cap_press_bcs.flowCond(i) = FlowBC(FlowBC::Periodic, 0.0);
173 }
174 }
175
176 // Compute injection rates from residual.
177 std::vector<double> injection_rates_residual(num_cells);
178 residual_.computeResidual(saturation, gravity, pressure_sol, injection_rates,
179 method_viscous_, method_gravity_, false,
180 injection_rates_residual);
181 for (int i = 0; i < num_cells; ++i) {
182 injection_rates_residual[i] = -injection_rates_residual[i];
183 }
184
185 // Compute capillary pressure.
186 // Note that the saturation is just a dummy for this call, since the mobilities are fixed.
187 psolver_.solve(capillary_mobilities, saturation, cap_press_bcs, injection_rates_residual,
188 residual_tolerance_, linsolver_verbosity_, linsolver_type_);
189
190 // Solve for constant to change capillary pressure solution by.
191 std::vector<double> cap_press(num_cells);
192 const PressureSolution& pcapsol = psolver_.getSolution();
193 for (CIt c = residual_.grid().cellbegin(); c != residual_.grid().cellend(); ++c) {
194 cap_press[c->index()] = pcapsol.pressure(c);
195 }
196 MatchSaturatedVolumeFunctor<GI, RP> functor(residual_.grid(),
197 residual_.reservoirProperties(),
198 saturation,
199 cap_press);
200 double min_cap_press = *std::min_element(cap_press.begin(), cap_press.end());
201 double max_cap_press = *std::max_element(cap_press.begin(), cap_press.end());
202 double cap_press_range = max_cap_press - min_cap_press;
203 double mod_low = 1e100;
204 double mod_high = -1e100;
205 Opm::bracketZero(functor, 0.0, cap_press_range, mod_low, mod_high);
206 const int max_iter = 40;
207 const double nonlinear_tolerance = 1e-12;
208 int iterations_used = -1;
209 typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
210 double mod_correct = RootFinder::solve(functor, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
211 std::cout << "Moved capillary pressure solution by " << mod_correct << " after "
212 << iterations_used << " iterations." << std::endl;
213 // saturation = functor.lastSaturations();
214 const std::vector<double>& sat_new = functor.lastSaturations();
215 for (int i = 0; i < num_cells; ++i) {
216 saturation[i] = (1.0 - update_relaxation_)*saturation[i] + update_relaxation_*sat_new[i];
217 }
218
219 // Optionally check and/or clamp results.
220 if (check_sat_ || clamp_sat_) {
221 checkAndPossiblyClampSat(saturation);
222 }
223
224 // Stop timer and optionally print seconds taken.
225 clock.stop();
226#ifdef VERBOSE
227 std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
228#endif // VERBOSE
229 }
230
231
232
233 template <class GI, class RP, class BC, template <class, class> class IP>
234 inline void ImplicitCapillarity<GI, RP, BC, IP>::checkAndPossiblyClampSat(std::vector<double>& s) const
235 {
236 int num_cells = s.size();
237 for (int cell = 0; cell < num_cells; ++cell) {
238 if (s[cell] > 1.0 || s[cell] < 0.0) {
239 if (clamp_sat_) {
240 s[cell] = std::max(std::min(s[cell], 1.0), 0.0);
241 } else if (s[cell] > 1.001 || s[cell] < -0.001) {
242 OPM_THROW(std::runtime_error, "Saturation out of range in ImplicitCapillarity: Cell " << cell << " sat " << s[cell]);
243 }
244 }
245 }
246 }
247
248
249
250} // end namespace Opm
251
252
253
254#endif // OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
Opm::BCBase< double >::Periodic
@ Periodic
Definition: BoundaryConditions.hpp:58
Opm::FlowBC
A class for representing a flow boundary condition.
Definition: BoundaryConditions.hpp:122
Opm::ImplicitCapillarity::CIt
GridInterface::CellIterator CIt
Definition: ImplicitCapillarity.hpp:102
Opm::ImplicitCapillarity::initObj
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary)
Definition: ImplicitCapillarity_impl.hpp:111
Opm::ImplicitCapillarity::ImplicitCapillarity
ImplicitCapillarity()
Definition: ImplicitCapillarity_impl.hpp:58
Opm::ImplicitCapillarity::init
void init(const Opm::parameter::ParameterGroup &param)
Definition: ImplicitCapillarity_impl.hpp:89
Opm::ImplicitCapillarity::transportSolve
void transportSolve(std::vector< double > &saturation, const double time, const typename GridInterface::Vector &gravity, const PressureSolution &pressure_sol, const Opm::SparseVector< double > &injection_rates) const
Solve transport equation, evolving.
Opm::ImplicitCapillarity::psolver_
PressureSolver psolver_
Definition: ImplicitCapillarity.hpp:106
Opm::ImplicitCapillarity::checkAndPossiblyClampSat
void checkAndPossiblyClampSat(std::vector< double > &s) const
Definition: ImplicitCapillarity_impl.hpp:234
Opm::IncompFlowSolverHybrid::init
void init(const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
All-in-one initialization routine. Enumerates all grid connections, allocates sufficient space,...
Definition: IncompFlowSolverHybrid.hpp:477
Opm::ReservoirPropertyFixedMobility
Definition: ReservoirPropertyFixedMobility.hpp:46
Opm::ImplicitCapillarityDetails::thresholdMobility
void thresholdMobility(double &m, double threshold)
Opm
Definition: BlackoilFluid.hpp:32
Opm::MatchSaturatedVolumeFunctor
Definition: MatchSaturatedVolumeFunctor.hpp:69
Opm::MatchSaturatedVolumeFunctor::lastSaturations
const std::vector< double > & lastSaturations() const
Definition: MatchSaturatedVolumeFunctor.hpp:97