mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Computes s and c-residual and its Jacobian at a cell interface.
This commit is contained in:
Xavier Raynaud
parent
c888d995da
commit
54844f60c1
675
opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp
Normal file
675
opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp
Normal file
@ -0,0 +1,675 @@
|
||||
/*===========================================================================
|
||||
//
|
||||
// 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 spu_2p {
|
||||
class ModelParameterStorage {
|
||||
public:
|
||||
ModelParameterStorage(int nc, int totconn)
|
||||
: drho_(0.0), mob_(0), dmobds_(0), dmobwatdc_(0), mc_(0), dmcds_(0), dmcdc_(0),
|
||||
porevol_(0), deadporespace_(0), dg_(0), sw_(0), c_(0), ds_(0), dsc_(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; // dmcds_
|
||||
alloc_sz += nc; // dmcdc_
|
||||
alloc_sz += 1 * nc; // porevol_
|
||||
alloc_sz += 1 * nc; // deadporespave__
|
||||
alloc_sz += 1 * totconn; // dg_
|
||||
alloc_sz += 1 * nc; // sw_
|
||||
alloc_sz += 1 * nc; // c_
|
||||
alloc_sz += 1 * nc; // dc_
|
||||
alloc_sz += 1 * nc; // dsc_
|
||||
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_ = dmobwatds_ + (2 * nc );
|
||||
mc_ = dmobdc_ + (1 * nc );
|
||||
dmcds_ = mc_ + (1 * nc );
|
||||
dmcdc_ = dmcds_ + (1 * nc );
|
||||
porevol_ = dmcdc_ + (1 * nc );
|
||||
deadporespace_ = porevol_ + (1 * nc );
|
||||
dg_ = deadporespace_ + (1 * nc );
|
||||
sw_ = dg_ + (1 * totconn);
|
||||
c_ = sw_ + (1 * nc );
|
||||
ds_ = c_ + (1 * nc );
|
||||
dsc_ = ds_ + (1 * nc );
|
||||
pc_ = dsc_ + (1 * nc );
|
||||
dpc_ = pc_ + (1 * nc );
|
||||
trans_ = dpc_ + (1 * nc );
|
||||
}
|
||||
|
||||
double& drho () { return drho_ ; }
|
||||
double drho () const { return drho_ ; }
|
||||
|
||||
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 dmobdc_[cell]; }
|
||||
double dmobwatdc (int cell) const { return dmobdc_[cell]; }
|
||||
|
||||
double& mc (int cell) { return mc_[cell]; }
|
||||
double mc (int cell) const { return mc_[cell]; }
|
||||
|
||||
double& dmcds (int cell) { return dmcds_[cell]; }
|
||||
double dmcds (int cell) const { return dmcds_[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* deadporespace() { return deadporespace_ ; }
|
||||
double deadporespace(int cell) const { return deadporespace_[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& 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& 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 *mob_ ;
|
||||
double *dmobds_ ;
|
||||
double *dmobdc_ ;
|
||||
double *mc_ ;
|
||||
double *dmcds_ ;
|
||||
double *dmcdc_ ;
|
||||
double *porevol_ ;
|
||||
double *dg_ ;
|
||||
double *ds_ ;
|
||||
double *dsc_ ;
|
||||
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());
|
||||
}
|
||||
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, dmcds, dmcdc;
|
||||
upwindMobility(dflux, gflux, n, pix, m, dmds, dmobwatdc, mc, dmcds, dmcdc);
|
||||
|
||||
assert (! ((m[0] < 0) || (m[1] < 0)));
|
||||
|
||||
double mt = m[0] + m[1];
|
||||
assert (mt > 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] = dFsds2; dFsds[1] = dFsds1;
|
||||
dFsdc[0] = dFsdc2; dFsdc[1] = dFsdc1;
|
||||
dFcds[0] = dFcds2; dFcds[1] = dFcds1;
|
||||
dFcdc[0] = dFcdc2; dFcdc[1] = dFcdc1;
|
||||
// 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];
|
||||
*dFcds[ pix[0] ] += dt * f1 * v1 * dmcds; // Polymer is only carried by water.
|
||||
|
||||
// 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 * dmobwatdc; // 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 = store_.deadporespace(cell);
|
||||
|
||||
F[0] += pv * store_.ds(cell);
|
||||
F[1] += pv * (1 - dps) * store_.dsc(cell);
|
||||
dF[0] += pv;
|
||||
dF[1] += 0.;
|
||||
dF[2] += pv * (1 - dps) * store_.c(cell);
|
||||
dF[3] += pv * (1 - dps) * store_.sw(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 c = src->cell[i];
|
||||
const double* m = store_.mob (c);
|
||||
const double* dm = store_.dmob(c);
|
||||
|
||||
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) (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();
|
||||
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];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class ReservoirState,
|
||||
class Grid ,
|
||||
class JacobianSystem>
|
||||
bool
|
||||
initIteration(const ReservoirState& state,
|
||||
const Grid& g ,
|
||||
JacobianSystem& sys) {
|
||||
|
||||
double s[2], mob[2], dmobds[2 * 2], dmobdc[2 * 2];
|
||||
double c;
|
||||
double mc, dmcds, 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();
|
||||
|
||||
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;
|
||||
store_.dsc(cell) = s[0]*c - sat[cell*2 + 0]*cpoly[cell];
|
||||
|
||||
|
||||
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, mob, dmobds, dmobwatdc);
|
||||
fluid_.mc(cell, s, mc, dmcds, 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_.dmcds(cell) = dmcds;
|
||||
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];
|
||||
|
||||
for (int cell = 0; cell < g.number_of_cells; ++cell, s += 2) {
|
||||
s[0] += x[2*cell + 0];
|
||||
c[0] += x[2*cell + 1];
|
||||
c += 1;
|
||||
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& dmcds,
|
||||
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] ]);
|
||||
dmcds = store_.dmcds(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);
|
||||
}
|
||||
|
||||
TwophaseFluid fluid_ ;
|
||||
const double* gravity_;
|
||||
std::vector<int> f2hf_ ;
|
||||
spu_2p::ModelParameterStorage store_ ;
|
||||
bool init_step_use_previous_sol_;
|
||||
double sat_tol_;
|
||||
};
|
||||
}
|
||||
#endif /* OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER */
|
||||
Loading…
Reference in New Issue
Block a user