36#ifndef OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
37#define OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
38#include <boost/date_time/posix_time/posix_time.hpp>
44#include <opm/input/eclipse/Units/Units.hpp>
55 template <
class Traits>
61 print_inoutflows_(false),
62 simulation_steps_(10),
64 relperm_threshold_(1.0e-8),
65 maximum_mobility_contrast_(1.0e9),
66 sat_change_year_(1.0e-5),
77 template <
class Traits>
80 Super::initImpl(param);
81 use_gravity_ = param.getDefault(
"use_gravity", use_gravity_);
82 output_vtk_ = param.getDefault(
"output_vtk", output_vtk_);
83 output_ecl_ = param.getDefault(
"output_ecl", output_ecl_);
85 grid_adapter_.init(Super::grid());
87 print_inoutflows_ = param.getDefault(
"print_inoutflows", print_inoutflows_);
88 simulation_steps_ = param.getDefault(
"simulation_steps", simulation_steps_);
89 init_stepsize_ = Opm::unit::convert::from(param.getDefault(
"init_stepsize", init_stepsize_),
91 relperm_threshold_ = param.getDefault(
"relperm_threshold", relperm_threshold_);
92 maximum_mobility_contrast_ = param.getDefault(
"maximum_mobility_contrast", maximum_mobility_contrast_);
93 sat_change_year_ = param.getDefault(
"sat_change_year", sat_change_year_);
94 dt_sat_tol_ = param.getDefault(
"dt_sat_tol", dt_sat_tol_);
95 max_it_ = param.getDefault(
"max_it", max_it_);
96 max_stepsize_ = Opm::unit::convert::from(param.getDefault(
"max_stepsize", max_stepsize_),Opm::unit::year);
97 use_maxdiff_ = param.getDefault(
"use_maxdiff", use_maxdiff_);
98 transport_solver_.init(param);
100 double v1_default = this->res_prop_.viscosityFirstPhase();
101 double v2_default = this->res_prop_.viscositySecondPhase();
102 this->res_prop_.setViscosities(param.getDefault(
"viscosity1", v1_default), param.getDefault(
"viscosity2", v2_default));
103 double d1_default = this->res_prop_.densityFirstPhase();
104 double d2_default = this->res_prop_.densitySecondPhase();
105 this->res_prop_.setDensities(param.getDefault(
"density1", d1_default), param.getDefault(
"density2", d2_default));
112 double maxMobility(
double m1,
double m2)
114 return std::max(m1, m2);
117 template <
class SomeMatrixType>
118 double maxMobility(
double m1, SomeMatrixType& m2)
121 for (
int i = 0; i <
std::min(m2.numRows(), m2.numCols()); ++i) {
122 m = std::max(m, m2(i,i));
128 m = std::max(m, threshold);
131 template <
class SomeMatrixType>
134 for (
int i = 0; i <
std::min(m.numRows(), m.numCols()); ++i) {
135 m(i, i) = std::max(m(i, i), threshold);
140 inline std::vector< double > destripe(
const std::vector< double >& v,
144 std::vector< double > dst( v.size() / stride );
146 for(
size_t i = offset; i < v.size(); i += stride ) {
147 dst[ di++ ] = v[ i ];
154 void waterSatToBothSat(
const std::vector<double>& sw,
155 std::vector<double>& sboth)
159 for (
int i = 0; i < num; ++i) {
161 sboth[2*i + 1] = 1.0 - sw[i];
171 template <
class Traits>
172 inline std::pair<typename SteadyStateUpscalerImplicit<Traits>::permtensor_t,
176 const std::vector<double>& initial_saturation,
177 const double boundary_saturation,
178 const double pressure_drop,
182 static int count = 0;
184 int num_cells = this->ginterf_.numberOfCells();
186 std::vector<double> src(num_cells, 0.0);
187 Opm::SparseVector<double> injection(num_cells);
189 Dune::FieldVector<double, 3> gravity(0.0);
191 gravity[2] = Opm::unit::gravity;
194 if (gravity.two_norm() > 0.0) {
195 OPM_MESSAGE(
"Warning: Gravity is experimental for flow solver.");
199 std::vector<double> pore_vol;
200 pore_vol.reserve(num_cells);
201 double tot_pore_vol = 0.0;
203 for (CellIterT c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
204 double cell_pore_vol = c->
volume()*this->res_prop_.porosity(c->index());
205 pore_vol.push_back(cell_pore_vol);
206 tot_pore_vol += cell_pore_vol;
210 std::vector<double> saturation = initial_saturation;
214 pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
217 if (flow_direction == 0) {
218 this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
220 transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
223 this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
224 this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_);
225 double max_mod = this->flow_solver_.postProcessFluxes();
226 std::cout <<
"Max mod = " << max_mod << std::endl;
229 std::vector<double> saturation_old = saturation;
230 bool stationary =
false;
232 double stepsize = init_stepsize_;
233 double ecl_time = 0.0;
234 std::vector<double> ecl_sat;
235 std::vector<double> ecl_press;
236 std::vector<double> init_saturation(saturation);
237 while ((!stationary) && (it_count < max_it_)) {
239 std::cout <<
"Running transport step " << it_count <<
" with stepsize "
240 << stepsize/Opm::unit::year <<
" years." << std::endl;
241 bool converged = transport_solver_.transportSolve(saturation, stepsize, gravity,
242 this->flow_solver_.getSolution(), injection);
245 init_saturation = saturation;
250 this->flow_solver_.postProcessFluxes();
252 if (print_inoutflows_) {
253 std::pair<double, double> w_io, o_io;
254 computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
255 std::cout <<
"Pressure step " << it_count
256 <<
"\nWater flow [in] " << w_io.first
257 <<
" [out] " << w_io.second
258 <<
"\nOil flow [in] " << o_io.first
259 <<
" [out] " << o_io.second
266 this->flow_solver_.getSolution(),
268 std::string(
"output-steadystate")
269 +
'-' + std::to_string(count)
270 +
'-' + std::to_string(flow_direction)
271 +
'-' + std::to_string(it_count));
274 const char* fd =
"xyz";
275 std::string basename = std::string(
"ecldump-steadystate")
276 +
'-' + std::to_string(boundary_saturation)
277 +
'-' + std::to_string(fd[flow_direction])
278 +
'-' + std::to_string(pressure_drop);
279 waterSatToBothSat(saturation, ecl_sat);
280 getCellPressure(ecl_press, this->ginterf_, this->flow_solver_.getSolution());
281 data::Solution solution;
282 solution.insert(
"PRESSURE" , UnitSystem::measure::pressure , ecl_press , data::TargetType::RESTART_SOLUTION );
283 solution.insert(
"SWAT" , UnitSystem::measure::identity , destripe( ecl_sat, 2, 0) , data::TargetType::RESTART_SOLUTION );
284 ecl_time += stepsize;
285 boost::posix_time::ptime ecl_startdate( boost::gregorian::date(2012, 1, 1) );
286 boost::posix_time::ptime ecl_curdate = ecl_startdate + boost::posix_time::seconds(
int(ecl_time));
287 boost::posix_time::ptime epoch( boost::gregorian::date( 1970, 1, 1 ) );
288 auto ecl_posix_time = ( ecl_curdate - epoch ).total_seconds();
289 const auto* cgrid = grid_adapter_.c_grid();
291 cgrid->cartdims[ 1 ],
292 cgrid->cartdims[ 2 ],
293 cgrid->number_of_cells,
295 ecl_time, ecl_posix_time,
299 double maxdiff = 0.0;
300 double euclidean_diff = 0.0;
301 for (
int i = 0; i < num_cells; ++i) {
302 const double sat_diff_cell = saturation[i] - saturation_old[i];
303 maxdiff = std::max(maxdiff, std::fabs(sat_diff_cell));
304 euclidean_diff += sat_diff_cell * sat_diff_cell * pore_vol[i];
306 euclidean_diff = std::sqrt(euclidean_diff / tot_pore_vol);
309 ds_year = maxdiff*Opm::unit::year/stepsize;
310 std::cout <<
"Maximum saturation change/year: " << ds_year << std::endl;
311 if (maxdiff < dt_sat_tol_) {
312 stepsize=
std::min(max_stepsize_,2*stepsize);
316 ds_year = euclidean_diff*Opm::unit::year/stepsize;
317 std::cout <<
"Euclidean saturation change/year: " << ds_year << std::endl;
318 if (euclidean_diff < dt_sat_tol_) {
319 stepsize=
std::min(max_stepsize_,2*stepsize);
322 if (ds_year < sat_change_year_) {
326 std::cerr <<
"Cutting time step\n";
327 init_saturation = saturation_old;
328 stepsize=stepsize/2.0;
332 saturation_old = saturation;
334 success = stationary;
338 typedef typename Super::ResProp::Mobility Mob;
342 for (
int c = 0; c < num_cells; ++c) {
343 this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
344 m1max = maxMobility(m1max, m.mob);
345 this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
346 m2max = maxMobility(m2max, m.mob);
349 const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
350 const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
351 const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
352 const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
353 const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
354 const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
356 std::vector<Mob> mob1(num_cells);
357 std::vector<Mob> mob2(num_cells);
358 for (
int c = 0; c < num_cells; ++c) {
359 this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
361 this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
367 permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
369 permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
372 last_saturation_state_.swap(saturation);
383 k_rw *= this->res_prop_.viscosityFirstPhase();
385 k_ro *= this->res_prop_.viscositySecondPhase();
386 return std::make_pair(k_rw, k_ro);
392 template <
class Traits>
393 inline const std::vector<double>&
396 return last_saturation_state_;
402 template <
class Traits>
406 double pore_vol = 0.0;
407 double sat_vol = 0.0;
408 for (CellIterT c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
409 double cell_pore_vol = c->
volume()*this->res_prop_.porosity(c->index());
410 pore_vol += cell_pore_vol;
411 sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
414 return sat_vol/pore_vol;
418 template <
class Traits>
421 for (
int c = 0; c < int (s.size()); c++ ) {
422 double s_min = this->res_prop_.s_min(c);
423 double s_max = this->res_prop_.s_max(c);
424 s[c] = std::max(s_min+1e-4, s[c]);
429 template <
class Traits>
432 int num_cells = this->ginterf_.numberOfCells();
433 std::vector<double> s_orig(num_cells, average_s);
435 initSatLimits(s_orig);
436 std::vector<double> cap_press(num_cells, 0.0);
439 double cap_press_range = 1e2;
440 double mod_low = 1e100;
441 double mod_high = -1e100;
442 Opm::bracketZero(func, 0.0, cap_press_range, mod_low, mod_high);
443 const int max_iter = 40;
444 const double nonlinear_tolerance = 1e-12;
445 int iterations_used = -1;
446 typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
447 double mod_correct = RootFinder::solve(func, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
448 std::cout <<
"Moved capillary pressure solution by " << mod_correct <<
" after "
449 << iterations_used <<
" iterations." << std::endl;
454 template <
class Traits>
455 template <
class FlowSol>
457 std::pair<double, double>& oil_inout,
458 const FlowSol& flow_solution,
459 const std::vector<double>& saturations)
const
462 typedef typename CellIterT::FaceIterator FaceIterT;
464 double side1_flux = 0.0;
465 double side2_flux = 0.0;
466 double side1_flux_oil = 0.0;
467 double side2_flux_oil = 0.0;
468 std::map<int, double> frac_flow_by_bid;
469 int num_cells = this->ginterf_.numberOfCells();
470 std::vector<double> cell_inflows_w(num_cells, 0.0);
471 std::vector<double> cell_outflows_w(num_cells, 0.0);
476 for (
int pass = 0; pass < 2; ++pass) {
477 for (CellIterT c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
478 for (FaceIterT f = c->facebegin(); f != c->faceend(); ++f) {
480 double flux = flow_solution.outflux(f);
481 const SatBC& sc = this->bcond_.satCond(f);
482 if (flux < 0.0 && pass == 1) {
484 double frac_flow = 0.0;
487 int partner_bid = this->bcond_.getPeriodicPartner(f->boundaryId());
488 std::map<int, double>::const_iterator it = frac_flow_by_bid.find(partner_bid);
489 if (it == frac_flow_by_bid.end()) {
490 OPM_THROW(std::runtime_error,
491 "Could not find periodic partner fractional flow. "
492 "Face bid = " + std::to_string(f->boundaryId()) +
493 " and partner bid = " + std::to_string(partner_bid));
495 frac_flow = it->second;
498 frac_flow = this->res_prop_.fractionalFlow(c->index(), sc.
saturation());
500 cell_inflows_w[c->index()] += flux*frac_flow;
501 side1_flux += flux*frac_flow;
502 side1_flux_oil += flux*(1.0 - frac_flow);
503 }
else if (flux >= 0.0 && pass == 0) {
505 double frac_flow = this->res_prop_.fractionalFlow(c->index(), saturations[c->index()]);
507 frac_flow_by_bid[f->boundaryId()] = frac_flow;
510 cell_outflows_w[c->index()] += flux*frac_flow;
511 side2_flux += flux*frac_flow;
512 side2_flux_oil += flux*(1.0 - frac_flow);
518 water_inout = std::make_pair(side1_flux, side2_flux);
519 oil_inout = std::make_pair(side1_flux_oil, side2_flux_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
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: SteadyStateUpscalerImplicit_impl.hpp:456
void initSatLimits(std::vector< double > &s) const
Ensure saturations are not outside table.
Definition: SteadyStateUpscalerImplicit_impl.hpp:419
Super::permtensor_t permtensor_t
Definition: SteadyStateUpscalerImplicit.hpp:58
virtual void initImpl(const Opm::ParameterGroup ¶m)
Override from superclass.
Definition: SteadyStateUpscalerImplicit_impl.hpp:78
double lastSaturationUpscaled() const
Definition: SteadyStateUpscalerImplicit_impl.hpp:403
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, bool &success)
Definition: SteadyStateUpscalerImplicit_impl.hpp:175
void setToCapillaryLimit(double average_s, std::vector< double > &s) const
Definition: SteadyStateUpscalerImplicit_impl.hpp:430
SteadyStateUpscalerImplicit()
Default constructor.
Definition: SteadyStateUpscalerImplicit_impl.hpp:56
const std::vector< double > & lastSaturationState() const
Definition: SteadyStateUpscalerImplicit_impl.hpp:394
A base class for upscaling.
Definition: UpscalerBase.hpp:56
Traits::template ResProp< Dimension >::Type ResProp
Definition: UpscalerBase.hpp:63
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
void writeECLData(int nx, int ny, int nz, int number_of_cells, data::Solution, const int current_step, const double current_time, time_t current_posix_time, const std::string &output_dir, const std::string &base_name)
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
void getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition: SimulatorUtilities.hpp:184
Definition: MatchSaturatedVolumeFunctor.hpp:68
const std::vector< double > & lastSaturations() const
Definition: MatchSaturatedVolumeFunctor.hpp:96