OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merged branches: profiling -> default.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-08-23 15:02:56 +02:00
commit 6ad3e6c050
7 changed files with 564 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
examples/Makefile.am
View File

@ -11,7 +11,8 @@ $(BOOST_SYSTEM_LIB)
noinst_PROGRAMS = \
polymer_reorder \
sim_poly2p_incomp_reorder
sim_poly2p_incomp_reorder \
test_singlecellsolves
polymer_reorder_SOURCES = \
polymer_reorder.cpp
@ -19,3 +20,6 @@ polymer_reorder.cpp
sim_poly2p_incomp_reorder_SOURCES = \
sim_poly2p_incomp_reorder.cpp
test_singlecellsolves_SOURCES = \
test_singlecellsolves.cpp

75
examples/polymer_reorder.cpp
View File

@ -43,6 +43,7 @@
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
//#include <opm/core/linalg/LinearSolverAGMG.hpp>
#include <opm/core/transport/transport_source.h>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
@ -107,6 +108,7 @@ static void outputState(const UnstructuredGrid& grid,
Opm::writeVtkData(grid, dm, vtkfile);
// Write data (not grid) in Matlab format
dm["faceflux"] = &state.faceflux();
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
@ -118,6 +120,7 @@ static void outputState(const UnstructuredGrid& grid,
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
#if PROFILING
std::ostringstream fname;
fname << output_dir << "/" << "residualcounts" << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
std::ofstream file(fname.str().c_str());
@ -132,6 +135,7 @@ static void outputState(const UnstructuredGrid& grid,
file << it->res_s << "," << it->cell << "," << std::setprecision(15) << it->s << "," << std::setprecision(15) << it->c << "\n";
}
file.close();
#endif
}
@ -395,7 +399,8 @@ main(int argc, char** argv)
const int num_transport_substeps = param.getDefault("num_transport_substeps", 1);
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
bool use_deck = param.getDefault("use_deck", true);
use_deck = param.has("deck_filename") && use_deck;
boost::scoped_ptr<Opm::GridManager> grid;
boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
boost::scoped_ptr<Opm::WellsManager> wells;
@ -456,8 +461,6 @@ main(int argc, char** argv)
wells.reset(new Opm::WellsManager());
// Timer init.
simtimer.init(param);
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
@ -481,34 +484,39 @@ main(int argc, char** argv)
}
// Init polymer properties.
// Setting defaults to provide a simple example case.
double c_max = param.getDefault("c_max_limit", 5.0);
double mix_param = param.getDefault("mix_param", 1.0);
double rock_density = param.getDefault("rock_density", 1000.0);
double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15);
double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability
double c_max_ads = param.getDefault("c_max_ads", 1.);
int ads_index = param.getDefault<int>("ads_index", Opm::PolymerProperties::NoDesorption);
std::vector<double> c_vals_visc(2, -1e100);
c_vals_visc[0] = 0.0;
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);
c_vals_ads[0] = 0.0;
c_vals_ads[1] = 2.0;
c_vals_ads[2] = 8.0;
// Here we set up adsorption equal to zero.
std::vector<double> ads_vals(3, -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);
bool use_deck_fluid = param.getDefault("use_deck_fluid", false);
if(!use_deck_fluid){
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(param));
double c_max = param.getDefault("c_max_limit", 5.0);
double mix_param = param.getDefault("mix_param", 1.0);
double rock_density = param.getDefault("rock_density", 1000.0);
double dead_pore_vol = param.getDefault("dead_pore_vol", 0.1);
double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability
double c_max_ads = param.getDefault("c_max_ads", 1.);
int ads_index = param.getDefault<int>("ads_index", Opm::PolymerProperties::NoDesorption);
std::vector<double> c_vals_visc(2, -1e100);
c_vals_visc[0] = 0.0;
c_vals_visc[1] = c_max;
std::vector<double> visc_mult_vals(2, -1e100);
visc_mult_vals[0] = 1.0;
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] = 8.0;
// Here we set up adsorption equal to zero.
std::vector<double> ads_vals(2, -1e100);
ads_vals[0] = 0.0;
ads_vals[1] = 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);
}else{
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename);
rock_comp.reset(new Opm::RockCompressibility(deck));
polyprop.readFromDeck(deck);
}
}
// Initialize polymer inflow function.
@ -590,6 +598,7 @@ main(int argc, char** argv)
// Solvers init.
// Linear solver.
Opm::LinearSolverFactory linsolver(param);
//Opm::LinearSolverAGMG linsolver;
// Pressure solver.
const double *grav = use_gravity ? &gravity[0] : 0;
Opm::IncompTpfaPolymer psolver(*grid->c_grid(), *props, rock_comp.get(), polyprop, linsolver,
@ -607,6 +616,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);
}

