mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-26 03:00:17 -06:00
714 lines
28 KiB
C++
714 lines
28 KiB
C++
/*===========================================================================
|
|
//
|
|
// File: SinglePointUpwindTwoPhase.hpp
|
|
//
|
|
// Created: 2011-09-28 14:21:34+0200
|
|
//
|
|
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
|
|
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
|
|
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
|
|
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
|
|
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
|
|
//
|
|
//==========================================================================*/
|
|
|
|
|
|
/*
|
|
Copyright 2011 SINTEF ICT, Applied Mathematics.
|
|
Copyright 2011 Statoil ASA.
|
|
|
|
This file is part of the Open Porous Media Project (OPM).
|
|
|
|
OPM 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.
|
|
|
|
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#ifndef OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
|
|
#define OPM_SINGLEPOINTUPWINDTWOPHASEPOLYMER_HPP_HEADER
|
|
|
|
#include <cassert>
|
|
#include <cstddef>
|
|
|
|
#include <algorithm>
|
|
#include <vector>
|
|
#include <iostream>
|
|
|
|
namespace Opm {
|
|
namespace polymer_reorder {
|
|
class ModelParameterStorage {
|
|
public:
|
|
ModelParameterStorage(int nc, int totconn)
|
|
: drho_(0.0), rockdensity_(0.0), mob_(0),
|
|
dmobds_(0), dmobwatdc_(0), mc_(0),
|
|
dmcdc_(0), porevol_(0), porosity_(0), dg_(0), sw_(0), c_(0), cmax_(0),
|
|
ds_(0), dsc_(0), dcads_(0), dcadsdc_(0), pc_(0), dpc_(0),
|
|
trans_(0), data_()
|
|
{
|
|
size_t alloc_sz;
|
|
|
|
alloc_sz = 2 * nc; // mob_
|
|
alloc_sz += 2 * nc; // dmobds_
|
|
alloc_sz += nc; // dmobwatdc_
|
|
alloc_sz += nc; // mc_
|
|
alloc_sz += nc; // dmcdc_
|
|
alloc_sz += 1 * nc; // porevol_
|
|
alloc_sz += 1 * nc; // porosity_
|
|
alloc_sz += 1 * totconn; // dg_
|
|
alloc_sz += 1 * nc; // sw_
|
|
alloc_sz += 1 * nc; // c_
|
|
alloc_sz += 1 * nc; // cmax_
|
|
alloc_sz += 1 * nc; // ds_
|
|
alloc_sz += 1 * nc; // dsc_
|
|
alloc_sz += 1 * nc; // dcads_
|
|
alloc_sz += 1 * nc; // dcadsdc_
|
|
alloc_sz += 1 * nc; // pc_
|
|
alloc_sz += 1 * nc; // dpc_
|
|
alloc_sz += 1 * totconn; // trans_
|
|
data_.resize(alloc_sz);
|
|
|
|
mob_ = &data_[0];
|
|
dmobds_ = mob_ + (2 * nc );
|
|
dmobwatdc_ = dmobds_ + (2 * nc );
|
|
mc_ = dmobwatdc_ + (1 * nc );
|
|
dmcdc_ = mc_ + (1 * nc );
|
|
porevol_ = dmcdc_ + (1 * nc );
|
|
porosity_ = porevol_ + (1 * nc );
|
|
dg_ = porosity_ + (1 * nc );
|
|
sw_ = dg_ + (1 * totconn);
|
|
c_ = sw_ + (1 * nc );
|
|
cmax_ = c_ + (1 * nc );
|
|
ds_ = cmax_ + (1 * nc );
|
|
dsc_ = ds_ + (1 * nc );
|
|
dcads_ = dsc_ + (1 * nc );
|
|
dcadsdc_ = dcads_ + (1 * nc );
|
|
pc_ = dcadsdc_ + (1 * nc );
|
|
dpc_ = pc_ + (1 * nc );
|
|
trans_ = dpc_ + (1 * nc );
|
|
}
|
|
|
|
double& drho () { return drho_ ; }
|
|
double drho () const { return drho_ ; }
|
|
|
|
double& rockdensity() { return rockdensity_ ; }
|
|
double rockdensity() const { return rockdensity_ ; }
|
|
|
|
double* mob (int cell) { return mob_ + (2*cell + 0); }
|
|
const double* mob (int cell) const { return mob_ + (2*cell + 0); }
|
|
|
|
double* dmobds (int cell) { return dmobds_ + (2*cell + 0); }
|
|
const double* dmobds (int cell) const { return dmobds_ + (2*cell + 0); }
|
|
|
|
double& dmobwatdc (int cell) { return dmobwatdc_[cell]; }
|
|
double dmobwatdc (int cell) const { return dmobwatdc_[cell]; }
|
|
|
|
double& mc (int cell) { return mc_[cell]; }
|
|
double mc (int cell) const { return mc_[cell]; }
|
|
|
|
double& dmcdc (int cell) { return dmcdc_[cell]; }
|
|
double dmcdc (int cell) const { return dmcdc_[cell]; }
|
|
|
|
double* porevol() { return porevol_ ; }
|
|
double porevol(int cell) const { return porevol_[cell] ; }
|
|
|
|
double* porosity() { return porosity_ ; }
|
|
double porosity(int cell) const { return porosity_[cell] ; }
|
|
|
|
|
|
double& dg(int i) { return dg_[i] ; }
|
|
double dg(int i) const { return dg_[i] ; }
|
|
|
|
double& sw(int cell) { return sw_[cell] ; }
|
|
double sw(int cell) const { return sw_[cell] ; }
|
|
|
|
double& c(int cell) { return c_[cell] ; }
|
|
double c(int cell) const { return c_[cell] ; }
|
|
|
|
double& cmax(int cell) { return cmax_[cell] ; }
|
|
double cmax(int cell) const { return cmax_[cell] ; }
|
|
|
|
double& ds(int cell) { return ds_[cell] ; }
|
|
double ds(int cell) const { return ds_[cell] ; }
|
|
|
|
double& dsc(int cell) { return dsc_[cell] ; }
|
|
double dsc(int cell) const { return dsc_[cell] ; }
|
|
|
|
double& dcads(int cell) { return dcads_[cell] ; }
|
|
double dcads(int cell) const { return dcads_[cell] ; }
|
|
|
|
double& dcadsdc(int cell) { return dcadsdc_[cell] ; }
|
|
double dcadsdc(int cell) const { return dcadsdc_[cell] ; }
|
|
|
|
double& pc(int cell) { return pc_[cell] ; }
|
|
double pc(int cell) const { return pc_[cell] ; }
|
|
|
|
double& dpc(int cell) { return dpc_[cell] ; }
|
|
double dpc(int cell) const { return dpc_[cell] ; }
|
|
|
|
double& trans(int f) { return trans_[f] ; }
|
|
double trans(int f) const { return trans_[f] ; }
|
|
|
|
private:
|
|
double drho_ ;
|
|
double rockdensity_ ;
|
|
double *mob_ ;
|
|
double *dmobds_ ;
|
|
double *dmobwatdc_ ;
|
|
double *mc_ ;
|
|
double *dmcdc_ ;
|
|
double *porevol_ ;
|
|
double *porosity_ ;
|
|
double *dg_ ;
|
|
double *sw_ ;
|
|
double *c_ ;
|
|
double *cmax_ ;
|
|
double *ds_ ;
|
|
double *dsc_ ;
|
|
double *dcads_ ; // difference of cads to compute residual
|
|
double *dcadsdc_ ; // derivative of cads
|
|
double *pc_ ;
|
|
double *dpc_ ;
|
|
double *trans_ ;
|
|
|
|
std::vector<double> data_;
|
|
};
|
|
}
|
|
|
|
|
|
template <class TwophaseFluidPolymer>
|
|
class SinglePointUpwindTwoPhasePolymer {
|
|
public:
|
|
template <class Grid>
|
|
SinglePointUpwindTwoPhasePolymer(const TwophaseFluidPolymer& fluid ,
|
|
const Grid& g ,
|
|
const std::vector<double>& porevol ,
|
|
const double* grav = 0,
|
|
const bool guess_previous = true)
|
|
: fluid_ (fluid) ,
|
|
gravity_ (grav) ,
|
|
f2hf_ (2 * g.number_of_faces, -1) ,
|
|
store_ (g.number_of_cells,
|
|
g.cell_facepos[ g.number_of_cells ]),
|
|
init_step_use_previous_sol_(guess_previous) ,
|
|
sat_tol_ (1e-5)
|
|
{
|
|
|
|
if (gravity_) {
|
|
store_.drho() = fluid_.density(0) - fluid_.density(1);
|
|
}
|
|
|
|
for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
|
|
for (; i < g.cell_facepos[c + 1]; ++i) {
|
|
const int f = g.cell_faces[i];
|
|
const int p = 1 - (g.face_cells[2*f + 0] == c);
|
|
f2hf_[2*f + p] = i;
|
|
}
|
|
}
|
|
|
|
std::copy(porevol.begin(), porevol.end(), store_.porevol());
|
|
const double* poro = fluid.porosity();
|
|
std::copy(poro, poro + g.number_of_cells, store_.porosity());
|
|
store_.rockdensity() = fluid.rockdensity();
|
|
}
|
|
|
|
void makefhfQPeriodic( const std::vector<int>& p_faces,const std::vector<int>& hf_faces,
|
|
const std::vector<int>& nb_faces)
|
|
{
|
|
std::vector<int> nbhf(hf_faces.size());
|
|
for(unsigned int i=0; i<p_faces.size(); ++i){
|
|
int nbf = nb_faces[i];
|
|
if(f2hf_[2*nbf] == -1){
|
|
nbhf[i] = f2hf_[2*nbf+1];
|
|
}else{
|
|
assert(f2hf_[2*nbf+1]==-1);
|
|
nbhf[i] = f2hf_[2*nbf];
|
|
}
|
|
}
|
|
for(unsigned int i=0; i<p_faces.size(); ++i){
|
|
|
|
int f = p_faces[i];
|
|
int hf = hf_faces[i];
|
|
bool changed=false;
|
|
|
|
if(f2hf_[2*f] == hf){
|
|
assert(f2hf_[2*f+1]==-1);
|
|
}else{
|
|
assert(f2hf_[2*f]==-1);
|
|
f2hf_[2*f]=nbhf[i];
|
|
changed=true;
|
|
}
|
|
if(!changed){
|
|
if(f2hf_[2*f+1]== hf){
|
|
assert(f2hf_[2*f]==-1);
|
|
}else{
|
|
assert(f2hf_[2*f+1]==-1);
|
|
f2hf_[2*f+1]=nbhf[i];
|
|
changed=true;
|
|
}
|
|
}
|
|
assert(changed);
|
|
}
|
|
}
|
|
|
|
// -----------------------------------------------------------------
|
|
// System assembly innards
|
|
// -----------------------------------------------------------------
|
|
|
|
enum { DofPerCell = 1 };
|
|
|
|
void
|
|
initResidual(const int c, double* Fs, double* Fc) const {
|
|
(void) c; // Suppress 'unused' warning
|
|
*Fs = 0.0;
|
|
*Fc = 0.0;
|
|
}
|
|
|
|
template <class ReservoirState,
|
|
class Grid >
|
|
void
|
|
fluxConnection(const ReservoirState& state ,
|
|
const Grid& g ,
|
|
const double dt ,
|
|
const int cell ,
|
|
const int f ,
|
|
double* F , // F[0] = s-residual, F[1] = c-residual
|
|
double* dFd1 , //Jacobi matrix for residual with respect to variables in cell
|
|
double* dFd2 //Jacobi matrix for residual with respect to variables in OTHER cell
|
|
//dFd1[0]= d(F[0])/d(s1), dFd1[1]= d(F[0])/d(c1), dFd1[2]= d(F[1])/d(s1), dFd1[3]= d(F[1])/d(c1),
|
|
//dFd2[0]= d(F[0])/d(s2), dFd2[1]= d(F[0])/d(c2), dFd2[2]= d(F[1])/d(s2), dFd2[3]= d(F[1])/d(c2).
|
|
|
|
) const {
|
|
|
|
const int *n = g.face_cells + (2 * f);
|
|
double dflux = state.faceflux()[f];
|
|
double gflux = gravityFlux(f);
|
|
double pcflux, dpcflux[2];
|
|
capFlux(f, n, pcflux, dpcflux);
|
|
gflux += pcflux;
|
|
|
|
int pix[2];
|
|
double m[2], dmds[2], dmobwatdc;
|
|
double mc, dmcdc;
|
|
upwindMobility(dflux, gflux, n, pix, m, dmds, dmobwatdc, mc, dmcdc);
|
|
|
|
assert ((m[0] >= 0.0) && (m[1] >= 0.0));
|
|
|
|
double mt = m[0] + m[1];
|
|
assert (mt >= 0.0);
|
|
|
|
double sgn = 2.0*(n[0] == cell) - 1.0;
|
|
dflux *= sgn;
|
|
gflux *= sgn;
|
|
|
|
|
|
double f1 = m[0] / mt;
|
|
const double v1 = dflux + m[1]*gflux;
|
|
|
|
// Assemble residual contributions
|
|
F[0] += dt * f1 * v1;
|
|
F[1] += dt * mc * f1 * v1;
|
|
|
|
// Assemble Jacobian (J1 <-> cell, J2 <-> other)
|
|
double *dFsds[2];
|
|
double *dFsdc[2];
|
|
double *dFcds[2];
|
|
double *dFcdc[2];
|
|
if (n[0] == cell) {
|
|
dFsds[0] = &dFd1[0]; dFsds[1] = &dFd2[0];
|
|
dFsdc[0] = &dFd1[1]; dFsdc[1] = &dFd2[1];
|
|
dFcds[0] = &dFd1[2]; dFcds[1] = &dFd2[2];
|
|
dFcdc[0] = &dFd1[3]; dFcdc[1] = &dFd2[3];
|
|
// sign is positive
|
|
dFd1[0] += sgn*dt * f1 * dpcflux[0] * m[1];
|
|
dFd2[0] += sgn*dt * f1 * dpcflux[1] * m[1];
|
|
dFd1[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
|
|
dFd2[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
|
|
// We assume that the capillary pressure is independent of the polymer concentration.
|
|
// Hence, no more contributions.
|
|
} else {
|
|
dFsds[0] = &dFd2[0]; dFsds[1] = &dFd1[0];
|
|
dFsdc[0] = &dFd2[1]; dFsdc[1] = &dFd1[1];
|
|
dFcds[0] = &dFd2[2]; dFcds[1] = &dFd1[2];
|
|
dFcdc[0] = &dFd2[3]; dFcdc[1] = &dFd1[3];
|
|
// sign is negative
|
|
dFd1[0] += sgn*dt * f1 * dpcflux[1] * m[1];
|
|
dFd2[0] += sgn*dt * f1 * dpcflux[0] * m[1];
|
|
dFd1[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
|
|
dFd2[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
|
|
// We assume that the capillary pressure is independent of the polymer concentration.
|
|
// Hence, no more contributions.
|
|
}
|
|
|
|
// dFs/dm_1 \cdot dm_1/ds
|
|
*dFsds[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmds[0];
|
|
// dFc/dm_1 \cdot dm_1/ds
|
|
*dFcds[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmds[0];
|
|
|
|
// dFs/dm_2 \cdot dm_2/ds
|
|
*dFsds[ pix[1] ] -= dt * f1 / mt * v1 * dmds[1];
|
|
*dFsds[ pix[1] ] += dt * f1 * gflux * dmds[1];
|
|
// dFc/dm_2 \cdot dm_2/ds
|
|
*dFcds[ pix[1] ] -= dt * f1 / mt * v1 * mc * dmds[1];
|
|
*dFcds[ pix[1] ] += dt * f1 * gflux * mc * dmds[1];
|
|
|
|
// dFs/dm_1 \cdot dm_1/dc
|
|
*dFsdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * dmobwatdc;
|
|
// dFc/dm_1 \cdot dm_1/dc
|
|
*dFcdc[ pix[0] ] += dt * (1 - f1) / mt * v1 * mc * dmobwatdc;
|
|
*dFcdc[ pix[0] ] += dt * f1 * v1 * dmcdc; // Polymer is only carried by water.
|
|
}
|
|
|
|
template <class Grid>
|
|
void
|
|
accumulation(const Grid& g,
|
|
const int cell,
|
|
double* F, // Residual vector,
|
|
double* dF // Jacobian, same convention as for fluxConnection.
|
|
) const {
|
|
(void) g;
|
|
|
|
const double pv = store_.porevol(cell);
|
|
const double dps = fluid_.deadporespace();
|
|
const double rhor = fluid_.rockdensity();
|
|
const double poro = store_.porosity(cell);
|
|
|
|
F[0] += pv * store_.ds(cell);
|
|
F[1] += pv * (1 - dps) * store_.dsc(cell) + rhor*(1 - poro)/poro*pv*store_.dcads(cell);
|
|
dF[0] += pv;
|
|
dF[1] += 0.;
|
|
dF[2] += pv * (1 - dps) * store_.c(cell);
|
|
dF[3] += pv * (1 - dps) * store_.sw(cell) + rhor*(1 - poro)/poro*pv*store_.dcadsdc(cell);
|
|
}
|
|
|
|
template <class Grid ,
|
|
class SourceTerms>
|
|
void
|
|
sourceTerms(const Grid& g ,
|
|
const SourceTerms* src,
|
|
const int i ,
|
|
const double dt ,
|
|
double* J ,
|
|
double* F ) const {
|
|
|
|
(void) g;
|
|
|
|
double dflux = -src->flux[i]; // ->flux[] is rate of *inflow*
|
|
|
|
if (dflux < 0) {
|
|
// src -> cell, affects residual only.
|
|
*F += dt * dflux * src->saturation[2*i + 0];
|
|
} else {
|
|
// cell -> src
|
|
const int cell = src->cell[i];
|
|
const double* m = store_.mob (cell);
|
|
const double* dm = store_.dmobds(cell);
|
|
|
|
const double mt = m[0] + m[1];
|
|
|
|
assert (! ((m[0] < 0) || (m[1] < 0)));
|
|
assert (mt > 0);
|
|
|
|
const double f = m[0] / mt;
|
|
const double df = ((1 - f)*dm[0] - f*dm[1]) / mt;
|
|
|
|
*F += dt * dflux * f;
|
|
*J += dt * dflux * df;
|
|
}
|
|
}
|
|
template <class Grid>
|
|
void
|
|
initGravityTrans(const Grid& g ,
|
|
const std::vector<double> & htrans) {
|
|
|
|
assert (htrans.size() ==
|
|
static_cast<std::vector<double>::size_type>(g.cell_facepos[ g.number_of_cells ]));
|
|
|
|
for (int f = 0; f < g.number_of_faces; ++f) {
|
|
store_.trans(f) = 0.0;
|
|
}
|
|
|
|
for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
|
|
for (; i < g.cell_facepos[c + 1]; ++i) {
|
|
int f = g.cell_faces[i];
|
|
|
|
assert (htrans[i] > 0.0);
|
|
|
|
store_.trans(f) += 1.0 / htrans[i];
|
|
}
|
|
}
|
|
|
|
for (int f = 0; f < g.number_of_faces; ++f) {
|
|
store_.trans(f) = 1.0 / store_.trans(f);
|
|
}
|
|
|
|
if (gravity_) {
|
|
this->computeStaticGravity(g);
|
|
}
|
|
}
|
|
|
|
// -----------------------------------------------------------------
|
|
// Newton control
|
|
// -----------------------------------------------------------------
|
|
|
|
template <class ReservoirState,
|
|
class Grid ,
|
|
class JacobianSystem>
|
|
void
|
|
initStep(const ReservoirState& state,
|
|
const Grid& g ,
|
|
JacobianSystem& sys ) {
|
|
|
|
(void) state; // Suppress 'unused' warning.
|
|
|
|
typename JacobianSystem::vector_type& x =
|
|
sys.vector().writableSolution();
|
|
|
|
assert (x.size() == (::std::size_t) (2*g.number_of_cells));
|
|
|
|
if (init_step_use_previous_sol_) {
|
|
std::fill(x.begin(), x.end(), 0.0);
|
|
} else {
|
|
std::fill(x.begin(), x.end(), 0.0);
|
|
const std::vector<double>& s = state.saturation();
|
|
const std::vector<double>& c = state.concentration();
|
|
for (int cell = 0, ncell = g.number_of_cells; cell < ncell; ++cell) {
|
|
// Impose s=0.5 at next time level as an NR initial value.
|
|
x[2*cell + 0] = 0.5 - s[2*cell + 0];
|
|
x[2*cell + 1] = 1e-5 - c[cell];
|
|
}
|
|
}
|
|
}
|
|
|
|
template <class ReservoirState,
|
|
class Grid ,
|
|
class JacobianSystem>
|
|
bool
|
|
initIteration(const ReservoirState& state,
|
|
const Grid& g ,
|
|
JacobianSystem& sys) {
|
|
|
|
double s[2];
|
|
double mob[2];
|
|
double dmobds[4];
|
|
double dmobwatdc;
|
|
double c, cmax;
|
|
double mc, dmcdc;
|
|
double pc, dpc;
|
|
|
|
const typename JacobianSystem::vector_type& x =
|
|
sys.vector().solution();
|
|
const ::std::vector<double>& sat = state.saturation();
|
|
const ::std::vector<double>& cpoly = state.concentration();
|
|
const ::std::vector<double>& cmaxpoly = state.maxconcentration();
|
|
|
|
bool in_range = true;
|
|
for (int cell = 0; cell < g.number_of_cells; ++cell) {
|
|
// Store wat-sat, sat-change, cpoly, (sat * cpoly)-change for accumulation().
|
|
store_.ds(cell) = x[2*cell + 0];
|
|
s[0] = sat[cell*2 + 0] + x[2*cell + 0];
|
|
c = cpoly[cell] + x[2*cell + 1];
|
|
store_.sw(cell) = s[0];
|
|
store_.c(cell) = c;
|
|
cmax = std::max(c, cmaxpoly[cell]);
|
|
store_.cmax(cell) = cmax;
|
|
store_.dsc(cell) = s[0]*c - sat[cell*2 + 0]*cpoly[cell];
|
|
double dcadsdc;
|
|
double cads;
|
|
fluid_.adsorption(cpoly[cell], cmaxpoly[cell], cads, dcadsdc);
|
|
store_.dcads(cell) = -cads;
|
|
fluid_.adsorption(c, cmax, cads, dcadsdc);
|
|
store_.dcads(cell) += cads;
|
|
store_.dcadsdc(cell) = dcadsdc;
|
|
double s_min = fluid_.s_min(cell);
|
|
double s_max = fluid_.s_max(cell);
|
|
|
|
if ( s[0] < (s_min - sat_tol_) || s[0] > (s_max + sat_tol_) ) {
|
|
// if (s[0] < s_min){
|
|
// std::cout << "Warning: s out of range, s-s_min = " << s_min-s[0] << std::endl;
|
|
// }
|
|
// if (s[0] > s_max){
|
|
// std::cout << "Warning: s out of range, s-s_max = " << s[0]-s_max << std::endl;
|
|
// }
|
|
in_range = false; //line search fails
|
|
}
|
|
s[0] = std::max(s_min, s[0]);
|
|
s[0] = std::min(s_max, s[0]);
|
|
s[1] = 1 - s[0];
|
|
|
|
fluid_.mobility(cell, s, c, cmax, mob, dmobds, dmobwatdc);
|
|
fluid_.computeMc(c, mc, dmcdc);
|
|
fluid_.pc(cell, s, pc, dpc);
|
|
|
|
store_.mob (cell)[0] = mob [0];
|
|
store_.mob (cell)[1] = mob [1];
|
|
store_.dmobds(cell)[0] = dmobds[0*2 + 0];
|
|
store_.dmobds(cell)[1] = -dmobds[1*2 + 1];
|
|
store_.dmobwatdc(cell) = dmobwatdc;
|
|
store_.mc(cell) = mc;
|
|
store_.dmcdc(cell) = dmcdc;
|
|
store_.pc(cell) = pc;
|
|
store_.dpc(cell) = dpc;
|
|
}
|
|
if (!in_range) {
|
|
std::cout << "Warning: initIteration() - s was clamped in some cells.\n";
|
|
}
|
|
return in_range;
|
|
}
|
|
|
|
template <class ReservoirState,
|
|
class Grid ,
|
|
class NewtonIterate >
|
|
void
|
|
finishIteration(const ReservoirState& state,
|
|
const Grid& g ,
|
|
NewtonIterate& it ) {
|
|
// Nothing to do at end of iteration in this model.
|
|
(void) state; (void) g; (void) it;
|
|
typedef typename NewtonIterate::vector_type vector_t;
|
|
}
|
|
|
|
template <class Grid ,
|
|
class SolutionVector,
|
|
class ReservoirState>
|
|
void
|
|
finishStep(const Grid& g ,
|
|
const SolutionVector& x ,
|
|
ReservoirState& state) {
|
|
|
|
double *s = &state.saturation()[0*2 + 0];
|
|
double *c = &state.concentration()[0*1 + 0];
|
|
double *cmax = &state.maxconcentration()[0*1 + 0];
|
|
|
|
for (int cell = 0; cell < g.number_of_cells; ++cell, s += 2, c += 1, cmax +=1) {
|
|
s[0] += x[2*cell + 0];
|
|
c[0] += x[2*cell + 1];
|
|
cmax[0] = std::max(c[0], cmax[0]);
|
|
double s_min = fluid_.s_min(cell);
|
|
double s_max = fluid_.s_max(cell);
|
|
assert(s[0] >= s_min - sat_tol_);
|
|
assert(s[0] <= s_max + sat_tol_);
|
|
s[0] = std::max(s_min, s[0]);
|
|
s[0] = std::min(s_max, s[0]);
|
|
s[1] = 1.0 - s[0];
|
|
}
|
|
}
|
|
|
|
private:
|
|
void
|
|
upwindMobility(const double dflux,
|
|
const double gflux,
|
|
const int* n ,
|
|
int* pix ,
|
|
double* m ,
|
|
double* dmds ,
|
|
double& dmobwatdc ,
|
|
double& mc,
|
|
double& dmcdc) const {
|
|
bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) ||
|
|
( (! (dflux > 0)) && (! (gflux > 0)) );
|
|
|
|
if (equal_sign) {
|
|
|
|
if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
|
|
else { pix[0] = 1; }
|
|
|
|
m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
|
|
mc = store_.mc(n[ pix[0] ]);
|
|
|
|
if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
|
|
else { pix[1] = 1; }
|
|
|
|
m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
|
|
|
|
} else {
|
|
|
|
if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
|
|
else { pix[1] = 1; }
|
|
|
|
m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
|
|
|
|
if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
|
|
else { pix[0] = 1; }
|
|
|
|
m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
|
|
mc = store_.mc(n[ pix[0] ]);
|
|
}
|
|
|
|
dmds[0] = store_.dmobds(n[ pix[0] ]) [ 0 ];
|
|
dmds[1] = store_.dmobds(n[ pix[1] ]) [ 1 ];
|
|
dmobwatdc = store_.dmobwatdc(n[ pix[0] ]);
|
|
dmcdc = store_.dmcdc(n[ pix[0] ]);
|
|
}
|
|
|
|
template <class Grid>
|
|
void
|
|
computeStaticGravity(const Grid& g) {
|
|
|
|
const int d = g.dimensions;
|
|
|
|
for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
|
|
const double* cc = g.cell_centroids + (c * d);
|
|
|
|
for (; i < g.cell_facepos[c + 1]; ++i) {
|
|
const int f = g.cell_faces[i];
|
|
const double* fc = g.face_centroids + (f * d);
|
|
|
|
double dg = 0.0;
|
|
for (int j = 0; j < d; ++j) {
|
|
dg += gravity_[j] * (fc[j] - cc[j]);
|
|
}
|
|
|
|
store_.dg(i) = store_.trans(f) * dg;
|
|
}
|
|
}
|
|
}
|
|
|
|
double
|
|
gravityFlux(const int f) const {
|
|
double gflux;
|
|
|
|
if (gravity_) {
|
|
int i1 = f2hf_[2*f + 0];
|
|
int i2 = f2hf_[2*f + 1];
|
|
|
|
assert ((i1 >= 0) && (i2 >= 0));
|
|
|
|
gflux = store_.dg(i1) - store_.dg(i2);
|
|
gflux *= store_.drho();
|
|
} else {
|
|
gflux = 0.0;
|
|
}
|
|
|
|
return gflux;
|
|
}
|
|
void
|
|
capFlux(const int f,const int* n,double& pcflux, double* dpcflux) const {
|
|
//double capflux;
|
|
int i1 = n[0];
|
|
int i2 = n[1];
|
|
assert ((i1 >= 0) && (i2 >= 0));
|
|
//double sgn=-1.0;
|
|
pcflux = store_.trans(f)*(store_.pc(i2) - store_.pc(i1));
|
|
dpcflux[0] = -store_.trans(f)*store_.dpc(i1);
|
|
dpcflux[1] = store_.trans(f)*store_.dpc(i2);
|
|
}
|
|
|
|
TwophaseFluidPolymer fluid_ ;
|
|
const double* gravity_;
|
|
std::vector<int> f2hf_ ;
|
|
polymer_reorder::ModelParameterStorage store_ ;
|
|
bool init_step_use_previous_sol_;
|
|
double sat_tol_;
|
|
};
|
|
}
|
|
#endif /* OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER */
|