37#ifndef OPENRS_EULERUPSTREAM_IMPL_HEADER
38#define OPENRS_EULERUPSTREAM_IMPL_HEADER
42#include <opm/common/ErrorMacros.hpp>
43#include <opm/core/utility/Average.hpp>
44#include <opm/core/utility/Units.hpp>
45#include <dune/grid/common/Volumes.hpp>
47#include <opm/core/utility/StopWatch.hpp>
59 template <
class GI,
class RP,
class BC>
61 : method_viscous_(true),
62 method_gravity_(true),
63 method_capillary_(true),
64 use_cfl_viscous_(true),
65 use_cfl_gravity_(true),
66 use_cfl_capillary_(true),
68 minimum_small_steps_(1),
69 maximum_small_steps_(10000),
75 template <
class GI,
class RP,
class BC>
77 : method_viscous_(true),
78 method_gravity_(true),
79 method_capillary_(true),
80 use_cfl_viscous_(true),
81 use_cfl_gravity_(true),
82 use_cfl_capillary_(true),
84 minimum_small_steps_(1),
85 maximum_small_steps_(10000),
94 template <
class GI,
class RP,
class BC>
97 courant_number_ = param.getDefault(
"courant_number", courant_number_);
98 method_viscous_ = param.getDefault(
"method_viscous", method_viscous_);
99 method_gravity_ = param.getDefault(
"method_gravity", method_gravity_);
100 method_capillary_ = param.getDefault(
"method_capillary", method_capillary_);
101 use_cfl_viscous_ = param.getDefault(
"use_cfl_viscous", use_cfl_viscous_);
102 use_cfl_gravity_ = param.getDefault(
"use_cfl_gravity", use_cfl_gravity_);
103 use_cfl_capillary_ = param.getDefault(
"use_cfl_capillary", use_cfl_capillary_);
104 minimum_small_steps_ = param.getDefault(
"minimum_small_steps", minimum_small_steps_);
105 maximum_small_steps_ = param.getDefault(
"maximum_small_steps", maximum_small_steps_);
106 check_sat_ = param.getDefault(
"check_sat", check_sat_);
107 clamp_sat_ = param.getDefault(
"clamp_sat", clamp_sat_);
110 template <
class GI,
class RP,
class BC>
112 const GI& g,
const RP& r,
const BC& b)
119 template <
class GI,
class RP,
class BC>
122 residual_computer_.initObj(g, r, b);
123 porevol_.resize(g.numberOfCells());
124 for (
CIt c = g.cellbegin(); c != g.cellend(); ++c) {
125 porevol_[c->index()] = c->volume()*r.porosity(c->index());
131 template <
class GI,
class RP,
class BC>
136 cout <<
"Displaying some members of EulerUpstream" << endl;
138 cout <<
"courant_number = " << courant_number_ << endl;
143 template <
class GI,
class RP,
class BC>
146 courant_number_ = cn;
151 template <
class GI,
class RP,
class BC>
152 template <
class PressureSolution>
155 const typename GI::Vector& gravity,
156 const PressureSolution& pressure_sol,
157 const Opm::SparseVector<double>& injection_rates)
const
160 double cfl_dt = computeCflTime(saturation, time, gravity, pressure_sol);
163 int nr_transport_steps;
165 nr_transport_steps = minimum_small_steps_;
167 double steps = std::min<double>(std::ceil(time/cfl_dt), std::numeric_limits<int>::max());
168 nr_transport_steps = int(steps);
169 nr_transport_steps = std::max(nr_transport_steps, minimum_small_steps_);
170 nr_transport_steps = std::min(nr_transport_steps, maximum_small_steps_);
172 double dt_transport = time/nr_transport_steps;
181 std::vector<double> saturation_initial(saturation);
182 bool finished =
false;
184 const int max_repeats = 10;
185 Opm::time::StopWatch clock;
190 std::cout <<
"Doing " << nr_transport_steps
191 <<
" steps for saturation equation with stepsize "
192 << dt_transport <<
" in seconds." << std::endl;
194 for (
int q = 0; q < nr_transport_steps; ++q) {
195 smallTimeStep(saturation,
205 if (repeats > max_repeats) {
208 OPM_MESSAGE(
"Warning: Transport failed, retrying with more steps.");
209 nr_transport_steps *= 2;
210 dt_transport = time/nr_transport_steps;
211 saturation = saturation_initial;
216 std::cout <<
"Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
263 template <
class GI,
class RP,
class BC>
264 template <
class PressureSolution>
271 const typename GI::Vector& gravity,
272 const PressureSolution& pressure_sol)
const
275 double cfl_dt_v = 1e99;
276 double cfl_dt_g = 1e99;
277 double cfl_dt_c = 1e99;
280 if (method_viscous_ && use_cfl_viscous_) {
282 residual_computer_.reservoirProperties(),
285 std::cout <<
"CFL dt for velocity is "
286 << cfl_dt_v <<
" seconds ("
287 << Opm::unit::convert::to(cfl_dt_v, Opm::unit::day)
288 <<
" days)." << std::endl;
293 if (method_gravity_ && use_cfl_gravity_) {
295 residual_computer_.reservoirProperties(),
298 std::cout <<
"CFL dt for gravity is "
299 << cfl_dt_g <<
" seconds ("
300 << Opm::unit::convert::to(cfl_dt_g, Opm::unit::day)
301 <<
" days)." << std::endl;
306 if (method_capillary_ && use_cfl_capillary_) {
308 residual_computer_.reservoirProperties());
310 std::cout <<
"CFL dt for capillary term is "
311 << cfl_dt_c <<
" seconds ("
312 << Opm::unit::convert::to(cfl_dt_c, Opm::unit::day)
313 <<
" days)." << std::endl;
317 double cfl_dt = std::min(std::min(cfl_dt_v, cfl_dt_g), cfl_dt_c);
318 cfl_dt *= courant_number_;
320 std::cout <<
"Total impes time is "
321 << time <<
" seconds ("
322 << Opm::unit::convert::to(time, Opm::unit::day)
323 <<
" days)." << std::endl;
325 std::cout <<
"Final modified CFL dt is "
326 << cfl_dt <<
" seconds ("
327 << Opm::unit::convert::to(cfl_dt, Opm::unit::day)
328 <<
" days)." << std::endl;
336 template <
class GI,
class RP,
class BC>
339 int num_cells = s.size();
340 for (
int cell = 0; cell < num_cells; ++cell) {
341 if (s[cell] > 1.0 || s[cell] < 0.0) {
343 s[cell] = std::max(std::min(s[cell], 1.0), 0.0);
344 }
else if (s[cell] > 1.001 || s[cell] < -0.001) {
345 OPM_THROW(std::runtime_error,
"Saturation out of range in EulerUpstream: Cell " << cell <<
" sat " << s[cell]);
355 template <
class GI,
class RP,
class BC>
356 template <
class PressureSolution>
359 const typename GI::Vector& gravity,
360 const PressureSolution& pressure_sol,
361 const Opm::SparseVector<double>& injection_rates)
const
363 if (method_capillary_) {
364 residual_computer_.computeCapPressures(saturation);
366 residual_computer_.computeResidual(saturation, gravity, pressure_sol, injection_rates,
367 method_viscous_, method_gravity_, method_capillary_,
371 int num_cells = saturation.size();
372 for (
int i = 0; i < num_cells; ++i) {
373 const double sat_change = dt*residual_[i]/porevol_[i];
374 saturation[i] += sat_change;
382 if (check_sat_ || clamp_sat_) {
383 checkAndPossiblyClampSat(saturation);
void transportSolve(std::vector< double > &saturation, const double time, const typename GridInterface::Vector &gravity, const PressureSolution &pressure_sol, const Opm::SparseVector< double > &injection_rates) const
Solve transport equation, evolving.
void setCourantNumber(double cn)
Set the Courant number. That is dt = dt_cfl*courant_number. For this explicit method it should be < 1...
Definition: EulerUpstream_impl.hpp:144
double computeCflTime(const std::vector< double > &saturation, const double time, const typename GridInterface::Vector &gravity, const PressureSolution &pressure_sol) const
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary)
Definition: EulerUpstream_impl.hpp:120
EulerUpstream()
Definition: EulerUpstream_impl.hpp:60
void init(const Opm::parameter::ParameterGroup ¶m)
Definition: EulerUpstream_impl.hpp:95
GridInterface::CellIterator CIt
Definition: EulerUpstream.hpp:105
void display()
Definition: EulerUpstream_impl.hpp:132
void smallTimeStep(std::vector< double > &saturation, const double time, const typename GridInterface::Vector &gravity, const PressureSolution &pressure_sol, const Opm::SparseVector< double > &injection_rates) const
EulerUpstreamResidual< GridInterface, ReservoirProperties, BoundaryConditions > residual_computer_
Definition: EulerUpstream.hpp:128
void checkAndPossiblyClampSat(std::vector< double > &s) const
Definition: EulerUpstream_impl.hpp:337
double findCFLtimeVelocity(const Grid &grid, const ReservoirProperties &resprop, const PressureSolution &pressure_sol)
Definition: CflCalculator.hpp:55
double findCFLtimeGravity(const Grid &grid, const ReservoirProperties &resprop, const typename Grid::Vector &gravity)
Definition: CflCalculator.hpp:91
double findCFLtimeCapillary(const Grid &grid, const ReservoirProperties &resprop)
Definition: CflCalculator.hpp:143
Definition: BlackoilFluid.hpp:32