Go to the documentation of this file.
36#ifndef OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
37#define OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
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>
57 template < class GI, class RP, class BC, template < class, class> class IP>
59 : method_viscous_(true),
60 method_gravity_(true),
63 residual_tolerance_(1e-8),
64 linsolver_verbosity_(1),
66 update_relaxation_(1.0)
70 template < class GI, class RP, class BC, template < class, class> class IP>
73 method_viscous_(true),
74 method_gravity_(true),
77 residual_tolerance_(1e-8),
78 linsolver_verbosity_(1),
80 update_relaxation_(1.0)
82 Dune::FieldVector<double, GI::Dimension> grav(0.0);
88 template < class GI, class RP, class BC, template < class, class> class IP>
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_);
101 template < class GI, class RP, class BC, template < class, class> class IP>
103 const GI& g, const RP& r, const BC& b)
110 template < class GI, class RP, class BC, template < class, class> class IP>
113 residual_.initObj(g, r, b);
114 Dune::FieldVector<double, GI::Dimension> grav(0.0);
115 psolver_.init(g, r, grav, b);
121 namespace ImplicitCapillarityDetails {
125 template < class SomeMatrixType>
128 for ( int i = 0; i < std::min(m.numRows(), m.numCols()); ++i) {
129 m(i, i) = std::max(m(i, i), threshold);
136 template < class GI, class RP, class BC, template < class, class> class IP>
137 template < class PressureSolution>
140 const typename GI::Vector& gravity,
141 const PressureSolution& pressure_sol,
142 const Opm::SparseVector<double>& injection_rates) const
145 Opm::time::StopWatch clock;
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) {
154 residual_.reservoirProperties().phaseMobility(0, c, saturation[c], m.mob);
156 residual_.reservoirProperties().phaseMobility(1, c, saturation[c], mob2.mob);
158 mob_tot.setToSum(m, mob2);
160 mob_totinv.setToInverse(mob_tot);
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()) {
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];
187 psolver_.solve(capillary_mobilities, saturation, cap_press_bcs, injection_rates_residual,
188 residual_tolerance_, linsolver_verbosity_, linsolver_type_);
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);
197 residual_.reservoirProperties(),
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;
215 for ( int i = 0; i < num_cells; ++i) {
216 saturation[i] = (1.0 - update_relaxation_)*saturation[i] + update_relaxation_*sat_new[i];
220 if (check_sat_ || clamp_sat_) {
221 checkAndPossiblyClampSat(saturation);
227 std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
233 template < class GI, class RP, class BC, template < class, class> class IP>
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) {
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]);
@ Periodic Definition: BoundaryConditions.hpp:58
A class for representing a flow boundary condition. Definition: BoundaryConditions.hpp:122
GridInterface::CellIterator CIt Definition: ImplicitCapillarity.hpp:102
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary) Definition: ImplicitCapillarity_impl.hpp:111
ImplicitCapillarity() Definition: ImplicitCapillarity_impl.hpp:58
void init(const Opm::parameter::ParameterGroup ¶m) Definition: ImplicitCapillarity_impl.hpp:89
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.
PressureSolver psolver_ Definition: ImplicitCapillarity.hpp:106
void checkAndPossiblyClampSat(std::vector< double > &s) const Definition: ImplicitCapillarity_impl.hpp:234
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
Definition: ReservoirPropertyFixedMobility.hpp:46
void thresholdMobility(double &m, double threshold)
Definition: BlackoilFluid.hpp:32
Definition: MatchSaturatedVolumeFunctor.hpp:69
const std::vector< double > & lastSaturations() const Definition: MatchSaturatedVolumeFunctor.hpp:97
|