45
examples/sim_poly2p_incomp_reorder.cpp
View File

@ -47,6 +47,7 @@
#include <opm/polymer/PolymerProperties.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
#include <algorithm>
#include <iostream>
@ -54,6 +55,18 @@
#include <numeric>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
// ----------------- Main program -----------------
@ -205,17 +218,20 @@ main(int argc, char** argv)
// Linear solver.
LinearSolverFactory linsolver(param);
// Warn if any parameters are unused.
// if (param.anyUnused()) {
// std::cout << "-------------------- Unused parameters: --------------------\n";
// param.displayUsage();
// std::cout << "----------------------------------------------------------------" << std::endl;
// }
// Write parameters used for later reference.
// if (output) {
// param.writeParam(output_dir + "/spu_2p.param");
// }
bool output = param.getDefault("output", true);
if (output) {
std::string output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
param.writeParam(output_dir + "/simulation.param");
}
std::cout << "\n\n================ Starting main simulation loop ===============\n"
@ -237,6 +253,7 @@ main(int argc, char** argv)
grav);
SimulatorTimer simtimer;
simtimer.init(param);
warnIfUnusedParams(param);
WellState well_state;
well_state.init(0, state);
rep = simulator.run(simtimer, state, well_state);
@ -270,7 +287,7 @@ main(int argc, char** argv)
<< "\n (number of steps: "
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
// Create new wells, well_satate
// Create new wells, well_state
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
@ -290,8 +307,10 @@ main(int argc, char** argv)
bcs.c_bcs(),
linsolver,
grav);
SimulatorReport epoch_rep
= simulator.run(simtimer, state, well_state);
if (epoch == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
// Update total timing report and remember step number.
rep += epoch_rep;

248
examples/test_singlecellsolves.cpp Normal file
View File

@ -0,0 +1,248 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/polymer/TransportModelPolymer.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <boost/scoped_ptr.hpp>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
using namespace Opm;
// std::cout << "\n================ Test program for single-cell solves with polymer ===============\n\n";
parameter::ParameterGroup param(argc, argv, false);
param.disableOutput();
// std::cout << "--------------- Reading parameters ---------------" << std::endl;
// If we have a "deck_filename", grid and props will be read from that.
boost::scoped_ptr<EclipseGridParser> deck;
boost::scoped_ptr<GridManager> grid;
boost::scoped_ptr<IncompPropertiesInterface> props;
PolymerState state;
Opm::PolymerProperties poly_props;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
// -------- Initialising section ----------
// Grid init.
grid.reset(new GridManager(2, 1, 1, 1.0, 1.0, 1.0));
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, 0.0, state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
double smin[2], smax[2];
props->satRange(1, &cell, smin, smax);
if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
state.concentration()[cell] = poly_init;
state.maxconcentration()[cell] = poly_init;
} else {
state.saturation()[2*cell + 0] = 0.;
state.saturation()[2*cell + 1] = 1.;
state.concentration()[cell] = 0.;
state.maxconcentration()[cell] = 0.;
}
}
}
// Init polymer properties.
// Setting defaults to provide a simple example case.
double c_max = param.getDefault("c_max_limit", 5.0);
double mix_param = param.getDefault("mix_param", 1.0);
double rock_density = param.getDefault("rock_density", 1000.0);
double dead_pore_vol = param.getDefault("dead_pore_vol", 0.0); // Note that we default to no dps here!
double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability
double c_max_ads = param.getDefault("c_max_ads", 1.);
int ads_index = param.getDefault<int>("ads_index", Opm::PolymerProperties::NoDesorption);
std::vector<double> c_vals_visc(2, -1e100);
c_vals_visc[0] = 0.0;
c_vals_visc[1] = 7.0;
std::vector<double> visc_mult_vals(2, -1e100);
visc_mult_vals[0] = 1.0;
// poly_props.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
visc_mult_vals[1] = 20.0;
std::vector<double> c_vals_ads(3, -1e100);
c_vals_ads[0] = 0.0;
c_vals_ads[1] = 2.0;
c_vals_ads[2] = 8.0;
std::vector<double> ads_vals(3, -1e100);
ads_vals[0] = 0.0;
ads_vals[1] = 0.0015;
ads_vals[2] = 0.0025;
// ads_vals[1] = 0.0;
// ads_vals[2] = 0.0;
poly_props.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);
// Initialising src
int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> src(num_cells, 0.0);
// Compute pore volumes, in order to enable specifying injection rate
// terms of total pore volume.
std::vector<double> porevol;
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
const double default_injection = 1.0;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_sec", default_injection)
*porevol[0];
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
// Boundary conditions.
FlowBCManager bcs;
// Linear solver.
LinearSolverFactory linsolver(param);
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportModelPolymer::SingleCellMethod method;
std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
if (method_string == "Bracketing") {
method = Opm::TransportModelPolymer::Bracketing;
} else if (method_string == "Newton") {
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);
}
Opm::TransportModelPolymer reorder_model(*grid->c_grid(), *props, poly_props,
method, nl_tolerance, nl_maxiter);
// Warn if any parameters are unused.
// if (param.anyUnused()) {
// std::cout << "-------------------- Unused parameters: --------------------\n";
// param.displayUsage();
// std::cout << "----------------------------------------------------------------" << std::endl;
// }
// Write parameters to file for later reference.
param.writeParam("test_singlecellsolves.param");
// Setting up a number of input (s, c) pairs and solving.
// HACK warning: we manipulate the source term,
// but the compressibility term in the solver
// assumes that all inflow is water inflow. Therefore
// one must zero the compressibility term in
// TransportModelPolymer line 365 before compiling this program.
// (To fix this we should add proper all-phase src terms.)
std::vector<double> transport_src = src;
const double dt = param.getDefault("dt", 1.0);
const int num_sats = 501;
const int num_concs = 501;
// Find the face between cell 0 and 1...
const UnstructuredGrid& ug = *grid->c_grid();
int face01 = -1;
for (int f = 0; f < ug.number_of_faces; ++f) {
if (ug.face_cells[2*f] == 0 && ug.face_cells[2*f+1] == 1) {
face01 = f;
break;
}
}
if (face01 == -1) {
THROW("Could not find face adjacent to cells [0 1]");
}
state.faceflux()[face01] = src[0];
for (int sats = 0; sats < num_sats; ++sats) {
const double s = double(sats)/double(num_sats - 1);
const double ff = s; // Simplified a lot...
for (int conc = 0; conc < num_concs; ++conc) {
const double c = poly_props.cMax()*double(conc)/double(num_concs - 1);
// std::cout << "(s, c) = (" << s << ", " << c << ")\n";
transport_src[0] = src[0]*ff;
// Resetting the state for next run.
state.saturation()[0] = 0.0;
state.saturation()[1] = 0.0;
state.concentration()[0] = 0.0;
state.concentration()[1] = 0.0;
state.maxconcentration()[0] = 0.0;
state.maxconcentration()[1] = 0.0;
reorder_model.solve(&state.faceflux()[0],
&porevol[0],
&transport_src[0],
dt,
c,
state.saturation(),
state.concentration(),
state.maxconcentration());
#if PROFILING
// Extract residual counts.
typedef std::list<Opm::TransportModelPolymer::Newton_Iter> ListRes;
const ListRes& res_counts = reorder_model.res_counts;
double counts[2] = { 0, 0 };
for (ListRes::const_iterator it = res_counts.begin(); it != res_counts.end(); ++it) {
if (it->cell == 0) {
++counts[it->res_s];
}
}
// std::cout << "c residual count: " << counts[0] << '\n';
// std::cout << "s residual count: " << counts[1] << '\n';
std::cout << counts[0] << ' ' << counts[1] << ' ' << s << ' ' << c << '\n';
#endif
}
}
}

