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));
133 const std::vector<double>& initial_saturation,
134 const double boundary_saturation,
135 const double pressure_drop,
136 const permtensor_t& upscaled_perm)
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
204 writeVtkOutput(this->ginterf_,
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);
258 thresholdMobility(mob1[c].mob, mob1_threshold);
259 this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
260 thresholdMobility(mob2[c].mob, mob2_threshold);
265 permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
267 permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
270 last_saturation_state_.swap(saturation);
275 permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
276 permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
280 permtensor_t k_rw(lambda_w);
281 k_rw *= this->res_prop_.viscosityFirstPhase();
282 permtensor_t k_ro(lambda_o);
283 k_ro *= this->res_prop_.viscositySecondPhase();
284 return std::make_pair(k_rw, k_ro);