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/core/utility/Units.hpp>
46 #include <opm/core/io/eclipse/writeECLData.hpp>
47 #include <opm/core/utility/miscUtilities.hpp>
56 template <
class Traits>
62 print_inoutflows_(false),
63 simulation_steps_(10),
65 relperm_threshold_(1.0e-8),
66 maximum_mobility_contrast_(1.0e9),
67 sat_change_year_(1.0e-5),
78 template <
class Traits>
81 Super::initImpl(param);
82 use_gravity_ = param.getDefault(
"use_gravity", use_gravity_);
83 output_vtk_ = param.getDefault(
"output_vtk", output_vtk_);
84 output_ecl_ = param.getDefault(
"output_ecl", output_ecl_);
86 grid_adapter_.init(Super::grid());
88 print_inoutflows_ = param.getDefault(
"print_inoutflows", print_inoutflows_);
89 simulation_steps_ = param.getDefault(
"simulation_steps", simulation_steps_);
90 init_stepsize_ = Opm::unit::convert::from(param.getDefault(
"init_stepsize", init_stepsize_),
92 relperm_threshold_ = param.getDefault(
"relperm_threshold", relperm_threshold_);
93 maximum_mobility_contrast_ = param.getDefault(
"maximum_mobility_contrast", maximum_mobility_contrast_);
94 sat_change_year_ = param.getDefault(
"sat_change_year", sat_change_year_);
95 dt_sat_tol_ = param.getDefault(
"dt_sat_tol", dt_sat_tol_);
96 max_it_ = param.getDefault(
"max_it", max_it_);
97 max_stepsize_ = Opm::unit::convert::from(param.getDefault(
"max_stepsize", max_stepsize_),Opm::unit::year);
98 use_maxdiff_ = param.getDefault(
"use_maxdiff", use_maxdiff_);
99 transport_solver_.init(param);
101 double v1_default = this->res_prop_.viscosityFirstPhase();
102 double v2_default = this->res_prop_.viscositySecondPhase();
103 this->res_prop_.setViscosities(param.getDefault(
"viscosity1", v1_default), param.getDefault(
"viscosity2", v2_default));
104 double d1_default = this->res_prop_.densityFirstPhase();
105 double d2_default = this->res_prop_.densitySecondPhase();
106 this->res_prop_.setDensities(param.getDefault(
"density1", d1_default), param.getDefault(
"density2", d2_default));
113 double maxMobility(
double m1,
double m2)
115 return std::max(m1, m2);
118 template <
class SomeMatrixType>
119 double maxMobility(
double m1, SomeMatrixType& m2)
122 for (
int i = 0; i <
std::min(m2.numRows(), m2.numCols()); ++i) {
123 m = std::max(m, m2(i,i));
127 void thresholdMobility(
double& m,
double threshold)
129 m = std::max(m, threshold);
132 template <
class SomeMatrixType>
133 void thresholdMobility(SomeMatrixType& m,
double threshold)
135 for (
int i = 0; i <
std::min(m.numRows(), m.numCols()); ++i) {
136 m(i, i) = std::max(m(i, i), threshold);
144 template <
class Traits>
145 inline std::pair<typename SteadyStateUpscalerImplicit<Traits>::permtensor_t,
146 typename SteadyStateUpscalerImplicit<Traits>::permtensor_t>
149 const std::vector<double>& initial_saturation,
150 const double boundary_saturation,
151 const double pressure_drop,
155 static int count = 0;
157 int num_cells = this->ginterf_.numberOfCells();
159 std::vector<double> src(num_cells, 0.0);
160 Opm::SparseVector<double> injection(num_cells);
162 Dune::FieldVector<double, 3> gravity(0.0);
164 gravity[2] = Opm::unit::gravity;
167 if (gravity.two_norm() > 0.0) {
168 OPM_MESSAGE(
"Warning: Gravity is experimental for flow solver.");
172 std::vector<double> pore_vol;
173 pore_vol.reserve(num_cells);
174 double tot_pore_vol = 0.0;
175 typedef typename GridInterface::CellIterator
CellIter;
176 for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
177 double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
178 pore_vol.push_back(cell_pore_vol);
179 tot_pore_vol += cell_pore_vol;
183 std::vector<double> saturation = initial_saturation;
186 setupUpscalingConditions(this->ginterf_, this->bctype_, flow_direction,
187 pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
190 if (flow_direction == 0) {
191 this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
193 transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
196 this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
197 this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_);
198 double max_mod = this->flow_solver_.postProcessFluxes();
199 std::cout <<
"Max mod = " << max_mod << std::endl;
202 std::vector<double> saturation_old = saturation;
203 bool stationary =
false;
205 double stepsize = init_stepsize_;
206 double ecl_time = 0.0;
207 std::vector<double> ecl_sat;
208 std::vector<double> ecl_press;
209 std::vector<double> init_saturation(saturation);
210 while ((!stationary) && (it_count < max_it_)) {
212 std::cout <<
"Running transport step " << it_count <<
" with stepsize "
213 << stepsize/Opm::unit::year <<
" years." << std::endl;
214 bool converged = transport_solver_.transportSolve(saturation, stepsize, gravity,
215 this->flow_solver_.getSolution(), injection);
218 init_saturation = saturation;
223 this->flow_solver_.postProcessFluxes();
225 if (print_inoutflows_) {
226 std::pair<double, double> w_io, o_io;
227 computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
228 std::cout <<
"Pressure step " << it_count
229 <<
"\nWater flow [in] " << w_io.first
230 <<
" [out] " << w_io.second
231 <<
"\nOil flow [in] " << o_io.first
232 <<
" [out] " << o_io.second
237 writeVtkOutput(this->ginterf_,
239 this->flow_solver_.getSolution(),
241 std::string(
"output-steadystate")
242 +
'-' + boost::lexical_cast<std::string>(count)
243 +
'-' + boost::lexical_cast<std::string>(flow_direction)
244 +
'-' + boost::lexical_cast<std::string>(it_count));
247 const char* fd =
"xyz";
248 std::string basename = std::string(
"ecldump-steadystate")
249 +
'-' + boost::lexical_cast<std::string>(boundary_saturation)
250 +
'-' + boost::lexical_cast<std::string>(fd[flow_direction])
251 +
'-' + boost::lexical_cast<std::string>(pressure_drop);
252 Opm::toBothSat(saturation, ecl_sat);
253 getCellPressure(ecl_press, this->ginterf_, this->flow_solver_.getSolution());
254 Opm::DataMap datamap;
255 datamap[
"saturation"] = &ecl_sat;
256 datamap[
"pressure"] = &ecl_press;
257 ecl_time += stepsize;
258 boost::posix_time::ptime ecl_startdate( boost::gregorian::date(2012, 1, 1) );
259 boost::posix_time::ptime ecl_curdate = ecl_startdate + boost::posix_time::seconds(
int(ecl_time));
260 Opm::writeECLData(*grid_adapter_.c_grid(), datamap, it_count, ecl_time, ecl_curdate,
"./", basename);
263 int num_cells = saturation.size();
264 double maxdiff = 0.0;
265 double euclidean_diff = 0.0;
266 for (
int i = 0; i < num_cells; ++i) {
267 const double sat_diff_cell = saturation[i] - saturation_old[i];
268 maxdiff = std::max(maxdiff, std::fabs(sat_diff_cell));
269 euclidean_diff += sat_diff_cell * sat_diff_cell * pore_vol[i];
271 euclidean_diff = std::sqrt(euclidean_diff / tot_pore_vol);
274 ds_year = maxdiff*Opm::unit::year/stepsize;
275 std::cout <<
"Maximum saturation change/year: " << ds_year << std::endl;
276 if (maxdiff < dt_sat_tol_) {
277 stepsize=
std::min(max_stepsize_,2*stepsize);
281 ds_year = euclidean_diff*Opm::unit::year/stepsize;
282 std::cout <<
"Euclidean saturation change/year: " << ds_year << std::endl;
283 if (euclidean_diff < dt_sat_tol_) {
284 stepsize=
std::min(max_stepsize_,2*stepsize);
287 if (ds_year < sat_change_year_) {
291 std::cerr <<
"Cutting time step\n";
292 init_saturation = saturation_old;
293 stepsize=stepsize/2.0;
297 saturation_old = saturation;
299 success = stationary;
303 typedef typename Super::ResProp::Mobility Mob;
307 for (
int c = 0; c < num_cells; ++c) {
308 this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
309 m1max = maxMobility(m1max, m.mob);
310 this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
311 m2max = maxMobility(m2max, m.mob);
314 const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
315 const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
316 const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
317 const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
318 const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
319 const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
321 std::vector<Mob> mob1(num_cells);
322 std::vector<Mob> mob2(num_cells);
323 for (
int c = 0; c < num_cells; ++c) {
324 this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
325 thresholdMobility(mob1[c].mob, mob1_threshold);
326 this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
327 thresholdMobility(mob2[c].mob, mob2_threshold);
331 ReservoirPropertyFixedMobility<Mob> fluid_first(mob1);
332 permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
333 ReservoirPropertyFixedMobility<Mob> fluid_second(mob2);
334 permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
337 last_saturation_state_.swap(saturation);
342 permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
343 permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
348 k_rw *= this->res_prop_.viscosityFirstPhase();
350 k_ro *= this->res_prop_.viscositySecondPhase();
351 return std::make_pair(k_rw, k_ro);
357 template <
class Traits>
358 inline const std::vector<double>&
361 return last_saturation_state_;
367 template <
class Traits>
370 typedef typename GridInterface::CellIterator
CellIter;
371 double pore_vol = 0.0;
372 double sat_vol = 0.0;
373 for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
374 double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
375 pore_vol += cell_pore_vol;
376 sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
379 return sat_vol/pore_vol;
383 template <
class Traits>
386 for (
int c = 0; c < int (s.size()); c++ ) {
387 double s_min = this->res_prop_.s_min(c);
388 double s_max = this->res_prop_.s_max(c);
389 s[c] = std::max(s_min+1e-4, s[c]);
394 template <
class Traits>
397 int num_cells = this->ginterf_.numberOfCells();
398 std::vector<double> s_orig(num_cells, average_s);
400 initSatLimits(s_orig);
401 std::vector<double> cap_press(num_cells, 0.0);
403 MatchSaturatedVolumeFunctor<GridInterface, ResProp> func(this->ginterf_, this->res_prop_, s_orig, cap_press);
404 double cap_press_range = 1e2;
405 double mod_low = 1e100;
406 double mod_high = -1e100;
407 Opm::bracketZero(func, 0.0, cap_press_range, mod_low, mod_high);
408 const int max_iter = 40;
409 const double nonlinear_tolerance = 1e-12;
410 int iterations_used = -1;
411 typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
412 double mod_correct = RootFinder::solve(func, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
413 std::cout <<
"Moved capillary pressure solution by " << mod_correct <<
" after "
414 << iterations_used <<
" iterations." << std::endl;
415 s = func.lastSaturations();
419 template <
class Traits>
420 template <
class FlowSol>
422 std::pair<double, double>& oil_inout,
423 const FlowSol& flow_solution,
424 const std::vector<double>& saturations)
const
426 typedef typename GridInterface::CellIterator
CellIter;
427 typedef typename CellIter::FaceIterator
FaceIter;
429 double side1_flux = 0.0;
430 double side2_flux = 0.0;
431 double side1_flux_oil = 0.0;
432 double side2_flux_oil = 0.0;
433 std::map<int, double> frac_flow_by_bid;
434 int num_cells = this->ginterf_.numberOfCells();
435 std::vector<double> cell_inflows_w(num_cells, 0.0);
436 std::vector<double> cell_outflows_w(num_cells, 0.0);
441 for (
int pass = 0; pass < 2; ++pass) {
442 for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
443 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
445 double flux = flow_solution.outflux(f);
446 const SatBC& sc = this->bcond_.satCond(f);
447 if (flux < 0.0 && pass == 1) {
449 double frac_flow = 0.0;
450 if (sc.isPeriodic()) {
451 assert(sc.saturationDifference() == 0.0);
452 int partner_bid = this->bcond_.getPeriodicPartner(f->boundaryId());
453 std::map<int, double>::const_iterator it = frac_flow_by_bid.find(partner_bid);
454 if (it == frac_flow_by_bid.end()) {
455 OPM_THROW(std::runtime_error,
"Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
456 <<
" and partner bid = " << partner_bid);
458 frac_flow = it->second;
460 assert(sc.isDirichlet());
461 frac_flow = this->res_prop_.fractionalFlow(c->index(), sc.saturation());
463 cell_inflows_w[c->index()] += flux*frac_flow;
464 side1_flux += flux*frac_flow;
465 side1_flux_oil += flux*(1.0 - frac_flow);
466 }
else if (flux >= 0.0 && pass == 0) {
468 double frac_flow = this->res_prop_.fractionalFlow(c->index(), saturations[c->index()]);
469 if (sc.isPeriodic()) {
470 frac_flow_by_bid[f->boundaryId()] = frac_flow;
473 cell_outflows_w[c->index()] += flux*frac_flow;
474 side2_flux += flux*frac_flow;
475 side2_flux_oil += flux*(1.0 - frac_flow);
481 water_inout = std::make_pair(side1_flux, side2_flux);
482 oil_inout = std::make_pair(side1_flux_oil, side2_flux_oil);
491 #endif // OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
virtual void initImpl(const Opm::parameter::ParameterGroup ¶m)
Override from superclass.
Definition: SteadyStateUpscalerImplicit_impl.hpp:79
Definition: applier.hpp:18
SteadyStateUpscalerImplicit()
Default constructor.
Definition: SteadyStateUpscalerImplicit_impl.hpp:57
A base class for upscaling.
Definition: UpscalerBase.hpp:52
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:148
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:421
const std::vector< double > & lastSaturationState() const
Definition: SteadyStateUpscalerImplicit_impl.hpp:359
CellIter::FaceIterator FaceIter
Definition: UpscalerBase.hpp:128
Traits::template ResProp< Dimension >::Type ResProp
Definition: UpscalerBase.hpp:60
min[0]
Definition: elasticity_upscale_impl.hpp:144
void setToCapillaryLimit(double average_s, std::vector< double > &s) const
Definition: SteadyStateUpscalerImplicit_impl.hpp:395
void initSatLimits(std::vector< double > &s) const
Ensure saturations are not outside table.
Definition: SteadyStateUpscalerImplicit_impl.hpp:384
double lastSaturationUpscaled() const
Definition: SteadyStateUpscalerImplicit_impl.hpp:368
GridInterface::CellIterator CellIter
Definition: UpscalerBase.hpp:127
Super::permtensor_t permtensor_t
Definition: SteadyStateUpscalerImplicit.hpp:58