36 #ifndef OPM_STEADYSTATEUPSCALER_IMPL_HEADER
37 #define OPM_STEADYSTATEUPSCALER_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/core/utility/Units.hpp>
53 template <
class Traits>
58 print_inoutflows_(false),
59 simulation_steps_(10),
61 relperm_threshold_(1.0e-8),
62 maximum_mobility_contrast_(1.0e9),
63 sat_change_threshold_(0.0)
70 template <
class Traits>
73 Super::initImpl(param);
74 use_gravity_ = param.getDefault(
"use_gravity", use_gravity_);
75 output_vtk_ = param.getDefault(
"output_vtk", output_vtk_);
76 print_inoutflows_ = param.getDefault(
"print_inoutflows", print_inoutflows_);
77 simulation_steps_ = param.getDefault(
"simulation_steps", simulation_steps_);
78 stepsize_ = Opm::unit::convert::from(param.getDefault(
"stepsize", stepsize_),
80 relperm_threshold_ = param.getDefault(
"relperm_threshold", relperm_threshold_);
81 maximum_mobility_contrast_ = param.getDefault(
"maximum_mobility_contrast", maximum_mobility_contrast_);
82 sat_change_threshold_ = param.getDefault(
"sat_change_threshold", sat_change_threshold_);
84 transport_solver_.init(param);
86 double v1_default = this->res_prop_.viscosityFirstPhase();
87 double v2_default = this->res_prop_.viscositySecondPhase();
88 this->res_prop_.setViscosities(param.getDefault(
"viscosity1", v1_default), param.getDefault(
"viscosity2", v2_default));
89 double d1_default = this->res_prop_.densityFirstPhase();
90 double d2_default = this->res_prop_.densitySecondPhase();
91 this->res_prop_.setDensities(param.getDefault(
"density1", d1_default), param.getDefault(
"density2", d2_default));
98 double maxMobility(
double m1,
double m2)
100 return std::max(m1, m2);
103 template <
class SomeMatrixType>
104 double maxMobility(
double m1, SomeMatrixType& m2)
107 for (
int i = 0; i <
std::min(m2.numRows(), m2.numCols()); ++i) {
108 m = std::max(m, m2(i,i));
112 void thresholdMobility(
double& m,
double threshold)
114 m = std::max(m, threshold);
117 template <
class SomeMatrixType>
118 void thresholdMobility(SomeMatrixType& m,
double threshold)
120 for (
int i = 0; i <
std::min(m.numRows(), m.numCols()); ++i) {
121 m(i, i) = std::max(m(i, i), threshold);
129 template <
class Traits>
130 inline std::pair<typename SteadyStateUpscaler<Traits>::permtensor_t,
131 typename SteadyStateUpscaler<Traits>::permtensor_t>
134 const std::vector<double>& initial_saturation,
135 const double boundary_saturation,
136 const double pressure_drop,
139 static int count = 0;
141 int num_cells = this->ginterf_.numberOfCells();
143 std::vector<double> src(num_cells, 0.0);
144 Opm::SparseVector<double> injection(num_cells);
146 Dune::FieldVector<double, 3> gravity(0.0);
148 gravity[2] = Opm::unit::gravity;
150 if (gravity.two_norm() > 0.0) {
151 OPM_MESSAGE(
"Warning: Gravity is experimental for flow solver.");
155 std::vector<double> saturation = initial_saturation;
158 setupUpscalingConditions(this->ginterf_, this->bctype_, flow_direction,
159 pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
162 if (flow_direction == 0) {
163 this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
165 transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
168 this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
169 this->residual_tolerance_, this->linsolver_verbosity_,
170 this->linsolver_type_,
false,
171 this->linsolver_maxit_, this->linsolver_prolongate_factor_,
172 this->linsolver_smooth_steps_);
173 double max_mod = this->flow_solver_.postProcessFluxes();
174 std::cout <<
"Max mod = " << max_mod << std::endl;
177 std::vector<double> saturation_old = saturation;
178 for (
int iter = 0; iter < simulation_steps_; ++iter) {
180 transport_solver_.transportSolve(saturation, stepsize_, gravity, this->flow_solver_.getSolution(), injection);
183 this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
184 this->residual_tolerance_, this->linsolver_verbosity_,
185 this->linsolver_type_,
false,
186 this->linsolver_maxit_, this->linsolver_prolongate_factor_,
187 this->linsolver_smooth_steps_);
188 max_mod = this->flow_solver_.postProcessFluxes();
189 std::cout <<
"Max mod = " << max_mod << std::endl;
192 if (print_inoutflows_) {
193 std::pair<double, double> w_io, o_io;
194 computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
195 std::cout <<
"Pressure step " << iter
196 <<
"\nWater flow [in] " << w_io.first
197 <<
" [out] " << w_io.second
198 <<
"\nOil flow [in] " << o_io.first
199 <<
" [out] " << o_io.second
205 writeVtkOutput(this->ginterf_,
207 this->flow_solver_.getSolution(),
209 std::string(
"output-steadystate")
210 +
'-' + boost::lexical_cast<std::string>(count)
211 +
'-' + boost::lexical_cast<std::string>(flow_direction)
212 +
'-' + boost::lexical_cast<std::string>(iter));
216 int num_cells = saturation.size();
217 double maxdiff = 0.0;
218 for (
int i = 0; i < num_cells; ++i) {
219 maxdiff = std::max(maxdiff, std::fabs(saturation[i] - saturation_old[i]));
222 std::cout <<
"Maximum saturation change: " << maxdiff << std::endl;
224 if (maxdiff < sat_change_threshold_) {
226 std::cout <<
"Maximum saturation change is under steady state threshold." << std::endl;
232 saturation_old = saturation;
237 typedef typename Super::ResProp::Mobility Mob;
241 for (
int c = 0; c < num_cells; ++c) {
242 this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
243 m1max = maxMobility(m1max, m.mob);
244 this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
245 m2max = maxMobility(m2max, m.mob);
248 const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
249 const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
250 const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
251 const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
252 const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
253 const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
255 std::vector<Mob> mob1(num_cells);
256 std::vector<Mob> mob2(num_cells);
257 for (
int c = 0; c < num_cells; ++c) {
258 this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
259 thresholdMobility(mob1[c].mob, mob1_threshold);
260 this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
261 thresholdMobility(mob2[c].mob, mob2_threshold);
265 ReservoirPropertyFixedMobility<Mob> fluid_first(mob1);
266 permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
267 ReservoirPropertyFixedMobility<Mob> fluid_second(mob2);
268 permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
271 last_saturation_state_.swap(saturation);
276 permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
277 permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
282 k_rw *= this->res_prop_.viscosityFirstPhase();
284 k_ro *= this->res_prop_.viscositySecondPhase();
285 return std::make_pair(k_rw, k_ro);
291 template <
class Traits>
292 inline const std::vector<double>&
295 return last_saturation_state_;
301 template <
class Traits>
304 typedef typename GridInterface::CellIterator
CellIter;
305 double pore_vol = 0.0;
306 double sat_vol = 0.0;
307 for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
308 double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
309 pore_vol += cell_pore_vol;
310 sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
313 return sat_vol/pore_vol;
318 std::map< int, double >& frac_flow_by_bid,
320 const ResProp& res_prop,
321 GridInterface::CellIterator c,
322 CellIterator::FaceIterator f ) {
324 const SatBC& sc = bcond.satCond( f );
326 if( sc.isPeriodic() ) {
327 assert( sc.isDirichlet() );
328 return res_prop.fractionalFlow( c->index(), sc.saturation() );
331 assert( sc.saturationDifference() == 0.0 );
332 const int partner_bid = bcond.getPeriodicPartner( f->boundaryId() );
333 const auto it = frac_flow_by_bid.find( partner_bid );
337 if( it == frac_flow_by_bid.end() ) {
338 OPM_THROW(std::runtime_error,
"Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
339 <<
" and partner bid = " << partner_bid);
345 template <
class Traits>
346 template <
class FlowSol>
348 std::pair<double, double>& oil_inout,
349 const FlowSol& flow_solution,
350 const std::vector<double>& saturations)
const
357 sideflux& operator+=(
const sideflux& other ) {
358 this->water += other.water;
359 this->oil += other.oil;
364 std::map< int, double > frac_flow_by_bid;
367 for(
auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
368 const double frac_flow = this->res_prop_.fractional_flow( c->index(), saturations[ c->index() ] );
370 for(
auto f = c->facebegin(); f != c->faceend(); ++f ) {
371 if( !f->boundary() )
continue;
373 const double flux = flow_solution.outflux( f );
374 if( flux < 0.0 )
continue;
376 const SatBC& sc = this->bcond_.satCond( f );
377 if( sc.isPeriodic() ) frac_flow_by_bid[ f->boundaryId() ] = frac_flow;
379 outflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
385 for(
auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
386 for(
auto f = c->facebegin(); f != c->faceend(); ++f ) {
387 if( !f->boundary() )
continue;
389 const double flux = flow_solution.outflux( f );
390 if( flux >= 0.0 )
continue;
392 const double frac_flow =
calc_frac_flow( frac_flow_by_bid, this->bcond_, this->res_prop_, c, f );
394 inflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
398 water_inout = std::make_pair( in.water, out.water );
399 oil_inout = std::make_pair( in.oil, out.oil );
405 #endif // OPM_STEADYSTATEUPSCALER_IMPL_HEADER
SteadyStateUpscaler()
Default constructor.
Definition: SteadyStateUpscaler_impl.hpp:54
Definition: applier.hpp:18
const std::vector< double > & lastSaturationState() const
Definition: SteadyStateUpscaler_impl.hpp:293
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)
Definition: SteadyStateUpscaler_impl.hpp:133
A base class for upscaling.
Definition: UpscalerBase.hpp:52
virtual void initImpl(const Opm::parameter::ParameterGroup ¶m)
Override from superclass.
Definition: SteadyStateUpscaler_impl.hpp:71
min[0]
Definition: elasticity_upscale_impl.hpp:144
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: SteadyStateUpscaler_impl.hpp:347
double lastSaturationUpscaled() const
Definition: SteadyStateUpscaler_impl.hpp:302
static double calc_frac_flow(std::map< int, double > &frac_flow_by_bid, const BCs &bcond, const ResProp &res_prop, GridInterface::CellIterator c, CellIterator::FaceIterator f)
Definition: SteadyStateUpscaler_impl.hpp:317
GridInterface::CellIterator CellIter
Definition: UpscalerBase.hpp:127
Super::permtensor_t permtensor_t
Definition: SteadyStateUpscaler.hpp:57