My Project
EulerUpstreamImplicit_impl.hpp
1 //===========================================================================
2 //
3 // File: EulerUpstreamImplicit_impl.hpp
4 //
5 // Created: Tue Jun 16 14:25:24 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // B�rd Skaflestad <bard.skaflestad@sintef.no>
9 // Halvor M Nilsen <hnil@sintef.no>
10 //
11 // $Date$
12 //
13 // $Revision$
14 //
15 //===========================================================================
16 
17 /*
18  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
19  Copyright 2009, 2010 Statoil ASA.
20 
21  This file is part of The Open Reservoir Simulator Project (OpenRS).
22 
23  OpenRS is free software: you can redistribute it and/or modify
24  it under the terms of the GNU General Public License as published by
25  the Free Software Foundation, either version 3 of the License, or
26  (at your option) any later version.
27 
28  OpenRS is distributed in the hope that it will be useful,
29  but WITHOUT ANY WARRANTY; without even the implied warranty of
30  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31  GNU General Public License for more details.
32 
33  You should have received a copy of the GNU General Public License
34  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
35 */
36 
37 #ifndef OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
38 #define OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
39 
40 #include <opm/common/ErrorMacros.hpp>
41 #include <opm/core/utility/Average.hpp>
42 #include <opm/input/eclipse/Units/Units.hpp>
43 #include <opm/grid/common/Volumes.hpp>
44 #include <opm/grid/utility/StopWatch.hpp>
45 #include <opm/porsol/common/ImplicitTransportDefs.hpp>
46 #include <opm/grid/transmissibility/trans_tpfa.h>
47 
48 #include <cmath>
49 #include <algorithm>
50 #include <limits>
51 #include <vector>
52 #include <array>
53 #include <iostream>
54 
55 namespace Opm
56 {
57 
58 
59  template <class GI, class RP, class BC>
61  {
62  check_sat_ = true;
63  clamp_sat_ = false;
64  }
65  template <class GI, class RP, class BC>
66  inline EulerUpstreamImplicit<GI, RP, BC>::EulerUpstreamImplicit(const GI& g, const RP& r, const BC& b)
67  {
68  check_sat_ = true;
69  clamp_sat_ = false;
70  initObj(g, r, b);
71  }
72 
73  template <class GI, class RP, class BC>
74  inline void EulerUpstreamImplicit<GI, RP, BC>::init(const Opm::ParameterGroup& param)
75  {
76  check_sat_ = param.getDefault("check_sat", check_sat_);
77  clamp_sat_ = param.getDefault("clamp_sat", clamp_sat_);
78  //Opm::ImplicitTransportDetails::NRControl ctrl_;
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);
86  }
87 
88  template <class GI, class RP, class BC>
89  inline void EulerUpstreamImplicit<GI, RP, BC>::init(const Opm::ParameterGroup& param,
90  const GI& g, const RP& r, const BC& b)
91  {
92  init(param);
93  initObj(g, r, b);
94  }
95 
96 
97  template <class GI, class RP, class BC>
98  inline void EulerUpstreamImplicit<GI, RP, BC>::initObj(const GI& g, const RP& r, const BC& b)
99  {
100  //residual_computer_.initObj(g, r, b);
101 
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);
106  }
107  // int numf=mygrid_.numFaces();
108  int num_cells = mygrid_.numCells();
109  int ngconn = mygrid_.c_grid()->cell_facepos[num_cells];
110  //std::vector<double> htrans_(ngconn);
111  htrans_.resize(ngconn);
112  const double* perm = &(r.permeability(0)(0,0));
113  tpfa_htrans_compute(mygrid_.c_grid(), perm, &htrans_[0]);
114  // int count = 0;
115 
116  myrp_= r;
117 
118  typedef typename GI::CellIterator CIt;
119  typedef typename CIt::FaceIterator FIt;
120  std::vector<FIt> bid_to_face;
121  int maxbid = 0;
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);
126  }
127  }
128 
129  bid_to_face.resize(maxbid + 1);
130  std::vector<int> egf_cf(mygrid_.numFaces());
131  int cix=0;
132  for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
133  int loc_fix=0;
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;
137  }
138  int egf=f->index();
139  int cf=mygrid_.cellFace(cix,loc_fix);
140  egf_cf[egf]=cf;
141  loc_fix+=1;
142  }
143  cix+=1;
144  }
145 
146 #ifndef NDEBUG
147  const UnstructuredGrid& c_grid=*mygrid_.c_grid();
148 #endif
149  int hf_ind=0;
150  periodic_cells_.resize(0);
151  periodic_faces_.resize(0);
152  periodic_hfaces_.resize(0);
153  periodic_nbfaces_.resize(0);
154  //cell1 = cell0;
155  direclet_cells_.resize(0);
156  direclet_sat_.resize(0);
157  direclet_sat_.resize(0);
158  direclet_hfaces_.resize(0);
159 
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) {
164  // Neighbour face, will be changed if on a periodic boundary.
165  // Compute cell[1], cell_sat[1]
166  FIt nbface = f;
167  if (f->boundary()) {
168  if (b.satCond(*f).isPeriodic()) {
169  nbface = bid_to_face[b.getPeriodicPartner(f->boundaryId())];
170  assert(nbface != f);
171  int cell1 = nbface->cellIndex();
172  assert(cell0 != cell1);
173 
174  int f_ind=f->index();
175 
176  int fn_ind=nbface->index();
177  // mapping face indices
178  f_ind=egf_cf[f_ind];
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)) {
190  //cell1 = cell0;
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());//only work for 2 phases
194  direclet_hfaces_.push_back(hf_ind);
195  }
196  }
197  hf_ind+=1;
198  }
199  }
200 
201  mygrid_.makeQPeriodic(periodic_hfaces_,periodic_cells_);
202  // use fractional flow instead of saturation as src
203  TwophaseFluid myfluid(myrp_);
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;
213  }
214  }
215 
216  template <class GI, class RP, class BC>
218  {
219  using namespace std;
220  cout << endl;
221  cout <<"Displaying some members of EulerUpstreamImplicit" << endl;
222  cout << endl;
223  }
224 
225  template <class GI, class RP, class BC>
226  template <class PressureSolution>
227  bool EulerUpstreamImplicit<GI, RP, BC>::transportSolve(std::vector<double>& saturation,
228  const double time,
229  const typename GI::Vector& /* gravity */,
230  const PressureSolution& pressure_sol,
231  const Opm::SparseVector<double>& /* injection_rates */) const
232  {
233 
234  Opm::ReservoirState<2> state(mygrid_.c_grid());
235  {
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];
240  }
241  }
242 
243  //int count=0;
244  const UnstructuredGrid* cgrid = mygrid_.c_grid();
245  int numhf = cgrid->cell_facepos[cgrid->number_of_cells];
246 
247  std::vector<double> faceflux(numhf);
248 
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;
255  }
256  }
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]);
261  }
262  state.faceflux()=faceflux;
263 
264  double dt_transport = time;
265  int nr_transport_steps = 1;
266  Opm::time::StopWatch clock;
267  int repeats = 0;
268  bool finished = false;
269  clock.start();
270 
271  TwophaseFluid myfluid(myrp_);
272  double* tmp_grav=0;
273  const UnstructuredGrid& c_grid=*mygrid_.c_grid();
274  TransportModel model(myfluid,c_grid,porevol_,tmp_grav);
275  model.makefhfQPeriodic(periodic_faces_,periodic_hfaces_, periodic_nbfaces_);
276  model.initGravityTrans(*mygrid_.c_grid(),htrans_);
277  TransportSolver tsolver(model);
278  LinearSolver linsolve_;
280 
281  Opm::TransportSource tsrc;//create_transport_source(0, 2);
282  // the input flux is assumed to be the saturation times the flux in the transport solver
283 
284  tsrc.nsrc =direclet_cells_.size();
285  tsrc.saturation = direclet_sat_;
286  tsrc.cell = direclet_cells_;
287  tsrc.flux = sflux;
288 
289  while (!finished) {
290  for (int q = 0; q < nr_transport_steps; ++q) {
291  tsolver.solve(*mygrid_.c_grid(), &tsrc, dt_transport, ctrl_, state, linsolve_, rpt_);
292  if(rpt_.flag<0){
293  break;
294  }
295  }
296  if(!(rpt_.flag<0) ){
297  finished =true;
298  }else{
299  if(repeats >max_repeats_){
300  finished=true;
301  }else{
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";
308  }
309 
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];
314  }
315  }
316  }
317  repeats +=1;
318  }
319  clock.stop();
320  std::cout << "EulerUpstreamImplicite used " << repeats
321  << " repeats and " << nr_transport_steps <<" steps"<< std::endl;
322 #ifdef VERBOSE
323  std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
324 #endif // VERBOSE
325  {
326  std::vector<double>& sat = state.saturation();
327  for (int i=0; i < mygrid_.numCells(); ++i){
328  saturation[i] = sat[2*i];
329  }
330  }
331  if((rpt_.flag<0)){
332  std::cerr << "EulerUpstreamImplicit did not converge" << std::endl;
333  return false;
334  }else{
335  return true;
336  }
337  }
338 
339  template <class GI, class RP, class BC>
340  inline void EulerUpstreamImplicit<GI, RP, BC>::checkAndPossiblyClampSat(std::vector<double>& s) const
341  {
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) {
345  if (clamp_sat_) {
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, "Saturation out of range in EulerUpstreamImplicit: Cell " << cell << " sat " << s[cell]);
349  }
350  }
351  }
352  }
353 
354 }
355 
356 
357 
358 #endif // OPENRS_EULERUPSTREAM_IMPL_HEADER
EulerUpstreamImplicit()
Definition: EulerUpstreamImplicit_impl.hpp:60
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 &param)
Definition: EulerUpstreamImplicit_impl.hpp:74
Definition: ImplicitTransportDefs.hpp:52
Definition: ImplicitTransportDefs.hpp:202
Definition: ImplicitTransportDefs.hpp:115
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
Definition: ImplicitTransport.hpp:63