37#ifndef OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
38#define OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
40#include <opm/common/ErrorMacros.hpp>
42#include <opm/input/eclipse/Units/Units.hpp>
43#include <opm/grid/common/Volumes.hpp>
44#include <opm/grid/utility/StopWatch.hpp>
46#include <opm/grid/transmissibility/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();
150 periodic_cells_.resize(0);
151 periodic_faces_.resize(0);
152 periodic_hfaces_.resize(0);
153 periodic_nbfaces_.resize(0);
155 direclet_cells_.resize(0);
156 direclet_sat_.resize(0);
157 direclet_sat_.resize(0);
158 direclet_hfaces_.resize(0);
160 assert(periodic_cells_.size()==0);
161 for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
162 int cell0 = c->index();
163 for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
168 if (b.satCond(*f).isPeriodic()) {
169 nbface = bid_to_face[b.getPeriodicPartner(f->boundaryId())];
171 int cell1 = nbface->cellIndex();
172 assert(cell0 != cell1);
174 int f_ind=f->index();
176 int fn_ind=nbface->index();
179 fn_ind=egf_cf[fn_ind];
180 assert((c_grid.face_cells[2*f_ind]==-1) || (c_grid.face_cells[2*f_ind+1]==-1));
181 assert((c_grid.face_cells[2*fn_ind]==-1) || (c_grid.face_cells[2*fn_ind+1]==-1));
182 assert((c_grid.face_cells[2*f_ind]==cell0) || (c_grid.face_cells[2*f_ind+1]==cell0));
183 assert((c_grid.face_cells[2*fn_ind]==cell1) || (c_grid.face_cells[2*fn_ind+1]==cell1));
184 periodic_cells_.push_back(cell0);
185 periodic_cells_.push_back(cell1);
186 periodic_faces_.push_back(f_ind);
187 periodic_hfaces_.push_back(hf_ind);
188 periodic_nbfaces_.push_back(fn_ind);
189 }
else if (!( b.flowCond(*f).isNeumann() && b.flowCond(*f).outflux() == 0.0)) {
191 direclet_cells_.push_back(cell0);
192 direclet_sat_.push_back(b.satCond(*f).saturation());
193 direclet_sat_.push_back(1-b.satCond(*f).saturation());
194 direclet_hfaces_.push_back(hf_ind);
201 mygrid_.makeQPeriodic(periodic_hfaces_,periodic_cells_);
204 int num_b=direclet_cells_.size();
205 for(
int i=0; i <num_b; ++i){
206 std::array<double,2> sat = {{direclet_sat_[2*i] ,direclet_sat_[2*i+1] }};
207 std::array<double,2> mob;
208 std::array<double,2*2> dmob;
209 myfluid.
mobility(direclet_cells_[i], sat, mob, dmob);
210 double fl = mob[0]/(mob[0]+mob[1]);
211 direclet_sat_[2*i] = fl;
212 direclet_sat_[2*i+1] = 1-fl;
216 template <
class GI,
class RP,
class BC>
221 cout <<
"Displaying some members of EulerUpstreamImplicit" << endl;
225 template <
class GI,
class RP,
class BC>
226 template <
class PressureSolution>
230 const PressureSolution& pressure_sol,
231 const Opm::SparseVector<double>& )
const
236 std::vector<double>& sat = state.
saturation();
237 for (
int i=0; i < mygrid_.numCells(); ++i){
238 sat[2*i] = saturation[i];
239 sat[2*i+1] = 1-saturation[i];
244 const UnstructuredGrid* cgrid = mygrid_.c_grid();
245 int numhf = cgrid->cell_facepos[cgrid->number_of_cells];
247 std::vector<double> faceflux(numhf);
249 for (
int c = 0, i = 0; c < cgrid->number_of_cells; ++c){
250 for (; i < cgrid->cell_facepos[c + 1]; ++i) {
251 int f= cgrid->cell_faces[i];
252 double outflux = pressure_sol.outflux(i);
253 double sgn = 2.0*(cgrid->face_cells[2*f + 0] == c) - 1;
254 faceflux[f] = sgn * outflux;
257 int num_db=direclet_hfaces_.size();
258 std::vector<double> sflux(num_db);
259 for (
int i=0; i < num_db;++i){
260 sflux[i]=-pressure_sol.outflux(direclet_hfaces_[i]);
264 double dt_transport = time;
265 int nr_transport_steps = 1;
266 Opm::time::StopWatch clock;
268 bool finished =
false;
273 const UnstructuredGrid& c_grid=*mygrid_.c_grid();
275 model.
makefhfQPeriodic(periodic_faces_,periodic_hfaces_, periodic_nbfaces_);
284 tsrc.
nsrc =direclet_cells_.size();
286 tsrc.
cell = direclet_cells_;
290 for (
int q = 0; q < nr_transport_steps; ++q) {
291 tsolver.
solve(*mygrid_.c_grid(), &tsrc, dt_transport, ctrl_, state, linsolve_, rpt_);
299 if(repeats >max_repeats_){
302 OPM_MESSAGE(
"Warning: Transport failed, retrying with more steps.");
303 nr_transport_steps *= 2;
304 dt_transport = time/nr_transport_steps;
305 if (ctrl_.verbosity){
306 std::cout <<
"Warning: Transport failed, retrying with more steps. dt = "
307 << dt_transport/Opm::unit::year <<
" year.\n";
310 std::vector<double>& sat = state.
saturation();
311 for (
int i=0; i < mygrid_.numCells(); ++i){
312 sat[2*i] = saturation[i];
313 sat[2*i+1] = 1-saturation[i];
320 std::cout <<
"EulerUpstreamImplicite used " << repeats
321 <<
" repeats and " << nr_transport_steps <<
" steps"<< std::endl;
323 std::cout <<
"Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
326 std::vector<double>& sat = state.
saturation();
327 for (
int i=0; i < mygrid_.numCells(); ++i){
328 saturation[i] = sat[2*i];
332 std::cerr <<
"EulerUpstreamImplicit did not converge" << std::endl;
339 template <
class GI,
class RP,
class BC>
342 int num_cells = s.size();
343 for (
int cell = 0; cell < num_cells; ++cell) {
344 if (s[cell] > 1.0 || s[cell] < 0.0) {
346 s[cell] = std::max(
std::min(s[cell], 1.0), 0.0);
347 }
else if (s[cell] > 1.001 || s[cell] < -0.001) {
348 OPM_THROW(std::runtime_error,
349 "Saturation out of range in EulerUpstreamImplicit: "
350 "Cell " + std::to_string(cell) +
351 " sat " + std::to_string(s[cell]));
EulerUpstreamImplicit()
Definition: EulerUpstreamImplicit_impl.hpp:60
void checkAndPossiblyClampSat(std::vector< double > &s) const
Definition: EulerUpstreamImplicit_impl.hpp:340
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:217
void init(const Opm::ParameterGroup ¶m)
Definition: EulerUpstreamImplicit_impl.hpp:74
Definition: ImplicitTransport.hpp:105
void solve(const Grid &g, const SourceTerms *src, const double dt, const ImplicitTransportDetails::NRControl &ctrl, ReservoirState &state, LinearSolver &linsolve, ImplicitTransportDetails::NRReport &rpt)
Definition: ImplicitTransport.hpp:116
Definition: ImplicitTransportDefs.hpp:162
Definition: ImplicitTransportDefs.hpp:52
::std::vector< double > & faceflux()
Definition: ImplicitTransportDefs.hpp:63
::std::vector< double > & saturation()
Definition: ImplicitTransportDefs.hpp:64
Definition: SinglePointUpwindTwoPhase.hpp:278
void initGravityTrans(const Grid &g, const std::vector< double > &htrans)
Definition: SinglePointUpwindTwoPhase.hpp:474
void makefhfQPeriodic(const std::vector< int > &p_faces, const std::vector< int > &hf_faces, const std::vector< int > &nb_faces)
Definition: SinglePointUpwindTwoPhase.hpp:311
Definition: ImplicitTransportDefs.hpp:196
::std::vector< double > saturation
Definition: ImplicitTransportDefs.hpp:205
::std::vector< int > cell
Definition: ImplicitTransportDefs.hpp:202
int nsrc
Definition: ImplicitTransportDefs.hpp:200
::std::vector< double > flux
Definition: ImplicitTransportDefs.hpp:204
Definition: ImplicitTransportDefs.hpp:115
void mobility(int c, const Sat &s, Mob &mob, DMob &dmob) const
Definition: ImplicitTransportDefs.hpp:125
min[0]
Definition: elasticity_upscale_impl.hpp:146
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
Definition: ImplicitAssembly.hpp:43
Definition: ImplicitTransport.hpp:63
int flag
Definition: ImplicitTransport.hpp:65