EulerUpstreamImplicit_impl.hpp
Go to the documentation of this file.
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/core/utility/Units.hpp>
43#include <dune/grid/common/Volumes.hpp>
44#include <opm/core/utility/StopWatch.hpp>
45#include <opm/porsol/common/ImplicitTransportDefs.hpp>
46#include <opm/core/pressure/tpfa/trans_tpfa.h>
47
48#include <cmath>
49#include <algorithm>
50#include <limits>
51#include <vector>
52#include <array>
53#include <iostream>
54
55namespace Opm
56{
57
58
59 template <class GI, class RP, class BC>
60 inline EulerUpstreamImplicit<GI, RP, BC>::EulerUpstreamImplicit()
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::parameter::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::parameter::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 int bf_ind=0;
151 periodic_cells_.resize(0);
152 periodic_faces_.resize(0);
153 periodic_hfaces_.resize(0);
154 periodic_nbfaces_.resize(0);
155 //cell1 = cell0;
156 direclet_cells_.resize(0);
157 direclet_sat_.resize(0);
158 direclet_sat_.resize(0);
159 direclet_hfaces_.resize(0);
160
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) {
165 // Neighbour face, will be changed if on a periodic boundary.
166 // Compute cell[1], cell_sat[1]
167 FIt nbface = f;
168 if (f->boundary()) {
169 bf_ind+=1;
170 if (b.satCond(*f).isPeriodic()) {
171 nbface = bid_to_face[b.getPeriodicPartner(f->boundaryId())];
172 assert(nbface != f);
173 int cell1 = nbface->cellIndex();
174 assert(cell0 != cell1);
175
176 int f_ind=f->index();
177
178 int fn_ind=nbface->index();
179 // mapping face indices
180 f_ind=egf_cf[f_ind];
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)) {
192 //cell1 = cell0;
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());//only work for 2 phases
196 direclet_hfaces_.push_back(hf_ind);
197 }
198 }
199 hf_ind+=1;
200 }
201 }
202
203 mygrid_.makeQPeriodic(periodic_hfaces_,periodic_cells_);
204 // use fractional flow instead of saturation as src
205 TwophaseFluid myfluid(myrp_);
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;
215 }
216 }
217
218 template <class GI, class RP, class BC>
219 inline void EulerUpstreamImplicit<GI, RP, BC>::display()
220 {
221 using namespace std;
222 cout << endl;
223 cout <<"Displaying some members of EulerUpstreamImplicit" << endl;
224 cout << endl;
225 }
226
227 template <class GI, class RP, class BC>
228 template <class PressureSolution>
229 bool EulerUpstreamImplicit<GI, RP, BC>::transportSolve(std::vector<double>& saturation,
230 const double time,
231 const typename GI::Vector& gravity,
232 const PressureSolution& pressure_sol,
233 const Opm::SparseVector<double>& injection_rates) const
234 {
235
236 Opm::ReservoirState<2> state(mygrid_.c_grid());
237 {
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];
242 }
243 }
244
245 //int count=0;
246 const UnstructuredGrid* cgrid = mygrid_.c_grid();
247 int numhf = cgrid->cell_facepos[cgrid->number_of_cells];
248
249 std::vector<double> faceflux(numhf);
250
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;
257 }
258 }
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]);
263 }
264 state.faceflux()=faceflux;
265
266 double dt_transport = time;
267 int nr_transport_steps = 1;
268 Opm::time::StopWatch clock;
269 int repeats = 0;
270 bool finished = false;
271 clock.start();
272
273 TwophaseFluid myfluid(myrp_);
274 double* tmp_grav=0;
275 const UnstructuredGrid& c_grid=*mygrid_.c_grid();
276 TransportModel model(myfluid,c_grid,porevol_,tmp_grav);
277 model.makefhfQPeriodic(periodic_faces_,periodic_hfaces_, periodic_nbfaces_);
278 model.initGravityTrans(*mygrid_.c_grid(),htrans_);
279 TransportSolver tsolver(model);
280 LinearSolver linsolve_;
281 Opm::ImplicitTransportDetails::NRReport rpt_;
282
283 Opm::TransportSource tsrc;//create_transport_source(0, 2);
284 // the input flux is assumed to be the saturation times the flux in the transport solver
285
286 tsrc.nsrc =direclet_cells_.size();
287 tsrc.saturation = direclet_sat_;
288 tsrc.cell = direclet_cells_;
289 tsrc.flux = sflux;
290
291 while (!finished) {
292 for (int q = 0; q < nr_transport_steps; ++q) {
293 tsolver.solve(*mygrid_.c_grid(), &tsrc, dt_transport, ctrl_, state, linsolve_, rpt_);
294 if(rpt_.flag<0){
295 break;
296 }
297 }
298 if(!(rpt_.flag<0) ){
299 finished =true;
300 }else{
301 if(repeats >max_repeats_){
302 finished=true;
303 }else{
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";
310 }
311
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];
316 }
317 }
318 }
319 repeats +=1;
320 }
321 clock.stop();
322 std::cout << "EulerUpstreamImplicite used " << repeats
323 << " repeats and " << nr_transport_steps <<" steps"<< std::endl;
324#ifdef VERBOSE
325 std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
326#endif // VERBOSE
327 {
328 std::vector<double>& sat = state.saturation();
329 for (int i=0; i < mygrid_.numCells(); ++i){
330 saturation[i] = sat[2*i];
331 }
332 }
333 if((rpt_.flag<0)){
334 std::cerr << "EulerUpstreamImplicit did not converge" << std::endl;
335 return false;
336 }else{
337 return true;
338 }
339 }
340
341 template <class GI, class RP, class BC>
342 inline void EulerUpstreamImplicit<GI, RP, BC>::checkAndPossiblyClampSat(std::vector<double>& s) const
343 {
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) {
347 if (clamp_sat_) {
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]);
351 }
352 }
353 }
354 }
355
356}
357
358
359
360#endif // OPENRS_EULERUPSTREAM_IMPL_HEADER
Opm::EulerUpstreamImplicit::EulerUpstreamImplicit
EulerUpstreamImplicit()
Definition: EulerUpstreamImplicit_impl.hpp:60
Opm::EulerUpstreamImplicit::checkAndPossiblyClampSat
void checkAndPossiblyClampSat(std::vector< double > &s) const
Definition: EulerUpstreamImplicit_impl.hpp:342
Opm::EulerUpstreamImplicit::TransportModel
Opm::SinglePointUpwindTwoPhase< TwophaseFluid > TransportModel
Definition: EulerUpstreamImplicit.hpp:123
Opm::EulerUpstreamImplicit::TransportSolver
Opm::ImplicitTransport< TransportModel, JacSys, Opm::MaxNormDune, Opm::ImplicitTransportDefault::VectorNegater, Opm::ImplicitTransportDefault::VectorZero, Opm::ImplicitTransportDefault::MatrixZero, Opm::ImplicitTransportDefault::VectorAssign > TransportSolver
Definition: EulerUpstreamImplicit.hpp:146
Opm::EulerUpstreamImplicit::initObj
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary)
Definition: EulerUpstreamImplicit_impl.hpp:98
Opm::EulerUpstreamImplicit::transportSolve
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.
Opm::EulerUpstreamImplicit::display
void display()
Definition: EulerUpstreamImplicit_impl.hpp:219
Opm::EulerUpstreamImplicit::init
void init(const Opm::parameter::ParameterGroup &param)
Definition: EulerUpstreamImplicit_impl.hpp:74
Opm::LinearSolverBICGSTAB
Definition: ImplicitTransportDefs.hpp:174
Opm::ReservoirState
Definition: ImplicitTransportDefs.hpp:57
Opm::ReservoirState::faceflux
::std::vector< double > & faceflux()
Definition: ImplicitTransportDefs.hpp:68
Opm::ReservoirState::saturation
::std::vector< double > & saturation()
Definition: ImplicitTransportDefs.hpp:69
Opm::TransportSource
Definition: ImplicitTransportDefs.hpp:237
Opm::TransportSource::saturation
::std::vector< double > saturation
Definition: ImplicitTransportDefs.hpp:246
Opm::TransportSource::cell
::std::vector< int > cell
Definition: ImplicitTransportDefs.hpp:243
Opm::TransportSource::nsrc
int nsrc
Definition: ImplicitTransportDefs.hpp:241
Opm::TransportSource::flux
::std::vector< double > flux
Definition: ImplicitTransportDefs.hpp:245
Opm::TwophaseFluidWrapper
Definition: ImplicitTransportDefs.hpp:120
Opm::TwophaseFluidWrapper::mobility
void mobility(int c, const Sat &s, Mob &mob, DMob &dmob) const
Definition: ImplicitTransportDefs.hpp:139
Opm
Definition: BlackoilFluid.hpp:32