36 #ifndef OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER 37 #define OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER 40 #include <boost/lexical_cast.hpp> 41 #include <opm/porsol/common/MatrixInverse.hpp> 42 #include <opm/porsol/common/SimulatorUtilities.hpp> 43 #include <opm/porsol/common/ReservoirPropertyFixedMobility.hpp> 44 #include <opm/porsol/euler/MatchSaturatedVolumeFunctor.hpp> 45 #include <opm/parser/eclipse/Units/Units.hpp> 47 #include <opm/upscaling/writeECLData.hpp> 48 #include <opm/core/utility/miscUtilities.hpp> 57 template <
class Traits>
63 print_inoutflows_(false),
64 simulation_steps_(10),
66 relperm_threshold_(1.0e-8),
67 maximum_mobility_contrast_(1.0e9),
68 sat_change_year_(1.0e-5),
79 template <
class Traits>
82 Super::initImpl(param);
83 use_gravity_ = param.getDefault(
"use_gravity", use_gravity_);
84 output_vtk_ = param.getDefault(
"output_vtk", output_vtk_);
85 output_ecl_ = param.getDefault(
"output_ecl", output_ecl_);
87 grid_adapter_.init(Super::grid());
89 print_inoutflows_ = param.getDefault(
"print_inoutflows", print_inoutflows_);
90 simulation_steps_ = param.getDefault(
"simulation_steps", simulation_steps_);
91 init_stepsize_ = Opm::unit::convert::from(param.getDefault(
"init_stepsize", init_stepsize_),
93 relperm_threshold_ = param.getDefault(
"relperm_threshold", relperm_threshold_);
94 maximum_mobility_contrast_ = param.getDefault(
"maximum_mobility_contrast", maximum_mobility_contrast_);
95 sat_change_year_ = param.getDefault(
"sat_change_year", sat_change_year_);
96 dt_sat_tol_ = param.getDefault(
"dt_sat_tol", dt_sat_tol_);
97 max_it_ = param.getDefault(
"max_it", max_it_);
98 max_stepsize_ = Opm::unit::convert::from(param.getDefault(
"max_stepsize", max_stepsize_),Opm::unit::year);
99 use_maxdiff_ = param.getDefault(
"use_maxdiff", use_maxdiff_);
100 transport_solver_.init(param);
102 double v1_default = this->res_prop_.viscosityFirstPhase();
103 double v2_default = this->res_prop_.viscositySecondPhase();
104 this->res_prop_.setViscosities(param.getDefault(
"viscosity1", v1_default), param.getDefault(
"viscosity2", v2_default));
105 double d1_default = this->res_prop_.densityFirstPhase();
106 double d2_default = this->res_prop_.densitySecondPhase();
107 this->res_prop_.setDensities(param.getDefault(
"density1", d1_default), param.getDefault(
"density2", d2_default));
114 double maxMobility(
double m1,
double m2)
116 return std::max(m1, m2);
119 template <
class SomeMatrixType>
120 double maxMobility(
double m1, SomeMatrixType& m2)
123 for (
int i = 0; i < std::min(m2.numRows(), m2.numCols()); ++i) {
124 m = std::max(m, m2(i,i));
128 void thresholdMobility(
double& m,
double threshold)
130 m = std::max(m, threshold);
133 template <
class SomeMatrixType>
134 void thresholdMobility(SomeMatrixType& m,
double threshold)
136 for (
int i = 0; i < std::min(m.numRows(), m.numCols()); ++i) {
137 m(i, i) = std::max(m(i, i), threshold);
142 inline std::vector< double > destripe(
const std::vector< double >& v,
146 std::vector< double > dst( v.size() / stride );
148 for(
size_t i = offset; i < v.size(); i += stride ) {
149 dst[ di++ ] = v[ i ];
159 template <
class Traits>
160 inline std::pair<typename SteadyStateUpscalerImplicit<Traits>::permtensor_t,
161 typename SteadyStateUpscalerImplicit<Traits>::permtensor_t>
164 const std::vector<double>& initial_saturation,
165 const double boundary_saturation,
166 const double pressure_drop,
167 const permtensor_t& upscaled_perm,
170 static int count = 0;
172 int num_cells = this->ginterf_.numberOfCells();
174 std::vector<double> src(num_cells, 0.0);
175 Opm::SparseVector<double> injection(num_cells);
177 Dune::FieldVector<double, 3> gravity(0.0);
179 gravity[2] = Opm::unit::gravity;
182 if (gravity.two_norm() > 0.0) {
183 OPM_MESSAGE(
"Warning: Gravity is experimental for flow solver.");
187 std::vector<double> pore_vol;
188 pore_vol.reserve(num_cells);
189 double tot_pore_vol = 0.0;
191 for (
CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
192 double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
193 pore_vol.push_back(cell_pore_vol);
194 tot_pore_vol += cell_pore_vol;
198 std::vector<double> saturation = initial_saturation;
202 pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
205 if (flow_direction == 0) {
206 this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
208 transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
211 this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
212 this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_);
213 double max_mod = this->flow_solver_.postProcessFluxes();
214 std::cout <<
"Max mod = " << max_mod << std::endl;
217 std::vector<double> saturation_old = saturation;
218 bool stationary =
false;
220 double stepsize = init_stepsize_;
221 double ecl_time = 0.0;
222 std::vector<double> ecl_sat;
223 std::vector<double> ecl_press;
224 std::vector<double> init_saturation(saturation);
225 while ((!stationary) && (it_count < max_it_)) {
227 std::cout <<
"Running transport step " << it_count <<
" with stepsize " 228 << stepsize/Opm::unit::year <<
" years." << std::endl;
229 bool converged = transport_solver_.transportSolve(saturation, stepsize, gravity,
230 this->flow_solver_.getSolution(), injection);
233 init_saturation = saturation;
238 this->flow_solver_.postProcessFluxes();
240 if (print_inoutflows_) {
241 std::pair<double, double> w_io, o_io;
242 computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
243 std::cout <<
"Pressure step " << it_count
244 <<
"\nWater flow [in] " << w_io.first
245 <<
" [out] " << w_io.second
246 <<
"\nOil flow [in] " << o_io.first
247 <<
" [out] " << o_io.second
252 writeVtkOutput(this->ginterf_,
254 this->flow_solver_.getSolution(),
256 std::string(
"output-steadystate")
257 +
'-' + boost::lexical_cast<std::string>(count)
258 +
'-' + boost::lexical_cast<std::string>(flow_direction)
259 +
'-' + boost::lexical_cast<std::string>(it_count));
262 const char* fd =
"xyz";
263 std::string basename = std::string(
"ecldump-steadystate")
264 +
'-' + boost::lexical_cast<std::string>(boundary_saturation)
265 +
'-' + boost::lexical_cast<std::string>(fd[flow_direction])
266 +
'-' + boost::lexical_cast<std::string>(pressure_drop);
267 Opm::toBothSat(saturation, ecl_sat);
268 getCellPressure(ecl_press, this->ginterf_, this->flow_solver_.getSolution());
269 data::Solution solution;
270 solution.insert(
"PRESSURE" , UnitSystem::measure::pressure , ecl_press , data::TargetType::RESTART_SOLUTION );
271 solution.insert(
"SWAT" , UnitSystem::measure::identity , destripe( ecl_sat, 2, 0) , data::TargetType::RESTART_SOLUTION );
272 ecl_time += stepsize;
273 boost::posix_time::ptime ecl_startdate( boost::gregorian::date(2012, 1, 1) );
274 boost::posix_time::ptime ecl_curdate = ecl_startdate + boost::posix_time::seconds(
int(ecl_time));
275 boost::posix_time::ptime epoch( boost::gregorian::date( 1970, 1, 1 ) );
276 auto ecl_posix_time = ( ecl_curdate - epoch ).total_seconds();
277 const auto* cgrid = grid_adapter_.c_grid();
278 Opm::writeECLData(cgrid->cartdims[ 0 ],
279 cgrid->cartdims[ 1 ],
280 cgrid->cartdims[ 2 ],
281 cgrid->number_of_cells,
283 ecl_time, ecl_posix_time,
287 double maxdiff = 0.0;
288 double euclidean_diff = 0.0;
289 for (
int i = 0; i < num_cells; ++i) {
290 const double sat_diff_cell = saturation[i] - saturation_old[i];
291 maxdiff = std::max(maxdiff, std::fabs(sat_diff_cell));
292 euclidean_diff += sat_diff_cell * sat_diff_cell * pore_vol[i];
294 euclidean_diff = std::sqrt(euclidean_diff / tot_pore_vol);
297 ds_year = maxdiff*Opm::unit::year/stepsize;
298 std::cout <<
"Maximum saturation change/year: " << ds_year << std::endl;
299 if (maxdiff < dt_sat_tol_) {
300 stepsize=std::min(max_stepsize_,2*stepsize);
304 ds_year = euclidean_diff*Opm::unit::year/stepsize;
305 std::cout <<
"Euclidean saturation change/year: " << ds_year << std::endl;
306 if (euclidean_diff < dt_sat_tol_) {
307 stepsize=std::min(max_stepsize_,2*stepsize);
310 if (ds_year < sat_change_year_) {
314 std::cerr <<
"Cutting time step\n";
315 init_saturation = saturation_old;
316 stepsize=stepsize/2.0;
320 saturation_old = saturation;
322 success = stationary;
326 typedef typename Super::ResProp::Mobility Mob;
330 for (
int c = 0; c < num_cells; ++c) {
331 this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
332 m1max = maxMobility(m1max, m.mob);
333 this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
334 m2max = maxMobility(m2max, m.mob);
337 const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
338 const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
339 const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
340 const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
341 const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
342 const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
344 std::vector<Mob> mob1(num_cells);
345 std::vector<Mob> mob2(num_cells);
346 for (
int c = 0; c < num_cells; ++c) {
347 this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
348 thresholdMobility(mob1[c].mob, mob1_threshold);
349 this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
350 thresholdMobility(mob2[c].mob, mob2_threshold);
355 permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
357 permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
360 last_saturation_state_.swap(saturation);
365 permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
366 permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
370 permtensor_t k_rw(lambda_w);
371 k_rw *= this->res_prop_.viscosityFirstPhase();
372 permtensor_t k_ro(lambda_o);
373 k_ro *= this->res_prop_.viscositySecondPhase();
374 return std::make_pair(k_rw, k_ro);
380 template <
class Traits>
381 inline const std::vector<double>&
384 return last_saturation_state_;
390 template <
class Traits>
394 double pore_vol = 0.0;
395 double sat_vol = 0.0;
396 for (
CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
397 double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
398 pore_vol += cell_pore_vol;
399 sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
402 return sat_vol/pore_vol;
406 template <
class Traits>
409 for (
int c = 0; c < int (s.size()); c++ ) {
410 double s_min = this->res_prop_.s_min(c);
411 double s_max = this->res_prop_.s_max(c);
412 s[c] = std::max(s_min+1e-4, s[c]);
413 s[c] = std::min(s_max-1e-4, s[c]);
417 template <
class Traits>
420 int num_cells = this->ginterf_.numberOfCells();
421 std::vector<double> s_orig(num_cells, average_s);
423 initSatLimits(s_orig);
424 std::vector<double> cap_press(num_cells, 0.0);
425 typedef typename UpscalerBase<Traits>::ResProp ResProp;
427 double cap_press_range = 1e2;
428 double mod_low = 1e100;
429 double mod_high = -1e100;
430 Opm::bracketZero(func, 0.0, cap_press_range, mod_low, mod_high);
431 const int max_iter = 40;
432 const double nonlinear_tolerance = 1e-12;
433 int iterations_used = -1;
434 typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
435 double mod_correct = RootFinder::solve(func, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
436 std::cout <<
"Moved capillary pressure solution by " << mod_correct <<
" after " 437 << iterations_used <<
" iterations." << std::endl;
438 s = func.lastSaturations();
442 template <
class Traits>
443 template <
class FlowSol>
444 void SteadyStateUpscalerImplicit<Traits>::computeInOutFlows(std::pair<double, double>& water_inout,
445 std::pair<double, double>& oil_inout,
446 const FlowSol& flow_solution,
447 const std::vector<double>& saturations)
const 449 typedef typename GridInterface::CellIterator CellIter;
450 typedef typename CellIter::FaceIterator FaceIter;
452 double side1_flux = 0.0;
453 double side2_flux = 0.0;
454 double side1_flux_oil = 0.0;
455 double side2_flux_oil = 0.0;
456 std::map<int, double> frac_flow_by_bid;
457 int num_cells = this->ginterf_.numberOfCells();
458 std::vector<double> cell_inflows_w(num_cells, 0.0);
459 std::vector<double> cell_outflows_w(num_cells, 0.0);
464 for (
int pass = 0; pass < 2; ++pass) {
465 for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
466 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
468 double flux = flow_solution.outflux(f);
469 const SatBC& sc = this->bcond_.satCond(f);
470 if (flux < 0.0 && pass == 1) {
472 double frac_flow = 0.0;
473 if (sc.isPeriodic()) {
474 assert(sc.saturationDifference() == 0.0);
475 int partner_bid = this->bcond_.getPeriodicPartner(f->boundaryId());
476 std::map<int, double>::const_iterator it = frac_flow_by_bid.find(partner_bid);
477 if (it == frac_flow_by_bid.end()) {
478 OPM_THROW(std::runtime_error,
"Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
479 <<
" and partner bid = " << partner_bid);
481 frac_flow = it->second;
483 assert(sc.isDirichlet());
484 frac_flow = this->res_prop_.fractionalFlow(c->index(), sc.saturation());
486 cell_inflows_w[c->index()] += flux*frac_flow;
487 side1_flux += flux*frac_flow;
488 side1_flux_oil += flux*(1.0 - frac_flow);
489 }
else if (flux >= 0.0 && pass == 0) {
491 double frac_flow = this->res_prop_.fractionalFlow(c->index(), saturations[c->index()]);
492 if (sc.isPeriodic()) {
493 frac_flow_by_bid[f->boundaryId()] = frac_flow;
496 cell_outflows_w[c->index()] += flux*frac_flow;
497 side2_flux += flux*frac_flow;
498 side2_flux_oil += flux*(1.0 - frac_flow);
504 water_inout = std::make_pair(side1_flux, side2_flux);
505 oil_inout = std::make_pair(side1_flux_oil, side2_flux_oil);
514 #endif // OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition: setupBoundaryConditions.hpp:99
Definition: MatchSaturatedVolumeFunctor.hpp:68
Definition: GridInterfaceEuler.hpp:376
A base class for upscaling.
Definition: UpscalerBase.hpp:55
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
double lastSaturationUpscaled() const
Computes the upscaled saturation corresponding to the saturation field returned by lastSaturationStat...
Definition: SteadyStateUpscalerImplicit_impl.hpp:391
void getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition: SimulatorUtilities.hpp:184
const std::vector< double > & lastSaturationState() const
Accessor for the steady state saturation field.
Definition: SteadyStateUpscalerImplicit_impl.hpp:382
SteadyStateUpscalerImplicit()
Default constructor.
Definition: SteadyStateUpscalerImplicit_impl.hpp:58
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)
Does a steady-state upscaling.
Definition: SteadyStateUpscalerImplicit_impl.hpp:163
virtual void initImpl(const Opm::ParameterGroup ¶m)
Override from superclass.
Definition: SteadyStateUpscalerImplicit_impl.hpp:80
void initSatLimits(std::vector< double > &s) const
Ensure saturations are not outside table.
Definition: SteadyStateUpscalerImplicit_impl.hpp:407
A class for doing steady state upscaling.
Definition: SteadyStateUpscalerImplicit.hpp:52
Definition: ReservoirPropertyFixedMobility.hpp:46