mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-15 02:53:26 -06:00
Made new single cell solvers.
This commit is contained in:
Halvor M. Nilsen
commit
b13416be86
@ -495,19 +495,14 @@ main(int argc, char** argv)
|
||||
c_vals_visc[1] = 7.0;
|
||||
std::vector<double> visc_mult_vals(2, -1e100);
|
||||
visc_mult_vals[0] = 1.0;
|
||||
// polyprop.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
|
||||
visc_mult_vals[1] = 1.0;
|
||||
std::vector<double> c_vals_ads(3, -1e100);
|
||||
visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
|
||||
std::vector<double> c_vals_ads(2, -1e100);
|
||||
c_vals_ads[0] = 0.0;
|
||||
c_vals_ads[1] = 2.0;
|
||||
c_vals_ads[2] = 8.0;
|
||||
c_vals_ads[1] = 8.0;
|
||||
// Here we set up adsorption equal to zero.
|
||||
std::vector<double> ads_vals(3, -1e100);
|
||||
std::vector<double> ads_vals(2, -1e100);
|
||||
ads_vals[0] = 0.0;
|
||||
ads_vals[1] = 0.0;
|
||||
ads_vals[2] = 0.0;
|
||||
// ads_vals[1] = 0.0;
|
||||
// ads_vals[2] = 0.0;
|
||||
polyprop.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads,
|
||||
static_cast<Opm::PolymerProperties::AdsorptionBehaviour>(ads_index),
|
||||
c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals);
|
||||
@ -609,6 +604,10 @@ main(int argc, char** argv)
|
||||
method = Opm::TransportModelPolymer::Newton;
|
||||
} else if (method_string == "Gradient") {
|
||||
method = Opm::TransportModelPolymer::Gradient;
|
||||
} else if (method_string == "NewtonSimpleSC") {
|
||||
method = Opm::TransportModelPolymer::NewtonSimpleSC;
|
||||
} else if (method_string == "NewtonSimpleC") {
|
||||
method = Opm::TransportModelPolymer::NewtonSimpleC;
|
||||
} else {
|
||||
THROW("Unknown method: " << method_string);
|
||||
}
|
||||
@ -767,6 +766,11 @@ main(int argc, char** argv)
|
||||
}
|
||||
} while (!well_control_passed);
|
||||
|
||||
// Update pore volumes if rock is compressible.
|
||||
if (rock_comp->isActive()) {
|
||||
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
|
||||
}
|
||||
|
||||
// Process transport sources (to include bdy terms and well flows).
|
||||
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
|
||||
wells->c_wells(), well_state.perfRates(), reorder_src);
|
||||
|
||||
@ -296,6 +296,11 @@ namespace Opm
|
||||
ptime += pt;
|
||||
} while (false);
|
||||
|
||||
// Update pore volumes if rock is compressible.
|
||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
||||
}
|
||||
|
||||
// Process transport sources (to include bdy terms and well flows).
|
||||
Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
|
||||
wells_, well_state.perfRates(), transport_src);
|
||||
|
||||
@ -26,7 +26,7 @@
|
||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
||||
#include <cmath>
|
||||
#include <list>
|
||||
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
// Choose error policy for scalar solves here.
|
||||
typedef Opm::RegulaFalsi<Opm::WarnAndContinueOnError> RootFinder;
|
||||
|
||||
@ -114,8 +114,11 @@ namespace
|
||||
return std::max(std::abs(res[0]), std::abs(res[1]));
|
||||
}
|
||||
|
||||
bool solveNewtonStep(const double* , const Opm::TransportModelPolymer::ResidualEquation&,
|
||||
const double*, double*);
|
||||
bool solveNewtonStepSC(const double* , const Opm::TransportModelPolymer::ResidualEquation&,
|
||||
const double*, double*);
|
||||
bool solveNewtonStepC(const double* , const Opm::TransportModelPolymer::ResidualEquation&,
|
||||
const double*, double*);
|
||||
|
||||
|
||||
|
||||
// Define a piecewise linear curve along which we will look for zero of the "s" or "r" residual.
|
||||
@ -333,7 +336,7 @@ namespace Opm
|
||||
comp_term = tm.source_[cell]; // Note: this assumes that all source flux is water.
|
||||
dtpv = tm.dt_/tm.porevolume_[cell];
|
||||
porosity = tm.porosity_[cell];
|
||||
s = -1e100;
|
||||
s = s0;
|
||||
|
||||
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
|
||||
int f = tm.grid_.cell_faces[i];
|
||||
@ -388,7 +391,7 @@ namespace Opm
|
||||
// Solve for s first.
|
||||
// s = modifiedRegulaFalsi(res_s, std::max(tm.smin_[2*cell], dps), tm.smax_[2*cell],
|
||||
// tm.maxit_, tm.tol_, iters_used);
|
||||
s = RootFinder::solve(res_s, s0, 0.0, 1.0,
|
||||
s = RootFinder::solve(res_s, s, 0.0, 1.0,
|
||||
tm.maxit_, tm.tol_, iters_used);
|
||||
double ff;
|
||||
tm.fracFlow(s, c, cmax0, cell, ff);
|
||||
@ -641,6 +644,12 @@ namespace Opm
|
||||
case Gradient:
|
||||
solveSingleCellGradient(cell);
|
||||
break;
|
||||
case NewtonSimpleSC:
|
||||
solveSingleCellNewtonSimple(cell,true);
|
||||
break;
|
||||
case NewtonSimpleC:
|
||||
solveSingleCellNewtonSimple(cell,false);
|
||||
break;
|
||||
default:
|
||||
THROW("Unknown method " << method_);
|
||||
}
|
||||
@ -857,8 +866,8 @@ namespace Opm
|
||||
// res_eq.computeResidual(x, res, mc, ff);
|
||||
}
|
||||
|
||||
double x_min[2] = { 0.0, 0.0 };
|
||||
double x_max[2] = { 1.0, polyprops_.cMax()*1.1 };
|
||||
const double x_min[2] = { 0.0, 0.0 };
|
||||
const double x_max[2] = { 1.0, polyprops_.cMax()*1.1 };
|
||||
bool successfull_newton_step = true;
|
||||
// initialize x_new to avoid warning
|
||||
double x_new[2] = {0.0, 0.0};
|
||||
@ -932,6 +941,147 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
|
||||
void TransportModelPolymer::solveSingleCellNewtonSimple(int cell,bool use_sc)
|
||||
{
|
||||
const int max_iters_split = maxit_;
|
||||
int iters_used_split = 0;
|
||||
|
||||
// Check if current state is an acceptable solution.
|
||||
ResidualEquation res_eq(*this, cell);
|
||||
double x[2] = {saturation_[cell], concentration_[cell]};
|
||||
double res[2];
|
||||
double mc;
|
||||
double ff;
|
||||
res_eq.computeResidual(x, res, mc, ff);
|
||||
if (norm(res) <= tol_) {
|
||||
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
|
||||
fractionalflow_[cell] = ff;
|
||||
mc_[cell] = mc;
|
||||
return;
|
||||
}else{
|
||||
x[0] = saturation_[cell]-res[0];
|
||||
if((x[0]>1) || (x[0]<0)){
|
||||
x[0] = 0.5;
|
||||
x[1] = x[1];
|
||||
}
|
||||
if(x[0]>0){
|
||||
x[1] = concentration_[cell]*saturation_[cell]-res[1];
|
||||
x[1] = x[1]/x[0];
|
||||
}else{
|
||||
x[1]=0;
|
||||
}
|
||||
//x[0]=0.5;x[1]=polyprops_.cMax()/2.0;
|
||||
res_eq.computeResidual(x, res, mc, ff);
|
||||
|
||||
}
|
||||
|
||||
|
||||
// double x_min[2] = { std::max(polyprops_.deadPoreVol(), smin_[2*cell]), 0.0 };
|
||||
double x_min[2] = { 0.0, 0.0 };
|
||||
double x_max[2] = { 1.0, polyprops_.cMax()*1.1 };
|
||||
bool successfull_newton_step = true;
|
||||
double x_new[2];
|
||||
double res_new[2];
|
||||
ResSOnCurve res_s_on_curve(res_eq);
|
||||
ResCOnCurve res_c_on_curve(res_eq);
|
||||
|
||||
while ((norm(res) > tol_) &&
|
||||
(iters_used_split < max_iters_split) &&
|
||||
successfull_newton_step) {
|
||||
// We first try a Newton step
|
||||
double dres_s_dsdc[2];
|
||||
double dres_c_dsdc[2];
|
||||
double dx=1e-6;
|
||||
double tmp_x[2];
|
||||
if(!(x[0]>0)){
|
||||
tmp_x[0]=dx;
|
||||
tmp_x[1]=0;
|
||||
}else{
|
||||
tmp_x[0]=x[0];
|
||||
tmp_x[1]=x[1];
|
||||
}
|
||||
res_eq.computeJacobiRes(tmp_x, dres_s_dsdc, dres_c_dsdc);
|
||||
double dFx_dx,dFx_dy,dFy_dx,dFy_dy;
|
||||
double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
|
||||
if(use_sc){
|
||||
dFx_dx=(dres_s_dsdc[0]-tmp_x[1]*dres_s_dsdc[1]/tmp_x[0]);
|
||||
dFx_dy= (dres_s_dsdc[1]/tmp_x[0]);
|
||||
dFy_dx=(dres_c_dsdc[0]-tmp_x[1]*dres_c_dsdc[1]/tmp_x[0]);
|
||||
dFy_dy= (dres_c_dsdc[1]/tmp_x[0]);
|
||||
}else{
|
||||
dFx_dx= dres_s_dsdc[0];
|
||||
dFx_dy= dres_s_dsdc[1];
|
||||
dFy_dx= dres_c_dsdc[0];
|
||||
dFy_dy= dres_c_dsdc[1];
|
||||
}
|
||||
det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
|
||||
if(det==0){
|
||||
THROW("det is zero");
|
||||
}
|
||||
|
||||
double alpha=1;
|
||||
int max_lin_it=10;
|
||||
int lin_it=0;
|
||||
/*
|
||||
std::cout << "Nonlinear " << iters_used_split
|
||||
<< " " << norm(res_new)
|
||||
<< " " << norm(res) << std::endl;*/
|
||||
|
||||
std::cout << "x" << std::endl;
|
||||
std::cout << " " << x[0] << " " << x[1] << std::endl;
|
||||
std::cout << "dF" << std::endl;
|
||||
std::cout << " " << dFx_dx << " " << dFx_dy << std::endl;
|
||||
std::cout << " " << dFy_dx << " " << dFy_dy << std::endl;
|
||||
std::cout << "dres" << std::endl;
|
||||
std::cout << " " << dres_s_dsdc[0] << " " << dres_s_dsdc[1] << std::endl;
|
||||
std::cout << " " << dres_c_dsdc[0] << " " << dres_c_dsdc[1] << std::endl;
|
||||
std::cout << std::endl;
|
||||
|
||||
res_new[0]=res[0]*2;
|
||||
res_new[1]=res[1]*2;
|
||||
double update[2]={(res[0]*dFy_dy - res[1]*dFx_dy)/det,
|
||||
(res[1]*dFx_dx - res[0]*dFy_dx)/det};
|
||||
while((norm(res_new)>norm(res)) && (lin_it<max_lin_it)) {
|
||||
x_new[0] = x[0] - alpha*update[0];
|
||||
if(use_sc){
|
||||
x_new[1] = x[1]*x[0] - alpha*update[1];
|
||||
x_new[1] = (x_new[0]>0) ? x_new[1]/x_new[0] : 0.0;
|
||||
}else{
|
||||
x_new[1] = x[1] - alpha*update[1];
|
||||
}
|
||||
check_interval(x_min, x_max, x_new);
|
||||
res_eq.computeResidual(x_new, res_new, mc, ff);
|
||||
// std::cout << " " << res_new[0] << " " << res_new[1] << std::endl;
|
||||
alpha=alpha/2.0;
|
||||
lin_it=lin_it+1;
|
||||
// std::cout << "Linear iterations" << lin_it << " " << norm(res_new) << std::endl;
|
||||
}
|
||||
if (lin_it>=max_lin_it) {
|
||||
successfull_newton_step = false;
|
||||
} else {
|
||||
x[0] = x_new[0];
|
||||
x[1] = x_new[1];
|
||||
res[0] = res_new[0];
|
||||
res[1] = res_new[1];
|
||||
iters_used_split += 1;
|
||||
successfull_newton_step = true;;
|
||||
}
|
||||
// std::cout << "Nonlinear " << iters_used_split << " " << norm(res) << std::endl;
|
||||
}
|
||||
|
||||
if ((iters_used_split >= max_iters_split) || (norm(res) > tol_)) {
|
||||
MESSAGE("Newton for single cell did not work in cell number " << cell);
|
||||
solveSingleCellBracketing(cell);
|
||||
} else {
|
||||
concentration_[cell] = x[1];
|
||||
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
|
||||
saturation_[cell] = x[0];
|
||||
fractionalflow_[cell] = ff;
|
||||
mc_[cell] = mc;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
void TransportModelPolymer::solveMultiCell(const int num_cells, const int* cells)
|
||||
{
|
||||
@ -1471,12 +1621,12 @@ namespace
|
||||
return res_eq_.computeResidualC(x_c);
|
||||
}
|
||||
|
||||
bool solveNewtonStep(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq,
|
||||
bool solveNewtonStepSC(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq,
|
||||
const double* res, double* x_new) {
|
||||
|
||||
double dres_s_dsdc[2];
|
||||
double dres_c_dsdc[2];
|
||||
double dx=1e-11;
|
||||
double dx=1e-8;
|
||||
double x[2];
|
||||
if(!(x[0]>0)){
|
||||
x[0] = dx;
|
||||
@ -1502,6 +1652,30 @@ namespace
|
||||
return true;
|
||||
}
|
||||
}
|
||||
bool solveNewtonStepC(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq,
|
||||
const double* res, double* x_new) {
|
||||
|
||||
double dres_s_dsdc[2];
|
||||
double dres_c_dsdc[2];
|
||||
double dx=1e-8;
|
||||
double x[2];
|
||||
if(!(x[0]>0)){
|
||||
x[0] = dx;
|
||||
x[1] = 0;
|
||||
} else {
|
||||
x[0] = xx[0];
|
||||
x[1] = xx[1];
|
||||
}
|
||||
res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
|
||||
double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1];
|
||||
if (std::abs(det) < 1e-8) {
|
||||
return false;
|
||||
} else {
|
||||
x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det;
|
||||
x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} // Anonymous namespace
|
||||
|
||||
@ -40,7 +40,7 @@ namespace Opm
|
||||
{
|
||||
public:
|
||||
|
||||
enum SingleCellMethod { Bracketing, Newton, Gradient };
|
||||
enum SingleCellMethod { Bracketing, Newton, Gradient, NewtonSimpleSC, NewtonSimpleC};
|
||||
enum GradientMethod { Analytic, FinDif }; // Analytic is chosen (hard-coded)
|
||||
|
||||
/// Construct solver.
|
||||
@ -106,6 +106,7 @@ namespace Opm
|
||||
void solveSingleCellBracketing(int cell);
|
||||
void solveSingleCellNewton(int cell);
|
||||
void solveSingleCellGradient(int cell);
|
||||
void solveSingleCellNewtonSimple(int cell,bool use_sc);
|
||||
class ResidualEquation;
|
||||
|
||||
void initGravity(const double* grav);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user