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));
176 const std::vector<double>& initial_saturation,
177 const double boundary_saturation,
178 const double pressure_drop,
179 const permtensor_t& upscaled_perm,
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 while ((!stationary) && (it_count < max_it_)) {
238 std::cout <<
"Running transport step " << it_count <<
" with stepsize "
239 << stepsize/Opm::unit::year <<
" years." << std::endl;
240 bool converged = transport_solver_.transportSolve(saturation, stepsize, gravity,
241 this->flow_solver_.getSolution(), injection);
248 this->flow_solver_.postProcessFluxes();
250 if (print_inoutflows_) {
251 std::pair<double, double> w_io, o_io;
252 computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
253 std::cout <<
"Pressure step " << it_count
254 <<
"\nWater flow [in] " << w_io.first
255 <<
" [out] " << w_io.second
256 <<
"\nOil flow [in] " << o_io.first
257 <<
" [out] " << o_io.second
262 writeVtkOutput(this->ginterf_,
264 this->flow_solver_.getSolution(),
266 std::string(
"output-steadystate")
267 +
'-' + std::to_string(count)
268 +
'-' + std::to_string(flow_direction)
269 +
'-' + std::to_string(it_count));
272 const char* fd =
"xyz";
273 std::string basename = std::string(
"ecldump-steadystate")
274 +
'-' + std::to_string(boundary_saturation)
275 +
'-' + std::to_string(fd[flow_direction])
276 +
'-' + std::to_string(pressure_drop);
277 waterSatToBothSat(saturation, ecl_sat);
278 getCellPressure(ecl_press, this->ginterf_, this->flow_solver_.getSolution());
279 data::Solution solution;
280 solution.insert(
"PRESSURE" , UnitSystem::measure::pressure , ecl_press , data::TargetType::RESTART_SOLUTION );
281 solution.insert(
"SWAT" , UnitSystem::measure::identity , destripe( ecl_sat, 2, 0) , data::TargetType::RESTART_SOLUTION );
282 ecl_time += stepsize;
283 boost::posix_time::ptime ecl_startdate( boost::gregorian::date(2012, 1, 1) );
284 boost::posix_time::ptime ecl_curdate = ecl_startdate + boost::posix_time::seconds(
int(ecl_time));
285 boost::posix_time::ptime epoch( boost::gregorian::date( 1970, 1, 1 ) );
286 auto ecl_posix_time = ( ecl_curdate - epoch ).total_seconds();
287 const auto* cgrid = grid_adapter_.c_grid();
288 Opm::writeECLData(cgrid->cartdims[ 0 ],
289 cgrid->cartdims[ 1 ],
290 cgrid->cartdims[ 2 ],
291 cgrid->number_of_cells,
293 ecl_time, ecl_posix_time,
297 double maxdiff = 0.0;
298 double euclidean_diff = 0.0;
299 for (
int i = 0; i < num_cells; ++i) {
300 const double sat_diff_cell = saturation[i] - saturation_old[i];
301 maxdiff = std::max(maxdiff, std::fabs(sat_diff_cell));
302 euclidean_diff += sat_diff_cell * sat_diff_cell * pore_vol[i];
304 euclidean_diff = std::sqrt(euclidean_diff / tot_pore_vol);
307 ds_year = maxdiff*Opm::unit::year/stepsize;
308 std::cout <<
"Maximum saturation change/year: " << ds_year << std::endl;
309 if (maxdiff < dt_sat_tol_) {
310 stepsize=std::min(max_stepsize_,2*stepsize);
314 ds_year = euclidean_diff*Opm::unit::year/stepsize;
315 std::cout <<
"Euclidean saturation change/year: " << ds_year << std::endl;
316 if (euclidean_diff < dt_sat_tol_) {
317 stepsize=std::min(max_stepsize_,2*stepsize);
320 if (ds_year < sat_change_year_) {
324 std::cerr <<
"Cutting time step\n";
325 stepsize=stepsize/2.0;
329 saturation_old = saturation;
331 success = stationary;
335 typedef typename Super::ResProp::Mobility Mob;
339 for (
int c = 0; c < num_cells; ++c) {
340 this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
341 m1max = maxMobility(m1max, m.mob);
342 this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
343 m2max = maxMobility(m2max, m.mob);
346 const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
347 const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
348 const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
349 const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
350 const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
351 const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
353 std::vector<Mob> mob1(num_cells);
354 std::vector<Mob> mob2(num_cells);
355 for (
int c = 0; c < num_cells; ++c) {
356 this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
357 thresholdMobility(mob1[c].mob, mob1_threshold);
358 this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
359 thresholdMobility(mob2[c].mob, mob2_threshold);
364 permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
366 permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
369 last_saturation_state_.swap(saturation);
374 permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
375 permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
379 permtensor_t k_rw(lambda_w);
380 k_rw *= this->res_prop_.viscosityFirstPhase();
381 permtensor_t k_ro(lambda_o);
382 k_ro *= this->res_prop_.viscositySecondPhase();
383 return std::make_pair(k_rw, k_ro);