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/input/eclipse/Units/Units.hpp>
42#include <opm/common/utility/numeric/RootFinders.hpp>
43#include <opm/grid/utility/StopWatch.hpp>
44#include <opm/grid/common/Volumes.hpp>
45#include <opm/porsol/common/ReservoirPropertyFixedMobility.hpp>
46#include <opm/porsol/euler/MatchSaturatedVolumeFunctor.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);
83 psolver_.
init(g, r, grav, b);
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 {
122 void thresholdMobility(
double& m,
double threshold);
125 template <
class SomeMatrixType>
126 void thresholdMobility(SomeMatrixType& m,
double threshold)
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);
163 ImplicitCapillarityDetails::thresholdMobility(m.mob, 1e-10);
166 ReservoirPropertyFixedMobility<Mob> capillary_mobilities(cap_mob);
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);
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);
196 MatchSaturatedVolumeFunctor<GI, RP> functor(residual_.grid(),
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;
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];
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>
234 inline void ImplicitCapillarity<GI, RP, BC, IP>::checkAndPossiblyClampSat(std::vector<double>& s)
const
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,
243 "Saturation out of range in ImplicitCapillarity: "
244 "Cell " + std::to_string(cell) +
245 " sat " + std::to_string(s[cell]));
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::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.
void init(const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
All-in-one initialization routine.
Definition IncompFlowSolverHybrid.hpp:476
Inverting small matrices.
Definition ImplicitAssembly.hpp:43