37#ifndef OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
38#define OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
40#include <opm/common/ErrorMacros.hpp>
41#include <opm/core/utility/Average.hpp>
42#include <opm/core/utility/Units.hpp>
43#include <dune/grid/common/Volumes.hpp>
44#include <opm/core/utility/StopWatch.hpp>
46#include <opm/core/pressure/tpfa/trans_tpfa.h>
59 template <
class GI,
class RP,
class BC>
65 template <
class GI,
class RP,
class BC>
73 template <
class GI,
class RP,
class BC>
76 check_sat_ = param.getDefault(
"check_sat", check_sat_);
77 clamp_sat_ = param.getDefault(
"clamp_sat", clamp_sat_);
79 ctrl_.max_it = param.getDefault(
"transport_nr_max_it", 10);
80 max_repeats_ = param.getDefault(
"transport_max_rep", 10);
81 ctrl_.atol = param.getDefault(
"transport_atol", 1.0e-6);
82 ctrl_.rtol = param.getDefault(
"transport_rtol", 5.0e-7);
83 ctrl_.max_it_ls = param.getDefault(
"transport_max_it_ls", 20);
84 ctrl_.dxtol = param.getDefault(
"transport_dxtol", 1e-6);
85 ctrl_.verbosity = param.getDefault(
"transport_verbosity", 0);
88 template <
class GI,
class RP,
class BC>
90 const GI& g,
const RP& r,
const BC& b)
97 template <
class GI,
class RP,
class BC>
102 mygrid_.init(g.grid());
103 porevol_.resize(mygrid_.numCells());
104 for (
int i = 0; i < mygrid_.numCells(); ++i){
105 porevol_[i]= mygrid_.cellVolume(i)*r.porosity(i);
108 int num_cells = mygrid_.numCells();
109 int ngconn = mygrid_.c_grid()->cell_facepos[num_cells];
111 htrans_.resize(ngconn);
112 const double* perm = &(r.permeability(0)(0,0));
113 tpfa_htrans_compute(mygrid_.c_grid(), perm, &htrans_[0]);
118 typedef typename GI::CellIterator CIt;
119 typedef typename CIt::FaceIterator FIt;
120 std::vector<FIt> bid_to_face;
122 for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
123 for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
124 int bid = f->boundaryId();
125 maxbid = std::max(maxbid, bid);
129 bid_to_face.resize(maxbid + 1);
130 std::vector<int> egf_cf(mygrid_.numFaces());
132 for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
134 for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
135 if (f->boundary() && b.satCond(*f).isPeriodic()) {
136 bid_to_face[f->boundaryId()] = f;
139 int cf=mygrid_.cellFace(cix,loc_fix);
147 const UnstructuredGrid& c_grid=*mygrid_.c_grid();
151 periodic_cells_.resize(0);
152 periodic_faces_.resize(0);
153 periodic_hfaces_.resize(0);
154 periodic_nbfaces_.resize(0);
156 direclet_cells_.resize(0);
157 direclet_sat_.resize(0);
158 direclet_sat_.resize(0);
159 direclet_hfaces_.resize(0);
161 assert(periodic_cells_.size()==0);
162 for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
163 int cell0 = c->index();
164 for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
170 if (b.satCond(*f).isPeriodic()) {
171 nbface = bid_to_face[b.getPeriodicPartner(f->boundaryId())];
173 int cell1 = nbface->cellIndex();
174 assert(cell0 != cell1);
176 int f_ind=f->index();
178 int fn_ind=nbface->index();
181 fn_ind=egf_cf[fn_ind];
182 assert((c_grid.face_cells[2*f_ind]==-1) || (c_grid.face_cells[2*f_ind+1]==-1));
183 assert((c_grid.face_cells[2*fn_ind]==-1) || (c_grid.face_cells[2*fn_ind+1]==-1));
184 assert((c_grid.face_cells[2*f_ind]==cell0) || (c_grid.face_cells[2*f_ind+1]==cell0));
185 assert((c_grid.face_cells[2*fn_ind]==cell1) || (c_grid.face_cells[2*fn_ind+1]==cell1));
186 periodic_cells_.push_back(cell0);
187 periodic_cells_.push_back(cell1);
188 periodic_faces_.push_back(f_ind);
189 periodic_hfaces_.push_back(hf_ind);
190 periodic_nbfaces_.push_back(fn_ind);
191 }
else if (!( b.flowCond(*f).isNeumann() && b.flowCond(*f).outflux() == 0.0)) {
193 direclet_cells_.push_back(cell0);
194 direclet_sat_.push_back(b.satCond(*f).saturation());
195 direclet_sat_.push_back(1-b.satCond(*f).saturation());
196 direclet_hfaces_.push_back(hf_ind);
203 mygrid_.makeQPeriodic(periodic_hfaces_,periodic_cells_);
206 int num_b=direclet_cells_.size();
207 for(
int i=0; i <num_b; ++i){
208 std::array<double,2> sat = {{direclet_sat_[2*i] ,direclet_sat_[2*i+1] }};
209 std::array<double,2> mob;
210 std::array<double,2*2> dmob;
211 myfluid.
mobility(direclet_cells_[i], sat, mob, dmob);
212 double fl = mob[0]/(mob[0]+mob[1]);
213 direclet_sat_[2*i] = fl;
214 direclet_sat_[2*i+1] = 1-fl;
218 template <
class GI,
class RP,
class BC>
223 cout <<
"Displaying some members of EulerUpstreamImplicit" << endl;
227 template <
class GI,
class RP,
class BC>
228 template <
class PressureSolution>
231 const typename GI::Vector& gravity,
232 const PressureSolution& pressure_sol,
233 const Opm::SparseVector<double>& injection_rates)
const
238 std::vector<double>& sat = state.
saturation();
239 for (
int i=0; i < mygrid_.numCells(); ++i){
240 sat[2*i] = saturation[i];
241 sat[2*i+1] = 1-saturation[i];
246 const UnstructuredGrid* cgrid = mygrid_.c_grid();
247 int numhf = cgrid->cell_facepos[cgrid->number_of_cells];
249 std::vector<double> faceflux(numhf);
251 for (
int c = 0, i = 0; c < cgrid->number_of_cells; ++c){
252 for (; i < cgrid->cell_facepos[c + 1]; ++i) {
253 int f= cgrid->cell_faces[i];
254 double outflux = pressure_sol.outflux(i);
255 double sgn = 2.0*(cgrid->face_cells[2*f + 0] == c) - 1;
256 faceflux[f] = sgn * outflux;
259 int num_db=direclet_hfaces_.size();
260 std::vector<double> sflux(num_db);
261 for (
int i=0; i < num_db;++i){
262 sflux[i]=-pressure_sol.outflux(direclet_hfaces_[i]);
266 double dt_transport = time;
267 int nr_transport_steps = 1;
268 Opm::time::StopWatch clock;
270 bool finished =
false;
275 const UnstructuredGrid& c_grid=*mygrid_.c_grid();
277 model.makefhfQPeriodic(periodic_faces_,periodic_hfaces_, periodic_nbfaces_);
278 model.initGravityTrans(*mygrid_.c_grid(),htrans_);
281 Opm::ImplicitTransportDetails::NRReport rpt_;
286 tsrc.
nsrc =direclet_cells_.size();
288 tsrc.
cell = direclet_cells_;
292 for (
int q = 0; q < nr_transport_steps; ++q) {
293 tsolver.solve(*mygrid_.c_grid(), &tsrc, dt_transport, ctrl_, state, linsolve_, rpt_);
301 if(repeats >max_repeats_){
304 OPM_MESSAGE(
"Warning: Transport failed, retrying with more steps.");
305 nr_transport_steps *= 2;
306 dt_transport = time/nr_transport_steps;
307 if (ctrl_.verbosity){
308 std::cout <<
"Warning: Transport failed, retrying with more steps. dt = "
309 << dt_transport/Opm::unit::year <<
" year.\n";
312 std::vector<double>& sat = state.
saturation();
313 for (
int i=0; i < mygrid_.numCells(); ++i){
314 sat[2*i] = saturation[i];
315 sat[2*i+1] = 1-saturation[i];
322 std::cout <<
"EulerUpstreamImplicite used " << repeats
323 <<
" repeats and " << nr_transport_steps <<
" steps"<< std::endl;
325 std::cout <<
"Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
328 std::vector<double>& sat = state.
saturation();
329 for (
int i=0; i < mygrid_.numCells(); ++i){
330 saturation[i] = sat[2*i];
334 std::cerr <<
"EulerUpstreamImplicit did not converge" << std::endl;
341 template <
class GI,
class RP,
class BC>
344 int num_cells = s.size();
345 for (
int cell = 0; cell < num_cells; ++cell) {
346 if (s[cell] > 1.0 || s[cell] < 0.0) {
348 s[cell] = std::max(std::min(s[cell], 1.0), 0.0);
349 }
else if (s[cell] > 1.001 || s[cell] < -0.001) {
350 OPM_THROW(std::runtime_error,
"Saturation out of range in EulerUpstreamImplicit: Cell " << cell <<
" sat " << s[cell]);
EulerUpstreamImplicit()
Definition: EulerUpstreamImplicit_impl.hpp:60
void checkAndPossiblyClampSat(std::vector< double > &s) const
Definition: EulerUpstreamImplicit_impl.hpp:342
Opm::SinglePointUpwindTwoPhase< TwophaseFluid > TransportModel
Definition: EulerUpstreamImplicit.hpp:123
Opm::ImplicitTransport< TransportModel, JacSys, Opm::MaxNormDune, Opm::ImplicitTransportDefault::VectorNegater, Opm::ImplicitTransportDefault::VectorZero, Opm::ImplicitTransportDefault::MatrixZero, Opm::ImplicitTransportDefault::VectorAssign > TransportSolver
Definition: EulerUpstreamImplicit.hpp:146
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary)
Definition: EulerUpstreamImplicit_impl.hpp:98
bool 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 display()
Definition: EulerUpstreamImplicit_impl.hpp:219
void init(const Opm::parameter::ParameterGroup ¶m)
Definition: EulerUpstreamImplicit_impl.hpp:74
Definition: ImplicitTransportDefs.hpp:174
Definition: ImplicitTransportDefs.hpp:57
::std::vector< double > & faceflux()
Definition: ImplicitTransportDefs.hpp:68
::std::vector< double > & saturation()
Definition: ImplicitTransportDefs.hpp:69
Definition: ImplicitTransportDefs.hpp:237
::std::vector< double > saturation
Definition: ImplicitTransportDefs.hpp:246
::std::vector< int > cell
Definition: ImplicitTransportDefs.hpp:243
int nsrc
Definition: ImplicitTransportDefs.hpp:241
::std::vector< double > flux
Definition: ImplicitTransportDefs.hpp:245
Definition: ImplicitTransportDefs.hpp:120
void mobility(int c, const Sat &s, Mob &mob, DMob &dmob) const
Definition: ImplicitTransportDefs.hpp:139
Definition: BlackoilFluid.hpp:32