Loading...
Searching...
No Matches
EulerUpstreamImplicit_impl.hpp
//===========================================================================
//
// File: EulerUpstreamImplicit_impl.hpp
//
// Created: Tue Jun 16 14:25:24 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B�rd Skaflestad <bard.skaflestad@sintef.no>
// Halvor M Nilsen <hnil@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
#define OPENRS_EULERUPSTREAMIMPLICIT_IMPL_HEADER
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/Average.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/grid/common/Volumes.hpp>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/porsol/common/ImplicitTransportDefs.hpp>
#include <opm/grid/transmissibility/trans_tpfa.h>
#include <cmath>
#include <algorithm>
#include <limits>
#include <vector>
#include <array>
#include <iostream>
namespace Opm
{
template <class GI, class RP, class BC>
inline EulerUpstreamImplicit<GI, RP, BC>::EulerUpstreamImplicit()
{
check_sat_ = true;
clamp_sat_ = false;
}
template <class GI, class RP, class BC>
inline EulerUpstreamImplicit<GI, RP, BC>::EulerUpstreamImplicit(const GI& g, const RP& r, const BC& b)
{
check_sat_ = true;
clamp_sat_ = false;
initObj(g, r, b);
}
template <class GI, class RP, class BC>
inline void EulerUpstreamImplicit<GI, RP, BC>::init(const Opm::ParameterGroup& param)
{
check_sat_ = param.getDefault("check_sat", check_sat_);
clamp_sat_ = param.getDefault("clamp_sat", clamp_sat_);
//Opm::ImplicitTransportDetails::NRControl ctrl_;
ctrl_.max_it = param.getDefault("transport_nr_max_it", 10);
max_repeats_ = param.getDefault("transport_max_rep", 10);
ctrl_.atol = param.getDefault("transport_atol", 1.0e-6);
ctrl_.rtol = param.getDefault("transport_rtol", 5.0e-7);
ctrl_.max_it_ls = param.getDefault("transport_max_it_ls", 20);
ctrl_.dxtol = param.getDefault("transport_dxtol", 1e-6);
ctrl_.verbosity = param.getDefault("transport_verbosity", 0);
}
template <class GI, class RP, class BC>
inline void EulerUpstreamImplicit<GI, RP, BC>::init(const Opm::ParameterGroup& param,
const GI& g, const RP& r, const BC& b)
{
init(param);
initObj(g, r, b);
}
template <class GI, class RP, class BC>
inline void EulerUpstreamImplicit<GI, RP, BC>::initObj(const GI& g, const RP& r, const BC& b)
{
//residual_computer_.initObj(g, r, b);
mygrid_.init(g.grid());
porevol_.resize(mygrid_.numCells());
for (int i = 0; i < mygrid_.numCells(); ++i){
porevol_[i]= mygrid_.cellVolume(i)*r.porosity(i);
}
// int numf=mygrid_.numFaces();
int num_cells = mygrid_.numCells();
int ngconn = mygrid_.c_grid()->cell_facepos[num_cells];
//std::vector<double> htrans_(ngconn);
htrans_.resize(ngconn);
const double* perm = &(r.permeability(0)(0,0));
tpfa_htrans_compute(mygrid_.c_grid(), perm, &htrans_[0]);
// int count = 0;
myrp_= r;
typedef typename GI::CellIterator CIt;
typedef typename CIt::FaceIterator FIt;
std::vector<FIt> bid_to_face;
int maxbid = 0;
for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
int bid = f->boundaryId();
maxbid = std::max(maxbid, bid);
}
}
bid_to_face.resize(maxbid + 1);
std::vector<int> egf_cf(mygrid_.numFaces());
int cix=0;
for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
int loc_fix=0;
for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
if (f->boundary() && b.satCond(*f).isPeriodic()) {
bid_to_face[f->boundaryId()] = f;
}
int egf=f->index();
int cf=mygrid_.cellFace(cix,loc_fix);
egf_cf[egf]=cf;
loc_fix+=1;
}
cix+=1;
}
#ifndef NDEBUG
const UnstructuredGrid& c_grid=*mygrid_.c_grid();
#endif
int hf_ind=0;
periodic_cells_.resize(0);
periodic_faces_.resize(0);
periodic_hfaces_.resize(0);
periodic_nbfaces_.resize(0);
//cell1 = cell0;
direclet_cells_.resize(0);
direclet_sat_.resize(0);
direclet_sat_.resize(0);
direclet_hfaces_.resize(0);
assert(periodic_cells_.size()==0);
for (CIt c = g.cellbegin(); c != g.cellend(); ++c) {
int cell0 = c->index();
for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
// Neighbour face, will be changed if on a periodic boundary.
// Compute cell[1], cell_sat[1]
FIt nbface = f;
if (f->boundary()) {
if (b.satCond(*f).isPeriodic()) {
nbface = bid_to_face[b.getPeriodicPartner(f->boundaryId())];
assert(nbface != f);
int cell1 = nbface->cellIndex();
assert(cell0 != cell1);
int f_ind=f->index();
int fn_ind=nbface->index();
// mapping face indices
f_ind=egf_cf[f_ind];
fn_ind=egf_cf[fn_ind];
assert((c_grid.face_cells[2*f_ind]==-1) || (c_grid.face_cells[2*f_ind+1]==-1));
assert((c_grid.face_cells[2*fn_ind]==-1) || (c_grid.face_cells[2*fn_ind+1]==-1));
assert((c_grid.face_cells[2*f_ind]==cell0) || (c_grid.face_cells[2*f_ind+1]==cell0));
assert((c_grid.face_cells[2*fn_ind]==cell1) || (c_grid.face_cells[2*fn_ind+1]==cell1));
periodic_cells_.push_back(cell0);
periodic_cells_.push_back(cell1);
periodic_faces_.push_back(f_ind);
periodic_hfaces_.push_back(hf_ind);
periodic_nbfaces_.push_back(fn_ind);
} else if (!( b.flowCond(*f).isNeumann() && b.flowCond(*f).outflux() == 0.0)) {
//cell1 = cell0;
direclet_cells_.push_back(cell0);
direclet_sat_.push_back(b.satCond(*f).saturation());
direclet_sat_.push_back(1-b.satCond(*f).saturation());//only work for 2 phases
direclet_hfaces_.push_back(hf_ind);
}
}
hf_ind+=1;
}
}
mygrid_.makeQPeriodic(periodic_hfaces_,periodic_cells_);
// use fractional flow instead of saturation as src
TwophaseFluid myfluid(myrp_);
int num_b=direclet_cells_.size();
for(int i=0; i <num_b; ++i){
std::array<double,2> sat = {{direclet_sat_[2*i] ,direclet_sat_[2*i+1] }};
std::array<double,2> mob;
std::array<double,2*2> dmob;
myfluid.mobility(direclet_cells_[i], sat, mob, dmob);
double fl = mob[0]/(mob[0]+mob[1]);
direclet_sat_[2*i] = fl;
direclet_sat_[2*i+1] = 1-fl;
}
}
template <class GI, class RP, class BC>
inline void EulerUpstreamImplicit<GI, RP, BC>::display()
{
using namespace std;
cout << endl;
cout <<"Displaying some members of EulerUpstreamImplicit" << endl;
cout << endl;
}
template <class GI, class RP, class BC>
template <class PressureSolution>
bool EulerUpstreamImplicit<GI, RP, BC>::transportSolve(std::vector<double>& saturation,
const double time,
const typename GI::Vector& /* gravity */,
const PressureSolution& pressure_sol,
const Opm::SparseVector<double>& /* injection_rates */) const
{
Opm::ReservoirState<2> state(mygrid_.c_grid());
{
std::vector<double>& sat = state.saturation();
for (int i=0; i < mygrid_.numCells(); ++i){
sat[2*i] = saturation[i];
sat[2*i+1] = 1-saturation[i];
}
}
//int count=0;
const UnstructuredGrid* cgrid = mygrid_.c_grid();
int numhf = cgrid->cell_facepos[cgrid->number_of_cells];
std::vector<double> faceflux(numhf);
{
auto i = grid_size_t(0);
for (int c = 0; c < cgrid->number_of_cells; ++c){
for (; i < cgrid->cell_facepos[c + 1]; ++i) {
int f= cgrid->cell_faces[i];
double outflux = pressure_sol.outflux(i);
double sgn = 2.0*(cgrid->face_cells[2*f + 0] == c) - 1;
faceflux[f] = sgn * outflux;
}
}
}
int num_db=direclet_hfaces_.size();
std::vector<double> sflux(num_db);
for (int i=0; i < num_db;++i){
sflux[i]=-pressure_sol.outflux(direclet_hfaces_[i]);
}
state.faceflux()=faceflux;
double dt_transport = time;
int nr_transport_steps = 1;
Opm::time::StopWatch clock;
int repeats = 0;
bool finished = false;
clock.start();
TwophaseFluid myfluid(myrp_);
double* tmp_grav=0;
const UnstructuredGrid& c_grid=*mygrid_.c_grid();
TransportModel model(myfluid,c_grid,porevol_,tmp_grav);
model.makefhfQPeriodic(periodic_faces_,periodic_hfaces_, periodic_nbfaces_);
model.initGravityTrans(*mygrid_.c_grid(),htrans_);
TransportSolver tsolver(model);
LinearSolver linsolve_;
Opm::ImplicitTransportDetails::NRReport rpt_;
Opm::TransportSource tsrc;//create_transport_source(0, 2);
// the input flux is assumed to be the saturation times the flux in the transport solver
tsrc.nsrc =direclet_cells_.size();
tsrc.saturation = direclet_sat_;
tsrc.cell = direclet_cells_;
tsrc.flux = sflux;
while (!finished) {
for (int q = 0; q < nr_transport_steps; ++q) {
tsolver.solve(*mygrid_.c_grid(), &tsrc, dt_transport, ctrl_, state, linsolve_, rpt_);
if(rpt_.flag<0){
break;
}
}
if(!(rpt_.flag<0) ){
finished =true;
}else{
if(repeats >max_repeats_){
finished=true;
}else{
OPM_MESSAGE("Warning: Transport failed, retrying with more steps.");
nr_transport_steps *= 2;
dt_transport = time/nr_transport_steps;
if (ctrl_.verbosity){
std::cout << "Warning: Transport failed, retrying with more steps. dt = "
<< dt_transport/Opm::unit::year << " year.\n";
}
std::vector<double>& sat = state.saturation();
for (int i=0; i < mygrid_.numCells(); ++i){
sat[2*i] = saturation[i];
sat[2*i+1] = 1-saturation[i];
}
}
}
repeats +=1;
}
clock.stop();
std::cout << "EulerUpstreamImplicite used " << repeats
<< " repeats and " << nr_transport_steps <<" steps"<< std::endl;
#ifdef VERBOSE
std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
#endif // VERBOSE
{
const std::vector<double>& sat = state.saturation();
for (int i=0; i < mygrid_.numCells(); ++i){
saturation[i] = sat[2*i];
}
}
if((rpt_.flag<0)){
std::cerr << "EulerUpstreamImplicit did not converge" << std::endl;
return false;
}else{
return true;
}
}
template <class GI, class RP, class BC>
inline void EulerUpstreamImplicit<GI, RP, BC>::checkAndPossiblyClampSat(std::vector<double>& s) const
{
int num_cells = s.size();
for (int cell = 0; cell < num_cells; ++cell) {
if (s[cell] > 1.0 || s[cell] < 0.0) {
if (clamp_sat_) {
s[cell] = std::max(std::min(s[cell], 1.0), 0.0);
} else if (s[cell] > 1.001 || s[cell] < -0.001) {
OPM_THROW(std::runtime_error,
"Saturation out of range in EulerUpstreamImplicit: "
"Cell " + std::to_string(cell) +
" sat " + std::to_string(s[cell]));
}
}
}
}
}
#endif // OPENRS_EULERUPSTREAM_IMPL_HEADER
Opm::EulerUpstreamImplicit::EulerUpstreamImplicit
EulerUpstreamImplicit()
Definition EulerUpstreamImplicit_impl.hpp:60
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.
Opm::EulerUpstreamImplicit::display
void display()
Definition EulerUpstreamImplicit_impl.hpp:217
Opm::EulerUpstreamImplicit::init
void init(const Opm::ParameterGroup &param)
Definition EulerUpstreamImplicit_impl.hpp:74
Opm::ReservoirState
Definition ImplicitTransportDefs.hpp:52
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43