36#ifndef OPM_STEADYSTATEUPSCALER_IMPL_HEADER
37#define OPM_STEADYSTATEUPSCALER_IMPL_HEADER
43#include <opm/input/eclipse/Units/Units.hpp>
52 template <
class Traits>
57 print_inoutflows_(false),
58 simulation_steps_(10),
60 relperm_threshold_(1.0e-8),
61 maximum_mobility_contrast_(1.0e9),
62 sat_change_threshold_(0.0)
69 template <
class Traits>
72 Super::initImpl(param);
73 use_gravity_ = param.getDefault(
"use_gravity", use_gravity_);
74 output_vtk_ = param.getDefault(
"output_vtk", output_vtk_);
75 print_inoutflows_ = param.getDefault(
"print_inoutflows", print_inoutflows_);
76 simulation_steps_ = param.getDefault(
"simulation_steps", simulation_steps_);
77 stepsize_ = Opm::unit::convert::from(param.getDefault(
"stepsize", stepsize_),
79 relperm_threshold_ = param.getDefault(
"relperm_threshold", relperm_threshold_);
80 maximum_mobility_contrast_ = param.getDefault(
"maximum_mobility_contrast", maximum_mobility_contrast_);
81 sat_change_threshold_ = param.getDefault(
"sat_change_threshold", sat_change_threshold_);
83 transport_solver_.init(param);
85 double v1_default = this->res_prop_.viscosityFirstPhase();
86 double v2_default = this->res_prop_.viscositySecondPhase();
87 this->res_prop_.setViscosities(param.getDefault(
"viscosity1", v1_default), param.getDefault(
"viscosity2", v2_default));
88 double d1_default = this->res_prop_.densityFirstPhase();
89 double d2_default = this->res_prop_.densitySecondPhase();
90 this->res_prop_.setDensities(param.getDefault(
"density1", d1_default), param.getDefault(
"density2", d2_default));
97 double maxMobility(
double m1,
double m2)
99 return std::max(m1, m2);
102 template <
class SomeMatrixType>
103 double maxMobility(
double m1, SomeMatrixType& m2)
106 for (
int i = 0; i <
std::min(m2.numRows(), m2.numCols()); ++i) {
107 m = std::max(m, m2(i,i));
113 m = std::max(m, threshold);
116 template <
class SomeMatrixType>
119 for (
int i = 0; i <
std::min(m.numRows(), m.numCols()); ++i) {
120 m(i, i) = std::max(m(i, i), threshold);
128 template <
class Traits>
129 inline std::pair<typename SteadyStateUpscaler<Traits>::permtensor_t,
133 const std::vector<double>& initial_saturation,
134 const double boundary_saturation,
135 const double pressure_drop,
138 static int count = 0;
140 int num_cells = this->ginterf_.numberOfCells();
142 std::vector<double> src(num_cells, 0.0);
143 Opm::SparseVector<double> injection(num_cells);
145 Dune::FieldVector<double, 3> gravity(0.0);
147 gravity[2] = Opm::unit::gravity;
149 if (gravity.two_norm() > 0.0) {
150 OPM_MESSAGE(
"Warning: Gravity is experimental for flow solver.");
154 std::vector<double> saturation = initial_saturation;
158 pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
161 if (flow_direction == 0) {
162 this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
164 transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
167 this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
168 this->residual_tolerance_, this->linsolver_verbosity_,
169 this->linsolver_type_,
false,
170 this->linsolver_maxit_, this->linsolver_prolongate_factor_,
171 this->linsolver_smooth_steps_);
172 double max_mod = this->flow_solver_.postProcessFluxes();
173 std::cout <<
"Max mod = " << max_mod << std::endl;
176 std::vector<double> saturation_old = saturation;
177 for (
int iter = 0; iter < simulation_steps_; ++iter) {
179 transport_solver_.transportSolve(saturation, stepsize_, gravity, this->flow_solver_.getSolution(), injection);
182 this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
183 this->residual_tolerance_, this->linsolver_verbosity_,
184 this->linsolver_type_,
false,
185 this->linsolver_maxit_, this->linsolver_prolongate_factor_,
186 this->linsolver_smooth_steps_);
187 max_mod = this->flow_solver_.postProcessFluxes();
188 std::cout <<
"Max mod = " << max_mod << std::endl;
191 if (print_inoutflows_) {
192 std::pair<double, double> w_io, o_io;
193 computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
194 std::cout <<
"Pressure step " << iter
195 <<
"\nWater flow [in] " << w_io.first
196 <<
" [out] " << w_io.second
197 <<
"\nOil flow [in] " << o_io.first
198 <<
" [out] " << o_io.second
206 this->flow_solver_.getSolution(),
208 std::string(
"output-steadystate")
209 +
'-' + std::to_string(count)
210 +
'-' + std::to_string(flow_direction)
211 +
'-' + std::to_string(iter));
215 int num_cells = saturation.size();
216 double maxdiff = 0.0;
217 for (
int i = 0; i < num_cells; ++i) {
218 maxdiff = std::max(maxdiff, std::fabs(saturation[i] - saturation_old[i]));
221 std::cout <<
"Maximum saturation change: " << maxdiff << std::endl;
223 if (maxdiff < sat_change_threshold_) {
225 std::cout <<
"Maximum saturation change is under steady state threshold." << std::endl;
231 saturation_old = saturation;
236 typedef typename Super::ResProp::Mobility Mob;
240 for (
int c = 0; c < num_cells; ++c) {
241 this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
242 m1max = maxMobility(m1max, m.mob);
243 this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
244 m2max = maxMobility(m2max, m.mob);
247 const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
248 const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
249 const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
250 const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
251 const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
252 const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
254 std::vector<Mob> mob1(num_cells);
255 std::vector<Mob> mob2(num_cells);
256 for (
int c = 0; c < num_cells; ++c) {
257 this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
259 this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
265 permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
267 permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
270 last_saturation_state_.swap(saturation);
281 k_rw *= this->res_prop_.viscosityFirstPhase();
283 k_ro *= this->res_prop_.viscositySecondPhase();
284 return std::make_pair(k_rw, k_ro);
290 template <
class Traits>
291 inline const std::vector<double>&
294 return last_saturation_state_;
300 template <
class Traits>
304 double pore_vol = 0.0;
305 double sat_vol = 0.0;
306 for (
CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
307 double cell_pore_vol = c->
volume()*this->res_prop_.porosity(c->index());
308 pore_vol += cell_pore_vol;
309 sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
312 return sat_vol/pore_vol;
317 std::map< int, double >& frac_flow_by_bid,
319 const ResProp& res_prop,
320 GridInterface::CellIterator c,
321 CellIterator::FaceIterator f ) {
323 const SatBC& sc = bcond.satCond( f );
327 return res_prop.fractionalFlow( c->index(), sc.
saturation() );
331 const int partner_bid = bcond.getPeriodicPartner( f->boundaryId() );
332 const auto it = frac_flow_by_bid.find( partner_bid );
336 if( it == frac_flow_by_bid.end() ) {
337 OPM_THROW(std::runtime_error,
"Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
338 <<
" and partner bid = " << partner_bid);
344 template <
class Traits>
345 template <
class FlowSol>
347 std::pair<double, double>& oil_inout,
348 const FlowSol& flow_solution,
349 const std::vector<double>& saturations)
const
356 sideflux& operator+=(
const sideflux& other ) {
357 this->water += other.water;
358 this->oil += other.oil;
363 std::map< int, double > frac_flow_by_bid;
366 for(
auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
367 const double frac_flow = this->res_prop_.fractional_flow( c->index(), saturations[ c->index() ] );
369 for(
auto f = c->facebegin(); f != c->faceend(); ++f ) {
370 if( !f->boundary() )
continue;
372 const double flux = flow_solution.outflux( f );
373 if( flux < 0.0 )
continue;
375 const SatBC& sc = this->bcond_.satCond( f );
376 if( sc.
isPeriodic() ) frac_flow_by_bid[ f->boundaryId() ] = frac_flow;
378 outflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
384 for(
auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
385 for(
auto f = c->facebegin(); f != c->faceend(); ++f ) {
386 if( !f->boundary() )
continue;
388 const double flux = flow_solution.outflux( f );
389 if( flux >= 0.0 )
continue;
391 const double frac_flow =
calc_frac_flow( frac_flow_by_bid, this->bcond_, this->res_prop_, c, f );
393 inflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
397 water_inout = std::make_pair( in.water, out.water );
398 oil_inout = std::make_pair( in.oil, out.oil );
bool isPeriodic() const
Type query.
Definition: BoundaryConditions.hpp:97
bool isDirichlet() const
Type query.
Definition: BoundaryConditions.hpp:85
Definition: GridInterfaceEuler.hpp:368
Scalar volume() const
Definition: GridInterfaceEuler.hpp:337
Definition: ReservoirPropertyFixedMobility.hpp:46
A class for representing a saturation boundary condition.
Definition: BoundaryConditions.hpp:176
double saturation() const
Query a Dirichlet condition.
Definition: BoundaryConditions.hpp:197
double saturationDifference() const
Query a Periodic condition.
Definition: BoundaryConditions.hpp:204
virtual void initImpl(const Opm::parameter::ParameterGroup ¶m)
Override from superclass.
Definition: SteadyStateUpscaler_impl.hpp:70
const std::vector< double > & lastSaturationState() const
Definition: SteadyStateUpscaler_impl.hpp:292
std::pair< permtensor_t, permtensor_t > upscaleSteadyState(const int flow_direction, const std::vector< double > &initial_saturation, const double boundary_saturation, const double pressure_drop, const permtensor_t &upscaled_perm)
Definition: SteadyStateUpscaler_impl.hpp:132
double lastSaturationUpscaled() const
Definition: SteadyStateUpscaler_impl.hpp:301
Super::permtensor_t permtensor_t
Definition: SteadyStateUpscaler.hpp:57
SteadyStateUpscaler()
Default constructor.
Definition: SteadyStateUpscaler_impl.hpp:53
void computeInOutFlows(std::pair< double, double > &water_inout, std::pair< double, double > &oil_inout, const FlowSol &flow_solution, const std::vector< double > &saturations) const
Definition: SteadyStateUpscaler_impl.hpp:346
A base class for upscaling.
Definition: UpscalerBase.hpp:56
min[0]
Definition: elasticity_upscale_impl.hpp:146
void thresholdMobility(double &m, double threshold)
Definition: ImplicitAssembly.hpp:43
M matprod(const M &m1, const M &m2)
Definition: MatrixInverse.hpp:69
M inverse3x3(const M &m)
Definition: MatrixInverse.hpp:86
static double calc_frac_flow(std::map< int, double > &frac_flow_by_bid, const BCs &bcond, const ResProp &res_prop, GridInterface::CellIterator c, CellIterator::FaceIterator f)
Definition: SteadyStateUpscaler_impl.hpp:316
void writeVtkOutput(const GridInterface &ginterf, const ReservoirProperties &rp, const FlowSol &flowsol, const std::vector< double > &saturation, const std::string &filename)
Definition: SimulatorUtilities.hpp:239
void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition: setupBoundaryConditions.hpp:99