6
opm/polymer/PolymerProperties.cpp
View File

@ -278,8 +278,10 @@ namespace Opm
dmobwat_dc = eff_relperm_wat*dinv_mu_w_eff_dc
+ deff_relperm_wat_dc*inv_mu_w_eff;
dmob_ds[0*2 + 0] = deff_relperm_wat_ds*inv_mu_w_eff;
dmob_ds[0*2 + 1] = (drelperm_ds[0*2 + 1] - drelperm_ds[1*2 + 2])/visc[1];
dmob_ds[1*2 + 0] = -deff_relperm_wat_ds*inv_mu_w_eff;
// one have to deside which variables to derive
// here the full derivative is written out
dmob_ds[0*2 + 1] = 0.0*(drelperm_ds[0*2 + 1] - drelperm_ds[1*2 + 1])/visc[1];
dmob_ds[1*2 + 0] = -0.0*deff_relperm_wat_ds*inv_mu_w_eff;
dmob_ds[1*2 + 1] = (drelperm_ds[1*2 + 1] - drelperm_ds[0*2 + 1])/visc[1];
}
}

245
opm/polymer/TransportModelPolymer.cpp
View File

@ -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.
@ -193,14 +196,17 @@ namespace Opm
cmax_(0),
fractionalflow_(grid.number_of_cells, -1.0),
mc_(grid.number_of_cells, -1.0),
method_(method)
method_(method),
adhoc_safety_(1.0)
{
if (props.numPhases() != 2) {
THROW("Property object must have 2 phases");
}
visc_ = props.viscosity();
#if PROFILING
res_counts.clear();
#endif
// Set up smin_ and smax_
int num_cells = props.numCells();
@ -241,8 +247,10 @@ namespace Opm
toWaterSat(saturation, saturation_);
concentration_ = &concentration[0];
cmax_ = &cmax[0];
#if PROFILING
res_counts.clear();
reorderAndTransport(grid_, darcyflux);
#endif
reorderAndTransport(grid_, darcyflux);
toBothSat(saturation_, saturation);
}
@ -288,7 +296,9 @@ namespace Opm
}
double operator()(double s) const
{
#if PROFILING
tm_.res_counts.push_back(Newton_Iter(true, cell_, s, c_));
#endif
double ff;
tm_.fracFlow(s, c_, cmax0_, cell_, ff);
return s - s0_ + dtpv_*(outflux_*ff + influx_ + s*comp_term_);
@ -333,7 +343,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];
@ -375,8 +385,10 @@ namespace Opm
+ rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0)
+ dtpv*(outflux*ff*mc + influx_polymer)
+ dtpv*(s_arg*c_arg*(1.0 - dps) - rhor*c_ads)*comp_term;
#if PROFILING
tm.res_counts.push_back(Newton_Iter(true, cell, s_arg, c_arg));
tm.res_counts.push_back(Newton_Iter(false, cell, s_arg, c_arg));
#endif
}
@ -399,8 +411,10 @@ namespace Opm
tm.polyprops_.adsorption(c0, cmax0, c_ads0);
double c_ads;
tm.polyprops_.adsorption(c, cmax0, c_ads);
#if PROFILING
tm.res_counts.push_back(Newton_Iter(false, cell, s, c));
double res = (1 - dps)*s*c - (1 - dps)*s0*c0
#endif
double res = (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0)
+ dtpv*(outflux*ff*mc + influx_polymer)
+ dtpv*(s*c*(1.0 - dps) - rhor*c_ads)*comp_term;
@ -559,14 +573,18 @@ namespace Opm
}
if (if_res_s) {
res[0] = s - s0 + dtpv*(outflux*ff + influx + s*comp_term);
#if PROFILING
tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1]));
#endif
}
if (if_res_c) {
res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
+ dtpv*(outflux*ff*mc + influx_polymer)
+ dtpv*(s*c*(1.0 - dps) - rhor*ads)*comp_term;
#if PROFILING
tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1]));
#endif
}
if (if_dres_s_dsdc) {
dres_s_dsdc[0] = 1 + dtpv*(outflux*dff_dsdc[0] + comp_term);
@ -586,7 +604,9 @@ namespace Opm
tm.fracFlow(s, c, cmax0, cell, ff);
if (if_res_s) {
res[0] = s - s0 + dtpv*(outflux*ff + influx + s*comp_term);
#if PROFILING
tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1]));
#endif
}
if (if_res_c) {
tm.computeMc(c, mc);
@ -596,7 +616,9 @@ namespace Opm
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
+ dtpv*(outflux*ff*mc + influx_polymer)
+ dtpv*(s*c*(1.0 - dps) - rhor*ads)*comp_term;
#if PROFILING
tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1]));
#endif
}
}
@ -641,6 +663,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_);
}
@ -651,7 +679,7 @@ namespace Opm
{
ResidualC res(*this, cell);
const double a = 0.0;
const double b = polyprops_.cMax()*1.1; // Add 10% to account for possible non-monotonicity of hyperbolic system.
const double b = polyprops_.cMax()*adhoc_safety_; // Add 10% to account for possible non-monotonicity of hyperbolic system.
int iters_used;
// Check if current state is an acceptable solution.
@ -700,7 +728,7 @@ namespace Opm
}
double x_min[2] = { 0.0, 0.0 };
double x_max[2] = { 1.0, polyprops_.cMax()*1.1 };
double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ };
double x_min_res_s[2] = { x_min[0], x_min[1] };
double x_max_res_s[2] = { x_max[0], x_min[0] };
double x_min_res_sc[2] = { x_min[0], x_min[1] };
@ -857,8 +885,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()*adhoc_safety_ };
bool successfull_newton_step = true;
// initialize x_new to avoid warning
double x_new[2] = {0.0, 0.0};
@ -932,6 +960,154 @@ 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];
if(x[1]> polyprops_.cMax()){
x[1]= polyprops_.cMax()/2.0;
}
if(x[1]<0){
x[1]=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()*adhoc_safety_ };
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=tol_;
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("NewtonSimple 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)
{
@ -1022,11 +1198,14 @@ namespace Opm
double dmobwat_dc;
polyprops_.effectiveMobilitiesBoth(c, cmax, visc_, relperm, drelperm_ds,
mob, dmob_ds, dmobwat_dc, if_with_der);
ff = mob[0]/(mob[0] + mob[1]);
if (if_with_der) {
dmob_dc[0] = dmobwat_dc;
dmob_dc[1] = 0.;
dff_dsdc[0] = (dmob_ds[0]*mob[1] + dmob_ds[3]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to s
//dff_dsdc[0] = (dmob_ds[0]*mob[1] + dmob_ds[3]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to s
// at the moment the dmob_ds only have diagonal elements since the saturation is derivated out in effectiveMobilitiesBouth
dff_dsdc[0] = ((dmob_ds[0]-dmob_ds[2])*mob[1] - (dmob_ds[1]-dmob_ds[3])*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to s
dff_dsdc[1] = (dmob_dc[0]*mob[1] - dmob_dc[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to c
}
}
@ -1221,7 +1400,7 @@ namespace Opm
}
const double a = 0.0;
const double b = polyprops_.cMax()*1.1; // Add 10% to account for possible non-monotonicity of hyperbolic system.
const double b = polyprops_.cMax()*adhoc_safety_; // Add 10% to account for possible non-monotonicity of hyperbolic system.
int iters_used;
concentration_[cell] = RootFinder::solve(res_c, a, b, maxit_, tol_, iters_used);
ResidualSGrav res_s(res_c);
@ -1293,11 +1472,11 @@ namespace Opm
void TransportModelPolymer::solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation,
std::vector<double>& concentration,
std::vector<double>& cmax)
const double* porevolume,
const double dt,
std::vector<double>& saturation,
std::vector<double>& concentration,
std::vector<double>& cmax)
{
// initialize variables.
porevolume_ = porevolume;
@ -1471,12 +1650,12 @@ namespace
return res_eq_.computeResidualC(x_c);
}
bool solveNewtonStep(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq,
const double* res, double* x_new) {
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 +1681,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

11
opm/polymer/TransportModelPolymer.hpp
View File

@ -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);
@ -113,8 +114,9 @@ namespace Opm
const int pos,
const double* gravflux);
int solveGravityColumn(const std::vector<int>& cells);
void scToc(const double* x, double* x_c) const;
// for testing
#if PROFILING
class Newton_Iter {
public:
bool res_s;
@ -131,8 +133,8 @@ namespace Opm
};
std::list<Newton_Iter> res_counts;
#endif
void scToc(const double* x, double* x_c) const;
private:
const UnstructuredGrid& grid_;
@ -156,7 +158,8 @@ namespace Opm
std::vector<double> mc_; // one per cell
const double* visc_;
SingleCellMethod method_;
double adhoc_safety_;
// For gravity segregation.
std::vector<double> gravflux_;
std::vector<double> mob_;