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

Merge pull request #17 from atgeirr/combined

Browse Source 
Adapt to changes in OPM/opm-core#203
This commit is contained in:
Bård Skaflestad 2013-03-18 12:14:46 -07:00
commit 8200b5dff7
14 changed files with 202 additions and 202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
examples/sim_poly2p_comp_reorder.cpp
View File

@ -25,19 +25,19 @@
#include <opm/core/pressure/FlowBCManager.hpp> #include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/GridManager.hpp> #include <opm/core/grid/GridManager.hpp>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp> #include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp> #include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp> #include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp> #include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp> #include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/BlackoilPropertiesBasic.hpp> #include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp> #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp> #include <opm/core/linalg/LinearSolverFactory.hpp>

12
examples/sim_poly2p_incomp_reorder.cpp
View File

@ -25,19 +25,19 @@
#include <opm/core/pressure/FlowBCManager.hpp> #include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/GridManager.hpp> #include <opm/core/grid/GridManager.hpp>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp> #include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp> #include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp> #include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp> #include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp> #include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp> #include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp> #include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp> #include <opm/core/linalg/LinearSolverFactory.hpp>

32
examples/test_singlecellsolves.cpp
View File

@ -25,26 +25,26 @@
#include <opm/core/pressure/FlowBCManager.hpp> #include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/GridManager.hpp> #include <opm/core/grid/GridManager.hpp>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp> #include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp> #include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp> #include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp> #include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp> #include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp> #include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp> #include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp> #include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/polymer/PolymerState.hpp> #include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp> #include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/IncompTpfaPolymer.hpp> #include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/polymer/TransportModelPolymer.hpp> #include <opm/polymer/TransportSolverTwophasePolymer.hpp>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <boost/scoped_ptr.hpp> #include <boost/scoped_ptr.hpp>
@ -154,22 +154,22 @@ main(int argc, char** argv)
// Reordering solver. // Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30); const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportModelPolymer::SingleCellMethod method; Opm::TransportSolverTwophasePolymer::SingleCellMethod method;
std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing")); std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
if (method_string == "Bracketing") { if (method_string == "Bracketing") {
method = Opm::TransportModelPolymer::Bracketing; method = Opm::TransportSolverTwophasePolymer::Bracketing;
} else if (method_string == "Newton") { } else if (method_string == "Newton") {
method = Opm::TransportModelPolymer::Newton; method = Opm::TransportSolverTwophasePolymer::Newton;
} else if (method_string == "Gradient") { } else if (method_string == "Gradient") {
method = Opm::TransportModelPolymer::Gradient; method = Opm::TransportSolverTwophasePolymer::Gradient;
} else if (method_string == "NewtonSimpleSC") { } else if (method_string == "NewtonSimpleSC") {
method = Opm::TransportModelPolymer::NewtonSimpleSC; method = Opm::TransportSolverTwophasePolymer::NewtonSimpleSC;
} else if (method_string == "NewtonSimpleC") { } else if (method_string == "NewtonSimpleC") {
method = Opm::TransportModelPolymer::NewtonSimpleC; method = Opm::TransportSolverTwophasePolymer::NewtonSimpleC;
} else { } else {
THROW("Unknown method: " << method_string); THROW("Unknown method: " << method_string);
} }
Opm::TransportModelPolymer reorder_model(*grid->c_grid(), *props, poly_props, Opm::TransportSolverTwophasePolymer reorder_model(*grid->c_grid(), *props, poly_props,
method, nl_tolerance, nl_maxiter); method, nl_tolerance, nl_maxiter);
// Warn if any parameters are unused. // Warn if any parameters are unused.
@ -187,7 +187,7 @@ main(int argc, char** argv)
// but the compressibility term in the solver // but the compressibility term in the solver
// assumes that all inflow is water inflow. Therefore // assumes that all inflow is water inflow. Therefore
// one must zero the compressibility term in // one must zero the compressibility term in
// TransportModelPolymer line 365 before compiling this program. // TransportSolverTwophasePolymer line 365 before compiling this program.
// (To fix this we should add proper all-phase src terms.) // (To fix this we should add proper all-phase src terms.)
std::vector<double> transport_src = src; std::vector<double> transport_src = src;
const double dt = param.getDefault("dt", 1.0); const double dt = param.getDefault("dt", 1.0);
@ -232,7 +232,7 @@ main(int argc, char** argv)
#ifdef PROFILING #ifdef PROFILING
// Extract residual counts. // Extract residual counts.
typedef std::list<Opm::TransportModelPolymer::Newton_Iter> ListRes; typedef std::list<Opm::TransportSolverTwophasePolymer::Newton_Iter> ListRes;
const ListRes& res_counts = reorder_model.res_counts; const ListRes& res_counts = reorder_model.res_counts;
double counts[2] = { 0, 0 }; double counts[2] = { 0, 0 };
for (ListRes::const_iterator it = res_counts.begin(); it != res_counts.end(); ++it) { for (ListRes::const_iterator it = res_counts.begin(); it != res_counts.end(); ++it) {

6
opm/polymer/CompressibleTpfaPolymer.cpp
View File

@ -20,8 +20,8 @@
#include <opm/polymer/PolymerBlackoilState.hpp> #include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/CompressibleTpfaPolymer.hpp> #include <opm/polymer/CompressibleTpfaPolymer.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp> #include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/pressure/tpfa/ifs_tpfa.h> #include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h> #include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/mimetic/mimetic.h> #include <opm/core/pressure/mimetic/mimetic.h>
@ -32,7 +32,7 @@
#include <opm/core/simulator/WellState.hpp> #include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp> #include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <iomanip> #include <iomanip>
#include <cmath> #include <cmath>
#include <algorithm> #include <algorithm>

2
opm/polymer/IncompPropertiesDefaultPolymer.hpp
View File

@ -20,7 +20,7 @@
#ifndef OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED #ifndef OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED
#define OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED #define OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED
#include <opm/core/fluid/IncompPropertiesBasic.hpp> #include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp> #include <opm/core/utility/linearInterpolation.hpp>

6
opm/polymer/IncompTpfaPolymer.cpp
View File

@ -20,8 +20,8 @@
#include <opm/polymer/IncompTpfaPolymer.hpp> #include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp> #include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/pressure/tpfa/ifs_tpfa.h> #include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h> #include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/mimetic/mimetic.h> #include <opm/core/pressure/mimetic/mimetic.h>
@ -33,7 +33,7 @@
#include <opm/core/simulator/WellState.hpp> #include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp> #include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <iomanip> #include <iomanip>
#include <cmath> #include <cmath>
#include <algorithm> #include <algorithm>

2
opm/polymer/PolymerInflow.cpp
View File

@ -19,7 +19,7 @@
#include <opm/polymer/PolymerInflow.hpp> #include <opm/polymer/PolymerInflow.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp> #include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <map> #include <map>
namespace Opm namespace Opm

20
opm/polymer/SimulatorCompressiblePolymer.cpp
View File

@ -28,7 +28,7 @@
#include <opm/polymer/CompressibleTpfaPolymer.hpp> #include <opm/polymer/CompressibleTpfaPolymer.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h> #include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp> #include <opm/core/simulator/SimulatorReport.hpp>
@ -40,14 +40,14 @@
#include <opm/core/wells/WellsManager.hpp> #include <opm/core/wells/WellsManager.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp> #include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/utility/ColumnExtract.hpp> #include <opm/core/grid/ColumnExtract.hpp>
#include <opm/core/utility/Units.hpp> #include <opm/core/utility/Units.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp> #include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp> #include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/TransportModelCompressiblePolymer.hpp> #include <opm/polymer/TransportSolverTwophaseCompressiblePolymer.hpp>
#include <opm/polymer/PolymerInflow.hpp> #include <opm/polymer/PolymerInflow.hpp>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/polymerUtilities.hpp> #include <opm/polymer/polymerUtilities.hpp>
@ -130,7 +130,7 @@ namespace Opm
const double* gravity_; const double* gravity_;
// Solvers // Solvers
CompressibleTpfaPolymer psolver_; CompressibleTpfaPolymer psolver_;
TransportModelCompressiblePolymer tsolver_; TransportSolverTwophaseCompressiblePolymer tsolver_;
// Needed by column-based gravity segregation solver. // Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_; std::vector< std::vector<int> > columns_;
// Misc. data // Misc. data
@ -198,7 +198,7 @@ namespace Opm
param.getDefault("nl_pressure_maxiter", 10), param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells()), gravity, wells_manager.c_wells()),
tsolver_(grid, props, poly_props, tsolver_(grid, props, poly_props,
TransportModelCompressiblePolymer::Bracketing, TransportSolverTwophaseCompressiblePolymer::Bracketing,
param.getDefault("nl_tolerance", 1e-9), param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30)) param.getDefault("nl_maxiter", 30))
{ {
@ -223,12 +223,12 @@ namespace Opm
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Transport related init. // Transport related init.
TransportModelCompressiblePolymer::SingleCellMethod method; TransportSolverTwophaseCompressiblePolymer::SingleCellMethod method;
std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing")); std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
if (method_string == "Bracketing") { if (method_string == "Bracketing") {
method = Opm::TransportModelCompressiblePolymer::Bracketing; method = Opm::TransportSolverTwophaseCompressiblePolymer::Bracketing;
} else if (method_string == "Newton") { } else if (method_string == "Newton") {
method = Opm::TransportModelCompressiblePolymer::Newton; method = Opm::TransportSolverTwophaseCompressiblePolymer::Newton;
} else { } else {
THROW("Unknown method: " << method_string); THROW("Unknown method: " << method_string);
} }

20
opm/polymer/SimulatorPolymer.cpp
View File

@ -29,7 +29,7 @@
#include <opm/polymer/IncompTpfaPolymer.hpp> #include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/newwells.h> #include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h> #include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp> #include <opm/core/simulator/SimulatorReport.hpp>
@ -40,14 +40,14 @@
#include <opm/core/wells/WellsManager.hpp> #include <opm/core/wells/WellsManager.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp> #include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/utility/ColumnExtract.hpp> #include <opm/core/grid/ColumnExtract.hpp>
#include <opm/core/utility/Units.hpp> #include <opm/core/utility/Units.hpp>
#include <opm/polymer/PolymerState.hpp> #include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp> #include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/TransportModelPolymer.hpp> #include <opm/polymer/TransportSolverTwophasePolymer.hpp>
#include <opm/polymer/PolymerInflow.hpp> #include <opm/polymer/PolymerInflow.hpp>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/polymerUtilities.hpp> #include <opm/polymer/polymerUtilities.hpp>
@ -141,7 +141,7 @@ namespace Opm
const double* gravity_; const double* gravity_;
// Solvers // Solvers
IncompTpfaPolymer psolver_; IncompTpfaPolymer psolver_;
TransportModelPolymer tsolver_; TransportSolverTwophasePolymer tsolver_;
// Needed by column-based gravity segregation solver. // Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_; std::vector< std::vector<int> > columns_;
// Misc. data // Misc. data
@ -208,7 +208,7 @@ namespace Opm
param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10), param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs), gravity, wells_manager.c_wells(), src, bcs),
tsolver_(grid, props, poly_props, TransportModelPolymer::Bracketing, tsolver_(grid, props, poly_props, TransportSolverTwophasePolymer::Bracketing,
param.getDefault("nl_tolerance", 1e-9), param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30)) param.getDefault("nl_maxiter", 30))
{ {
@ -239,12 +239,12 @@ namespace Opm
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Transport related init. // Transport related init.
TransportModelPolymer::SingleCellMethod method; TransportSolverTwophasePolymer::SingleCellMethod method;
std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing")); std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
if (method_string == "Bracketing") { if (method_string == "Bracketing") {
method = Opm::TransportModelPolymer::Bracketing; method = Opm::TransportSolverTwophasePolymer::Bracketing;
} else if (method_string == "Newton") { } else if (method_string == "Newton") {
method = Opm::TransportModelPolymer::Newton; method = Opm::TransportSolverTwophasePolymer::Newton;
} else { } else {
THROW("Unknown method: " << method_string); THROW("Unknown method: " << method_string);
} }

124
opm/polymer/TransportModelCompressiblePolymer.cpp → opm/polymer/TransportSolverTwophaseCompressiblePolymer.cpp
View File

@ -18,8 +18,8 @@
*/ */
#include <opm/polymer/TransportModelCompressiblePolymer.hpp> #include <opm/polymer/TransportSolverTwophaseCompressiblePolymer.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp> #include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/transport/reorder/reordersequence.h> #include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/utility/RootFinders.hpp> #include <opm/core/utility/RootFinders.hpp>
@ -33,7 +33,7 @@
typedef Opm::RegulaFalsi<Opm::WarnAndContinueOnError> RootFinder; typedef Opm::RegulaFalsi<Opm::WarnAndContinueOnError> RootFinder;
class Opm::TransportModelCompressiblePolymer::ResidualEquation class Opm::TransportSolverTwophaseCompressiblePolymer::ResidualEquation
{ {
public: public:
GradientMethod gradient_method; GradientMethod gradient_method;
@ -55,9 +55,9 @@ public:
double rhor; double rhor;
double ads0; double ads0;
TransportModelCompressiblePolymer& tm; TransportSolverTwophaseCompressiblePolymer& tm;
ResidualEquation(TransportModelCompressiblePolymer& tmodel, int cell_index); ResidualEquation(TransportSolverTwophaseCompressiblePolymer& tmodel, int cell_index);
void computeResidual(const double* x, double* res) const; void computeResidual(const double* x, double* res) const;
void computeResidual(const double* x, double* res, double& mc, double& ff) const; void computeResidual(const double* x, double* res, double& mc, double& ff) const;
double computeResidualS(const double* x) const; double computeResidualS(const double* x) const;
@ -75,9 +75,9 @@ private:
double* dres_c_dsdc, double& mc, double& ff) const; double* dres_c_dsdc, double& mc, double& ff) const;
}; };
class Opm::TransportModelCompressiblePolymer::ResidualCGrav { class Opm::TransportSolverTwophaseCompressiblePolymer::ResidualCGrav {
public: public:
const TransportModelCompressiblePolymer& tm; const TransportSolverTwophaseCompressiblePolymer& tm;
const int cell; const int cell;
const double s0; const double s0;
const double c0; const double c0;
@ -92,7 +92,7 @@ public:
int nbcell[2]; int nbcell[2];
mutable double last_s; mutable double last_s;
ResidualCGrav(const TransportModelCompressiblePolymer& tmodel, ResidualCGrav(const TransportSolverTwophaseCompressiblePolymer& tmodel,
const std::vector<int>& cells, const std::vector<int>& cells,
const int pos, const int pos,
const double* gravflux); const double* gravflux);
@ -103,7 +103,7 @@ public:
double lastSaturation() const; double lastSaturation() const;
}; };
class Opm::TransportModelCompressiblePolymer::ResidualSGrav { class Opm::TransportSolverTwophaseCompressiblePolymer::ResidualSGrav {
public: public:
const ResidualCGrav& res_c_eq_; const ResidualCGrav& res_c_eq_;
double c; double c;
@ -151,7 +151,7 @@ namespace
namespace Opm namespace Opm
{ {
TransportModelCompressiblePolymer::TransportModelCompressiblePolymer(const UnstructuredGrid& grid, TransportSolverTwophaseCompressiblePolymer::TransportSolverTwophaseCompressiblePolymer(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props, const BlackoilPropertiesInterface& props,
const PolymerProperties& polyprops, const PolymerProperties& polyprops,
const SingleCellMethod method, const SingleCellMethod method,
@ -202,7 +202,7 @@ namespace Opm
void TransportModelCompressiblePolymer::setPreferredMethod(SingleCellMethod method) void TransportSolverTwophaseCompressiblePolymer::setPreferredMethod(SingleCellMethod method)
{ {
method_ = method; method_ = method;
} }
@ -210,7 +210,7 @@ namespace Opm
void TransportModelCompressiblePolymer::solve(const double* darcyflux, void TransportSolverTwophaseCompressiblePolymer::solve(const double* darcyflux,
const std::vector<double>& initial_pressure, const std::vector<double>& initial_pressure,
const std::vector<double>& pressure, const std::vector<double>& pressure,
const double* porevolume0, const double* porevolume0,
@ -243,7 +243,7 @@ namespace Opm
// Check immiscibility requirement (only done for first cell). // Check immiscibility requirement (only done for first cell).
if (A_[1] != 0.0 || A_[2] != 0.0) { if (A_[1] != 0.0 || A_[2] != 0.0) {
THROW("TransportCompressibleModelCompressibleTwophase requires a property object without miscibility."); THROW("TransportCompressibleSolverTwophaseCompressibleTwophase requires a property object without miscibility.");
} }
std::vector<int> seq(grid_.number_of_cells); std::vector<int> seq(grid_.number_of_cells);
std::vector<int> comp(grid_.number_of_cells + 1); std::vector<int> comp(grid_.number_of_cells + 1);
@ -273,11 +273,11 @@ namespace Opm
// //
// where influx is water influx, outflux is total outflux. // where influx is water influx, outflux is total outflux.
// Influxes are negative, outfluxes positive. // Influxes are negative, outfluxes positive.
struct TransportModelCompressiblePolymer::ResidualS struct TransportSolverTwophaseCompressiblePolymer::ResidualS
{ {
TransportModelCompressiblePolymer::ResidualEquation& res_eq_; TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq_;
const double c_; const double c_;
explicit ResidualS(TransportModelCompressiblePolymer::ResidualEquation& res_eq, explicit ResidualS(TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq,
const double c) const double c)
: res_eq_(res_eq), : res_eq_(res_eq),
c_(c) c_(c)
@ -298,11 +298,11 @@ namespace Opm
// \TODO doc me // \TODO doc me
// where ... // where ...
// Influxes are negative, outfluxes positive. // Influxes are negative, outfluxes positive.
struct TransportModelCompressiblePolymer::ResidualC struct TransportSolverTwophaseCompressiblePolymer::ResidualC
{ {
mutable double s; // Mutable in order to change it with every operator() call to be the last computed s value. mutable double s; // Mutable in order to change it with every operator() call to be the last computed s value.
TransportModelCompressiblePolymer::ResidualEquation& res_eq_; TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq_;
explicit ResidualC(TransportModelCompressiblePolymer::ResidualEquation& res_eq) explicit ResidualC(TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq) : res_eq_(res_eq)
{} {}
@ -346,7 +346,7 @@ namespace Opm
// ResidualEquation gathers parameters to construct the residual, computes its // ResidualEquation gathers parameters to construct the residual, computes its
// value and the values of its derivatives. // value and the values of its derivatives.
TransportModelCompressiblePolymer::ResidualEquation::ResidualEquation(TransportModelCompressiblePolymer& tmodel, int cell_index) TransportSolverTwophaseCompressiblePolymer::ResidualEquation::ResidualEquation(TransportSolverTwophaseCompressiblePolymer& tmodel, int cell_index)
: tm(tmodel) : tm(tmodel)
{ {
gradient_method = Analytic; gradient_method = Analytic;
@ -399,7 +399,7 @@ namespace Opm
} }
void TransportModelCompressiblePolymer::ResidualEquation::computeResidual(const double* x, double* res) const void TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeResidual(const double* x, double* res) const
{ {
double dres_s_dsdc[2]; double dres_s_dsdc[2];
double dres_c_dsdc[2]; double dres_c_dsdc[2];
@ -408,7 +408,7 @@ namespace Opm
computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff); computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
} }
void TransportModelCompressiblePolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const void TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const
{ {
double dres_s_dsdc[2]; double dres_s_dsdc[2];
double dres_c_dsdc[2]; double dres_c_dsdc[2];
@ -416,7 +416,7 @@ namespace Opm
} }
double TransportModelCompressiblePolymer::ResidualEquation::computeResidualS(const double* x) const double TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeResidualS(const double* x) const
{ {
double res[2]; double res[2];
double dres_s_dsdc[2]; double dres_s_dsdc[2];
@ -427,7 +427,7 @@ namespace Opm
return res[0]; return res[0];
} }
double TransportModelCompressiblePolymer::ResidualEquation::computeResidualC(const double* x) const double TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeResidualC(const double* x) const
{ {
double res[2]; double res[2];
double dres_s_dsdc[2]; double dres_s_dsdc[2];
@ -438,7 +438,7 @@ namespace Opm
return res[1]; return res[1];
} }
void TransportModelCompressiblePolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const void TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const
// If gradient_method == FinDif, use finite difference // If gradient_method == FinDif, use finite difference
// If gradient_method == Analytic, use analytic expresions // If gradient_method == Analytic, use analytic expresions
{ {
@ -448,7 +448,7 @@ namespace Opm
computeResAndJacobi(x, true, true, true, false, res, gradient, dres_c_dsdc, mc, ff); computeResAndJacobi(x, true, true, true, false, res, gradient, dres_c_dsdc, mc, ff);
} }
void TransportModelCompressiblePolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const void TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const
// If gradient_method == FinDif, use finite difference // If gradient_method == FinDif, use finite difference
// If gradient_method == Analytic, use analytic expresions // If gradient_method == Analytic, use analytic expresions
{ {
@ -459,7 +459,7 @@ namespace Opm
} }
// Compute the Jacobian of the residual equations. // Compute the Jacobian of the residual equations.
void TransportModelCompressiblePolymer::ResidualEquation::computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const void TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const
{ {
double res[2]; double res[2];
double mc; double mc;
@ -467,7 +467,7 @@ namespace Opm
computeResAndJacobi(x, false, false, true, true, res, dres_s_dsdc, dres_c_dsdc, mc, ff); computeResAndJacobi(x, false, false, true, true, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
} }
void TransportModelCompressiblePolymer::ResidualEquation::computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c, void TransportSolverTwophaseCompressiblePolymer::ResidualEquation::computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c,
const bool if_dres_s_dsdc, const bool if_dres_c_dsdc, const bool if_dres_s_dsdc, const bool if_dres_c_dsdc,
double* res, double* dres_s_dsdc, double* res, double* dres_s_dsdc,
double* dres_c_dsdc, double& mc, double& ff) const double* dres_c_dsdc, double& mc, double& ff) const
@ -567,34 +567,34 @@ namespace Opm
// Compute the "s" residual along the curve "curve" for a given residual equation "res_eq". // Compute the "s" residual along the curve "curve" for a given residual equation "res_eq".
// The operator() is sent to a root solver. // The operator() is sent to a root solver.
class TransportModelCompressiblePolymer::ResSOnCurve class TransportSolverTwophaseCompressiblePolymer::ResSOnCurve
{ {
public: public:
ResSOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq); ResSOnCurve(const TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq);
double operator()(const double t) const; double operator()(const double t) const;
CurveInSCPlane curve; CurveInSCPlane curve;
private: private:
const TransportModelCompressiblePolymer::ResidualEquation& res_eq_; const TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq_;
}; };
// Compute the "c" residual along the curve "curve" for a given residual equation "res_eq". // Compute the "c" residual along the curve "curve" for a given residual equation "res_eq".
// The operator() is sent to a root solver. // The operator() is sent to a root solver.
class TransportModelCompressiblePolymer::ResCOnCurve class TransportSolverTwophaseCompressiblePolymer::ResCOnCurve
{ {
public: public:
ResCOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq); ResCOnCurve(const TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq);
double operator()(const double t) const; double operator()(const double t) const;
CurveInSCPlane curve; CurveInSCPlane curve;
private: private:
const TransportModelCompressiblePolymer::ResidualEquation& res_eq_; const TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq_;
}; };
TransportModelCompressiblePolymer::ResSOnCurve::ResSOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq) TransportSolverTwophaseCompressiblePolymer::ResSOnCurve::ResSOnCurve(const TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq) : res_eq_(res_eq)
{ {
} }
double TransportModelCompressiblePolymer::ResSOnCurve::operator()(const double t) const double TransportSolverTwophaseCompressiblePolymer::ResSOnCurve::operator()(const double t) const
{ {
double x_of_t[2]; double x_of_t[2];
double x_c[2]; double x_c[2];
@ -603,12 +603,12 @@ namespace Opm
return res_eq_.computeResidualS(x_c); return res_eq_.computeResidualS(x_c);
} }
TransportModelCompressiblePolymer::ResCOnCurve::ResCOnCurve(const TransportModelCompressiblePolymer::ResidualEquation& res_eq) TransportSolverTwophaseCompressiblePolymer::ResCOnCurve::ResCOnCurve(const TransportSolverTwophaseCompressiblePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq) : res_eq_(res_eq)
{ {
} }
double TransportModelCompressiblePolymer::ResCOnCurve::operator()(const double t) const double TransportSolverTwophaseCompressiblePolymer::ResCOnCurve::operator()(const double t) const
{ {
double x_of_t[2]; double x_of_t[2];
double x_c[2]; double x_c[2];
@ -619,7 +619,7 @@ namespace Opm
void TransportModelCompressiblePolymer::solveSingleCell(const int cell) void TransportSolverTwophaseCompressiblePolymer::solveSingleCell(const int cell)
{ {
switch (method_) { switch (method_) {
case Bracketing: case Bracketing:
@ -640,7 +640,7 @@ namespace Opm
} }
void TransportModelCompressiblePolymer::solveSingleCellBracketing(int cell) void TransportSolverTwophaseCompressiblePolymer::solveSingleCellBracketing(int cell)
{ {
ResidualEquation res_eq(*this, cell); ResidualEquation res_eq(*this, cell);
@ -672,7 +672,7 @@ namespace Opm
// Newton method, where we first try a Newton step. Then, if it does not work well, we look for // Newton method, where we first try a Newton step. Then, if it does not work well, we look for
// the zero of either the residual in s or the residual in c along a specified piecewise linear // the zero of either the residual in s or the residual in c along a specified piecewise linear
// curve. In these cases, we can use a robust 1d solver. // curve. In these cases, we can use a robust 1d solver.
void TransportModelCompressiblePolymer::solveSingleCellGradient(int cell) void TransportSolverTwophaseCompressiblePolymer::solveSingleCellGradient(int cell)
{ {
int iters_used_falsi = 0; int iters_used_falsi = 0;
const int max_iters_split = maxit_; const int max_iters_split = maxit_;
@ -829,7 +829,7 @@ namespace Opm
} }
} }
void TransportModelCompressiblePolymer::solveSingleCellNewton(int cell, bool use_sc, void TransportSolverTwophaseCompressiblePolymer::solveSingleCellNewton(int cell, bool use_sc,
bool use_explicit_step) bool use_explicit_step)
{ {
const int max_iters_split = maxit_; const int max_iters_split = maxit_;
@ -975,7 +975,7 @@ namespace Opm
} }
void TransportModelCompressiblePolymer::solveMultiCell(const int num_cells, const int* cells) void TransportSolverTwophaseCompressiblePolymer::solveMultiCell(const int num_cells, const int* cells)
{ {
double max_s_change = 0.0; double max_s_change = 0.0;
double max_c_change = 0.0; double max_c_change = 0.0;
@ -1032,21 +1032,21 @@ namespace Opm
// << num_iters << " iterations." << std::endl; // << num_iters << " iterations." << std::endl;
} }
void TransportModelCompressiblePolymer::fracFlow(double s, double c, double cmax, void TransportSolverTwophaseCompressiblePolymer::fracFlow(double s, double c, double cmax,
int cell, double& ff) const int cell, double& ff) const
{ {
double dummy[2]; double dummy[2];
fracFlowBoth(s, c, cmax, cell, ff, dummy, false); fracFlowBoth(s, c, cmax, cell, ff, dummy, false);
} }
void TransportModelCompressiblePolymer::fracFlowWithDer(double s, double c, double cmax, void TransportSolverTwophaseCompressiblePolymer::fracFlowWithDer(double s, double c, double cmax,
int cell, double& ff, int cell, double& ff,
double* dff_dsdc) const double* dff_dsdc) const
{ {
fracFlowBoth(s, c, cmax, cell, ff, dff_dsdc, true); fracFlowBoth(s, c, cmax, cell, ff, dff_dsdc, true);
} }
void TransportModelCompressiblePolymer::fracFlowBoth(double s, double c, double cmax, int cell, void TransportSolverTwophaseCompressiblePolymer::fracFlowBoth(double s, double c, double cmax, int cell,
double& ff, double* dff_dsdc, double& ff, double* dff_dsdc,
bool if_with_der) const bool if_with_der) const
{ {
@ -1077,12 +1077,12 @@ namespace Opm
} }
} }
void TransportModelCompressiblePolymer::computeMc(double c, double& mc) const void TransportSolverTwophaseCompressiblePolymer::computeMc(double c, double& mc) const
{ {
polyprops_.computeMc(c, mc); polyprops_.computeMc(c, mc);
} }
void TransportModelCompressiblePolymer::computeMcWithDer(double c, double& mc, void TransportSolverTwophaseCompressiblePolymer::computeMcWithDer(double c, double& mc,
double &dmc_dc) const double &dmc_dc) const
{ {
polyprops_.computeMcWithDer(c, mc, dmc_dc); polyprops_.computeMcWithDer(c, mc, dmc_dc);
@ -1090,14 +1090,14 @@ namespace Opm
TransportModelCompressiblePolymer::ResidualSGrav::ResidualSGrav(const ResidualCGrav& res_c_eq, TransportSolverTwophaseCompressiblePolymer::ResidualSGrav::ResidualSGrav(const ResidualCGrav& res_c_eq,
const double c_init) const double c_init)
: res_c_eq_(res_c_eq), : res_c_eq_(res_c_eq),
c(c_init) c(c_init)
{ {
} }
double TransportModelCompressiblePolymer::ResidualSGrav::operator()(double s) const double TransportSolverTwophaseCompressiblePolymer::ResidualSGrav::operator()(double s) const
{ {
return res_c_eq_.computeGravResidualS(s, c); return res_c_eq_.computeGravResidualS(s, c);
} }
@ -1113,7 +1113,7 @@ namespace Opm
// where ... // where ...
// Influxes are negative, outfluxes positive. // Influxes are negative, outfluxes positive.
TransportModelCompressiblePolymer::ResidualCGrav::ResidualCGrav(const TransportModelCompressiblePolymer& tmodel, TransportSolverTwophaseCompressiblePolymer::ResidualCGrav::ResidualCGrav(const TransportSolverTwophaseCompressiblePolymer& tmodel,
const std::vector<int>& cells, const std::vector<int>& cells,
const int pos, const int pos,
const double* gravflux) // Always oriented towards next in column. Size = colsize - 1. const double* gravflux) // Always oriented towards next in column. Size = colsize - 1.
@ -1146,7 +1146,7 @@ namespace Opm
tm.polyprops_.adsorption(c0, cmax0, c_ads0); tm.polyprops_.adsorption(c0, cmax0, c_ads0);
} }
double TransportModelCompressiblePolymer::ResidualCGrav::operator()(double c) const double TransportSolverTwophaseCompressiblePolymer::ResidualCGrav::operator()(double c) const
{ {
ResidualSGrav res_s(*this); ResidualSGrav res_s(*this);
@ -1159,7 +1159,7 @@ namespace Opm
} }
double TransportModelCompressiblePolymer::ResidualCGrav::computeGravResidualS(double s, double c) const double TransportSolverTwophaseCompressiblePolymer::ResidualCGrav::computeGravResidualS(double s, double c) const
{ {
double mobcell[2]; double mobcell[2];
@ -1185,7 +1185,7 @@ namespace Opm
return res; return res;
} }
double TransportModelCompressiblePolymer::ResidualCGrav::computeGravResidualC(double s, double c) const double TransportSolverTwophaseCompressiblePolymer::ResidualCGrav::computeGravResidualC(double s, double c) const
{ {
double mobcell[2]; double mobcell[2];
@ -1216,13 +1216,13 @@ namespace Opm
return res; return res;
} }
double TransportModelCompressiblePolymer::ResidualCGrav::lastSaturation() const double TransportSolverTwophaseCompressiblePolymer::ResidualCGrav::lastSaturation() const
{ {
return last_s; return last_s;
} }
void TransportModelCompressiblePolymer::mobility(double s, double c, int cell, double* mob) const void TransportSolverTwophaseCompressiblePolymer::mobility(double s, double c, int cell, double* mob) const
{ {
double sat[2] = { s, 1.0 - s }; double sat[2] = { s, 1.0 - s };
double relperm[2]; double relperm[2];
@ -1231,7 +1231,7 @@ namespace Opm
polyprops_.effectiveMobilities(c, cmax0_[cell], &visc_[np*cell], relperm, mob); polyprops_.effectiveMobilities(c, cmax0_[cell], &visc_[np*cell], relperm, mob);
} }
void TransportModelCompressiblePolymer::initGravity(const double* grav) void TransportSolverTwophaseCompressiblePolymer::initGravity(const double* grav)
{ {
// Set up transmissibilities. // Set up transmissibilities.
std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]); std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
@ -1245,7 +1245,7 @@ namespace Opm
gravity_ = grav; gravity_ = grav;
} }
void TransportModelCompressiblePolymer::initGravityDynamic() void TransportSolverTwophaseCompressiblePolymer::initGravityDynamic()
{ {
// Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f) // Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f)
// + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ] // + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ]
@ -1276,7 +1276,7 @@ namespace Opm
} }
void TransportModelCompressiblePolymer::solveSingleCellGravity(const std::vector<int>& cells, void TransportSolverTwophaseCompressiblePolymer::solveSingleCellGravity(const std::vector<int>& cells,
const int pos, const int pos,
const double* gravflux) const double* gravflux)
{ {
@ -1303,7 +1303,7 @@ namespace Opm
mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]); mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]);
} }
int TransportModelCompressiblePolymer::solveGravityColumn(const std::vector<int>& cells) int TransportSolverTwophaseCompressiblePolymer::solveGravityColumn(const std::vector<int>& cells)
{ {
// Set up column gravflux. // Set up column gravflux.
const int nc = cells.size(); const int nc = cells.size();
@ -1363,7 +1363,7 @@ namespace Opm
} }
void TransportModelCompressiblePolymer::solveGravity(const std::vector<std::vector<int> >& columns, void TransportSolverTwophaseCompressiblePolymer::solveGravity(const std::vector<std::vector<int> >& columns,
const double dt, const double dt,
std::vector<double>& saturation, std::vector<double>& saturation,
std::vector<double>& surfacevol, std::vector<double>& surfacevol,
@ -1409,7 +1409,7 @@ namespace Opm
} }
void TransportModelCompressiblePolymer::scToc(const double* x, double* x_c) const { void TransportSolverTwophaseCompressiblePolymer::scToc(const double* x, double* x_c) const {
x_c[0] = x[0]; x_c[0] = x[0];
if (x[0] < 1e-2*tol_) { if (x[0] < 1e-2*tol_) {
x_c[1] = 0.5*polyprops_.cMax(); x_c[1] = 0.5*polyprops_.cMax();

18
opm/polymer/TransportModelCompressiblePolymer.hpp → opm/polymer/TransportSolverTwophaseCompressiblePolymer.hpp
View File

@ -17,11 +17,11 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED #ifndef OPM_TRANSPORTSOLVERTWOPHASECOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#define OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED #define OPM_TRANSPORTSOLVERTWOPHASECOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/core/transport/reorder/TransportModelInterface.hpp> #include <opm/core/transport/reorder/ReorderSolverInterface.hpp>
#include <opm/core/utility/linearInterpolation.hpp> #include <opm/core/utility/linearInterpolation.hpp>
#include <vector> #include <vector>
#include <list> #include <list>
@ -41,7 +41,7 @@ namespace Opm
/// Implements a reordering transport solver for incompressible two-phase flow /// Implements a reordering transport solver for incompressible two-phase flow
/// with polymer in the water phase. /// with polymer in the water phase.
/// \TODO Include permeability reduction effect. /// \TODO Include permeability reduction effect.
class TransportModelCompressiblePolymer : public TransportModelInterface class TransportSolverTwophaseCompressiblePolymer : public ReorderSolverInterface
{ {
public: public:
@ -59,7 +59,7 @@ namespace Opm
/// (using gradient variant and bracketing as fallbacks). /// (using gradient variant and bracketing as fallbacks).
/// \param[in] tol Tolerance used in the solver. /// \param[in] tol Tolerance used in the solver.
/// \param[in] maxit Maximum number of non-linear iterations used. /// \param[in] maxit Maximum number of non-linear iterations used.
TransportModelCompressiblePolymer(const UnstructuredGrid& grid, TransportSolverTwophaseCompressiblePolymer(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props, const BlackoilPropertiesInterface& props,
const PolymerProperties& polyprops, const PolymerProperties& polyprops,
const SingleCellMethod method, const SingleCellMethod method,
@ -183,9 +183,9 @@ namespace Opm
class ResSOnCurve; class ResSOnCurve;
class ResCOnCurve; class ResCOnCurve;
friend class TransportModelCompressiblePolymer::ResidualEquation; friend class TransportSolverTwophaseCompressiblePolymer::ResidualEquation;
friend class TransportModelCompressiblePolymer::ResSOnCurve; friend class TransportSolverTwophaseCompressiblePolymer::ResSOnCurve;
friend class TransportModelCompressiblePolymer::ResCOnCurve; friend class TransportSolverTwophaseCompressiblePolymer::ResCOnCurve;
virtual void solveSingleCell(const int cell); virtual void solveSingleCell(const int cell);
@ -231,4 +231,4 @@ namespace Opm
} // namespace Opm } // namespace Opm
#endif // OPM_TRANSPORTMODELCOMPRESSIBLEgPOLYMER_HEADER_INCLUDED #endif // OPM_TRANSPORTSOLVERTWOPHASECOMPRESSIBLEPOLYMER_HEADER_INCLUDED

122
opm/polymer/TransportModelPolymer.cpp → opm/polymer/TransportSolverTwophasePolymer.cpp
View File

@ -18,8 +18,8 @@
*/ */
#include <opm/polymer/TransportModelPolymer.hpp> #include <opm/polymer/TransportSolverTwophasePolymer.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp> #include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/utility/RootFinders.hpp> #include <opm/core/utility/RootFinders.hpp>
#include <opm/core/utility/miscUtilities.hpp> #include <opm/core/utility/miscUtilities.hpp>
@ -31,7 +31,7 @@
typedef Opm::RegulaFalsi<Opm::WarnAndContinueOnError> RootFinder; typedef Opm::RegulaFalsi<Opm::WarnAndContinueOnError> RootFinder;
class Opm::TransportModelPolymer::ResidualEquation class Opm::TransportSolverTwophasePolymer::ResidualEquation
{ {
public: public:
int cell; int cell;
@ -49,9 +49,9 @@ public:
double ads0; double ads0;
GradientMethod gradient_method; GradientMethod gradient_method;
TransportModelPolymer& tm; TransportSolverTwophasePolymer& tm;
ResidualEquation(TransportModelPolymer& tmodel, int cell_index); ResidualEquation(TransportSolverTwophasePolymer& tmodel, int cell_index);
void computeResidual(const double* x, double* res) const; void computeResidual(const double* x, double* res) const;
void computeResidual(const double* x, double* res, double& mc, double& ff) const; void computeResidual(const double* x, double* res, double& mc, double& ff) const;
double computeResidualS(const double* x) const; double computeResidualS(const double* x) const;
@ -67,9 +67,9 @@ private:
double* dres_c_dsdc, double& mc, double& ff) const; double* dres_c_dsdc, double& mc, double& ff) const;
}; };
class Opm::TransportModelPolymer::ResidualCGrav { class Opm::TransportSolverTwophasePolymer::ResidualCGrav {
public: public:
const TransportModelPolymer& tm; const TransportSolverTwophasePolymer& tm;
const int cell; const int cell;
const double s0; const double s0;
const double c0; const double c0;
@ -83,7 +83,7 @@ public:
int nbcell[2]; int nbcell[2];
mutable double last_s; mutable double last_s;
ResidualCGrav(const TransportModelPolymer& tmodel, ResidualCGrav(const TransportSolverTwophasePolymer& tmodel,
const std::vector<int>& cells, const std::vector<int>& cells,
const int pos, const int pos,
const double* gravflux); const double* gravflux);
@ -94,7 +94,7 @@ public:
double lastSaturation() const; double lastSaturation() const;
}; };
class Opm::TransportModelPolymer::ResidualSGrav { class Opm::TransportSolverTwophasePolymer::ResidualSGrav {
public: public:
const ResidualCGrav& res_c_eq_; const ResidualCGrav& res_c_eq_;
double c; double c;
@ -113,9 +113,9 @@ namespace
return std::max(std::abs(res[0]), std::abs(res[1])); return std::max(std::abs(res[0]), std::abs(res[1]));
} }
bool solveNewtonStepSC(const double* , const Opm::TransportModelPolymer::ResidualEquation&, bool solveNewtonStepSC(const double* , const Opm::TransportSolverTwophasePolymer::ResidualEquation&,
const double*, double*); const double*, double*);
bool solveNewtonStepC(const double* , const Opm::TransportModelPolymer::ResidualEquation&, bool solveNewtonStepC(const double* , const Opm::TransportSolverTwophasePolymer::ResidualEquation&,
const double*, double*); const double*, double*);
@ -150,11 +150,11 @@ namespace
class ResSOnCurve class ResSOnCurve
{ {
public: public:
ResSOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq); ResSOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq);
double operator()(const double t) const; double operator()(const double t) const;
CurveInSCPlane curve; CurveInSCPlane curve;
private: private:
const Opm::TransportModelPolymer::ResidualEquation& res_eq_; const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
}; };
// Compute the "c" residual along the curve "curve" for a given residual equation "res_eq". // Compute the "c" residual along the curve "curve" for a given residual equation "res_eq".
@ -162,11 +162,11 @@ namespace
class ResCOnCurve class ResCOnCurve
{ {
public: public:
ResCOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq); ResCOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq);
double operator()(const double t) const; double operator()(const double t) const;
CurveInSCPlane curve; CurveInSCPlane curve;
private: private:
const Opm::TransportModelPolymer::ResidualEquation& res_eq_; const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
}; };
} }
@ -174,7 +174,7 @@ namespace
namespace Opm namespace Opm
{ {
TransportModelPolymer::TransportModelPolymer(const UnstructuredGrid& grid, TransportSolverTwophasePolymer::TransportSolverTwophasePolymer(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props, const IncompPropertiesInterface& props,
const PolymerProperties& polyprops, const PolymerProperties& polyprops,
const SingleCellMethod method, const SingleCellMethod method,
@ -221,7 +221,7 @@ namespace Opm
void TransportModelPolymer::setPreferredMethod(SingleCellMethod method) void TransportSolverTwophasePolymer::setPreferredMethod(SingleCellMethod method)
{ {
method_ = method; method_ = method;
} }
@ -229,7 +229,7 @@ namespace Opm
void TransportModelPolymer::solve(const double* darcyflux, void TransportSolverTwophasePolymer::solve(const double* darcyflux,
const double* porevolume, const double* porevolume,
const double* source, const double* source,
const double* polymer_inflow_c, const double* polymer_inflow_c,
@ -262,11 +262,11 @@ namespace Opm
// //
// where influx is water influx, outflux is total outflux. // where influx is water influx, outflux is total outflux.
// Influxes are negative, outfluxes positive. // Influxes are negative, outfluxes positive.
struct TransportModelPolymer::ResidualS struct TransportSolverTwophasePolymer::ResidualS
{ {
TransportModelPolymer::ResidualEquation& res_eq_; TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
const double c_; const double c_;
explicit ResidualS(TransportModelPolymer::ResidualEquation& res_eq, explicit ResidualS(TransportSolverTwophasePolymer::ResidualEquation& res_eq,
const double c) const double c)
: res_eq_(res_eq), : res_eq_(res_eq),
c_(c) c_(c)
@ -287,11 +287,11 @@ namespace Opm
// \TODO doc me // \TODO doc me
// where ... // where ...
// Influxes are negative, outfluxes positive. // Influxes are negative, outfluxes positive.
struct TransportModelPolymer::ResidualC struct TransportSolverTwophasePolymer::ResidualC
{ {
mutable double s; // Mutable in order to change it with every operator() call to be the last computed s value. mutable double s; // Mutable in order to change it with every operator() call to be the last computed s value.
TransportModelPolymer::ResidualEquation& res_eq_; TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
explicit ResidualC(TransportModelPolymer::ResidualEquation& res_eq) explicit ResidualC(TransportSolverTwophasePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq) : res_eq_(res_eq)
{} {}
@ -335,7 +335,7 @@ namespace Opm
// ResidualEquation gathers parameters to construct the residual, computes its // ResidualEquation gathers parameters to construct the residual, computes its
// value and the values of its derivatives. // value and the values of its derivatives.
TransportModelPolymer::ResidualEquation::ResidualEquation(TransportModelPolymer& tmodel, int cell_index) TransportSolverTwophasePolymer::ResidualEquation::ResidualEquation(TransportSolverTwophasePolymer& tmodel, int cell_index)
: tm(tmodel) : tm(tmodel)
{ {
gradient_method = Analytic; gradient_method = Analytic;
@ -382,7 +382,7 @@ namespace Opm
} }
void TransportModelPolymer::ResidualEquation::computeResidual(const double* x, double* res) const void TransportSolverTwophasePolymer::ResidualEquation::computeResidual(const double* x, double* res) const
{ {
double dres_s_dsdc[2]; double dres_s_dsdc[2];
double dres_c_dsdc[2]; double dres_c_dsdc[2];
@ -391,7 +391,7 @@ namespace Opm
computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff); computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
} }
void TransportModelPolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const void TransportSolverTwophasePolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const
{ {
double dres_s_dsdc[2]; double dres_s_dsdc[2];
double dres_c_dsdc[2]; double dres_c_dsdc[2];
@ -399,7 +399,7 @@ namespace Opm
} }
double TransportModelPolymer::ResidualEquation::computeResidualS(const double* x) const double TransportSolverTwophasePolymer::ResidualEquation::computeResidualS(const double* x) const
{ {
double res[2]; double res[2];
double dres_s_dsdc[2]; double dres_s_dsdc[2];
@ -410,7 +410,7 @@ namespace Opm
return res[0]; return res[0];
} }
double TransportModelPolymer::ResidualEquation::computeResidualC(const double* x) const double TransportSolverTwophasePolymer::ResidualEquation::computeResidualC(const double* x) const
{ {
double res[2]; double res[2];
double dres_s_dsdc[2]; double dres_s_dsdc[2];
@ -421,7 +421,7 @@ namespace Opm
return res[1]; return res[1];
} }
void TransportModelPolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const void TransportSolverTwophasePolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const
// If gradient_method == FinDif, use finite difference // If gradient_method == FinDif, use finite difference
// If gradient_method == Analytic, use analytic expresions // If gradient_method == Analytic, use analytic expresions
{ {
@ -431,7 +431,7 @@ namespace Opm
computeResAndJacobi(x, true, true, true, false, res, gradient, dres_c_dsdc, mc, ff); computeResAndJacobi(x, true, true, true, false, res, gradient, dres_c_dsdc, mc, ff);
} }
void TransportModelPolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const void TransportSolverTwophasePolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const
// If gradient_method == FinDif, use finite difference // If gradient_method == FinDif, use finite difference
// If gradient_method == Analytic, use analytic expresions // If gradient_method == Analytic, use analytic expresions
{ {
@ -442,7 +442,7 @@ namespace Opm
} }
// Compute the Jacobian of the residual equations. // Compute the Jacobian of the residual equations.
void TransportModelPolymer::ResidualEquation::computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const void TransportSolverTwophasePolymer::ResidualEquation::computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const
{ {
double res[2]; double res[2];
double mc; double mc;
@ -450,7 +450,7 @@ namespace Opm
computeResAndJacobi(x, false, false, true, true, res, dres_s_dsdc, dres_c_dsdc, mc, ff); computeResAndJacobi(x, false, false, true, true, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
} }
void TransportModelPolymer::ResidualEquation::computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c, void TransportSolverTwophasePolymer::ResidualEquation::computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c,
const bool if_dres_s_dsdc, const bool if_dres_c_dsdc, const bool if_dres_s_dsdc, const bool if_dres_c_dsdc,
double* res, double* dres_s_dsdc, double* res, double* dres_s_dsdc,
double* dres_c_dsdc, double& mc, double& ff) const double* dres_c_dsdc, double& mc, double& ff) const
@ -550,7 +550,7 @@ namespace Opm
} }
} }
void TransportModelPolymer::solveSingleCell(const int cell) void TransportSolverTwophasePolymer::solveSingleCell(const int cell)
{ {
switch (method_) { switch (method_) {
case Bracketing: case Bracketing:
@ -574,7 +574,7 @@ namespace Opm
} }
void TransportModelPolymer::solveSingleCellBracketing(int cell) void TransportSolverTwophasePolymer::solveSingleCellBracketing(int cell)
{ {
ResidualEquation res_eq(*this, cell); ResidualEquation res_eq(*this, cell);
@ -606,7 +606,7 @@ namespace Opm
// Newton method, where we first try a Newton step. Then, if it does not work well, we look for // Newton method, where we first try a Newton step. Then, if it does not work well, we look for
// the zero of either the residual in s or the residual in c along a specified piecewise linear // the zero of either the residual in s or the residual in c along a specified piecewise linear
// curve. In these cases, we can use a robust 1d solver. // curve. In these cases, we can use a robust 1d solver.
void TransportModelPolymer::solveSingleCellGradient(int cell) void TransportSolverTwophasePolymer::solveSingleCellGradient(int cell)
{ {
int iters_used_falsi = 0; int iters_used_falsi = 0;
const int max_iters_split = maxit_; const int max_iters_split = maxit_;
@ -763,7 +763,7 @@ namespace Opm
} }
} }
void TransportModelPolymer::solveSingleCellNewton(int cell) void TransportSolverTwophasePolymer::solveSingleCellNewton(int cell)
{ {
const int max_iters_split = maxit_; const int max_iters_split = maxit_;
int iters_used_split = 0; int iters_used_split = 0;
@ -861,7 +861,7 @@ namespace Opm
} }
} }
void TransportModelPolymer::solveSingleCellNewtonSimple(int cell,bool use_sc) void TransportSolverTwophasePolymer::solveSingleCellNewtonSimple(int cell,bool use_sc)
{ {
const int max_iters_split = maxit_; const int max_iters_split = maxit_;
int iters_used_split = 0; int iters_used_split = 0;
@ -1010,7 +1010,7 @@ namespace Opm
void TransportModelPolymer::solveMultiCell(const int num_cells, const int* cells) void TransportSolverTwophasePolymer::solveMultiCell(const int num_cells, const int* cells)
{ {
double max_s_change = 0.0; double max_s_change = 0.0;
double max_c_change = 0.0; double max_c_change = 0.0;
@ -1067,21 +1067,21 @@ namespace Opm
<< num_iters << " iterations." << std::endl; << num_iters << " iterations." << std::endl;
} }
void TransportModelPolymer::fracFlow(double s, double c, double cmax, void TransportSolverTwophasePolymer::fracFlow(double s, double c, double cmax,
int cell, double& ff) const int cell, double& ff) const
{ {
double dummy[2]; double dummy[2];
fracFlowBoth(s, c, cmax, cell, ff, dummy, false); fracFlowBoth(s, c, cmax, cell, ff, dummy, false);
} }
void TransportModelPolymer::fracFlowWithDer(double s, double c, double cmax, void TransportSolverTwophasePolymer::fracFlowWithDer(double s, double c, double cmax,
int cell, double& ff, int cell, double& ff,
double* dff_dsdc) const double* dff_dsdc) const
{ {
fracFlowBoth(s, c, cmax, cell, ff, dff_dsdc, true); fracFlowBoth(s, c, cmax, cell, ff, dff_dsdc, true);
} }
void TransportModelPolymer::fracFlowBoth(double s, double c, double cmax, int cell, void TransportSolverTwophasePolymer::fracFlowBoth(double s, double c, double cmax, int cell,
double& ff, double* dff_dsdc, double& ff, double* dff_dsdc,
bool if_with_der) const bool if_with_der) const
{ {
@ -1111,12 +1111,12 @@ namespace Opm
} }
} }
void TransportModelPolymer::computeMc(double c, double& mc) const void TransportSolverTwophasePolymer::computeMc(double c, double& mc) const
{ {
polyprops_.computeMc(c, mc); polyprops_.computeMc(c, mc);
} }
void TransportModelPolymer::computeMcWithDer(double c, double& mc, void TransportSolverTwophasePolymer::computeMcWithDer(double c, double& mc,
double &dmc_dc) const double &dmc_dc) const
{ {
polyprops_.computeMcWithDer(c, mc, dmc_dc); polyprops_.computeMcWithDer(c, mc, dmc_dc);
@ -1124,14 +1124,14 @@ namespace Opm
TransportModelPolymer::ResidualSGrav::ResidualSGrav(const ResidualCGrav& res_c_eq, TransportSolverTwophasePolymer::ResidualSGrav::ResidualSGrav(const ResidualCGrav& res_c_eq,
const double c_init) const double c_init)
: res_c_eq_(res_c_eq), : res_c_eq_(res_c_eq),
c(c_init) c(c_init)
{ {
} }
double TransportModelPolymer::ResidualSGrav::operator()(double s) const double TransportSolverTwophasePolymer::ResidualSGrav::operator()(double s) const
{ {
return res_c_eq_.computeGravResidualS(s, c); return res_c_eq_.computeGravResidualS(s, c);
} }
@ -1145,7 +1145,7 @@ namespace Opm
// where ... // where ...
// Influxes are negative, outfluxes positive. // Influxes are negative, outfluxes positive.
TransportModelPolymer::ResidualCGrav::ResidualCGrav(const TransportModelPolymer& tmodel, TransportSolverTwophasePolymer::ResidualCGrav::ResidualCGrav(const TransportSolverTwophasePolymer& tmodel,
const std::vector<int>& cells, const std::vector<int>& cells,
const int pos, const int pos,
const double* gravflux) // Always oriented towards next in column. Size = colsize - 1. const double* gravflux) // Always oriented towards next in column. Size = colsize - 1.
@ -1177,7 +1177,7 @@ namespace Opm
tm.polyprops_.adsorption(c0, cmax0, c_ads0); tm.polyprops_.adsorption(c0, cmax0, c_ads0);
} }
double TransportModelPolymer::ResidualCGrav::operator()(double c) const double TransportSolverTwophasePolymer::ResidualCGrav::operator()(double c) const
{ {
ResidualSGrav res_s(*this); ResidualSGrav res_s(*this);
@ -1190,7 +1190,7 @@ namespace Opm
} }
double TransportModelPolymer::ResidualCGrav::computeGravResidualS(double s, double c) const double TransportSolverTwophasePolymer::ResidualCGrav::computeGravResidualS(double s, double c) const
{ {
double mobcell[2]; double mobcell[2];
@ -1216,7 +1216,7 @@ namespace Opm
return res; return res;
} }
double TransportModelPolymer::ResidualCGrav::computeGravResidualC(double s, double c) const double TransportSolverTwophasePolymer::ResidualCGrav::computeGravResidualC(double s, double c) const
{ {
double mobcell[2]; double mobcell[2];
@ -1247,13 +1247,13 @@ namespace Opm
return res; return res;
} }
double TransportModelPolymer::ResidualCGrav::lastSaturation() const double TransportSolverTwophasePolymer::ResidualCGrav::lastSaturation() const
{ {
return last_s; return last_s;
} }
void TransportModelPolymer::mobility(double s, double c, int cell, double* mob) const void TransportSolverTwophasePolymer::mobility(double s, double c, int cell, double* mob) const
{ {
double sat[2] = { s, 1.0 - s }; double sat[2] = { s, 1.0 - s };
double relperm[2]; double relperm[2];
@ -1262,7 +1262,7 @@ namespace Opm
} }
void TransportModelPolymer::initGravity(const double* grav) void TransportSolverTwophasePolymer::initGravity(const double* grav)
{ {
// Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j) // Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j)
std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]); std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
@ -1285,7 +1285,7 @@ namespace Opm
} }
void TransportModelPolymer::solveSingleCellGravity(const std::vector<int>& cells, void TransportSolverTwophasePolymer::solveSingleCellGravity(const std::vector<int>& cells,
const int pos, const int pos,
const double* gravflux) const double* gravflux)
{ {
@ -1312,7 +1312,7 @@ namespace Opm
mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]); mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]);
} }
int TransportModelPolymer::solveGravityColumn(const std::vector<int>& cells) int TransportSolverTwophasePolymer::solveGravityColumn(const std::vector<int>& cells)
{ {
// Set up column gravflux. // Set up column gravflux.
const int nc = cells.size(); const int nc = cells.size();
@ -1372,7 +1372,7 @@ namespace Opm
} }
void TransportModelPolymer::solveGravity(const std::vector<std::vector<int> >& columns, void TransportSolverTwophasePolymer::solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume, const double* porevolume,
const double dt, const double dt,
std::vector<double>& saturation, std::vector<double>& saturation,
@ -1410,7 +1410,7 @@ namespace Opm
toBothSat(saturation_, saturation); toBothSat(saturation_, saturation);
} }
void TransportModelPolymer::scToc(const double* x, double* x_c) const { void TransportSolverTwophasePolymer::scToc(const double* x, double* x_c) const {
x_c[0] = x[0]; x_c[0] = x[0];
if (x[0] < 1e-2*tol_) { if (x[0] < 1e-2*tol_) {
x_c[1] = 0.5*polyprops_.cMax(); x_c[1] = 0.5*polyprops_.cMax();
@ -1523,7 +1523,7 @@ namespace
} }
ResSOnCurve::ResSOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq) ResSOnCurve::ResSOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq) : res_eq_(res_eq)
{ {
} }
@ -1537,7 +1537,7 @@ namespace
return res_eq_.computeResidualS(x_c); return res_eq_.computeResidualS(x_c);
} }
ResCOnCurve::ResCOnCurve(const Opm::TransportModelPolymer::ResidualEquation& res_eq) ResCOnCurve::ResCOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq)
: res_eq_(res_eq) : res_eq_(res_eq)
{ {
} }
@ -1551,7 +1551,7 @@ namespace
return res_eq_.computeResidualC(x_c); return res_eq_.computeResidualC(x_c);
} }
bool solveNewtonStepSC(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq, bool solveNewtonStepSC(const double* xx, const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq,
const double* res, double* x_new) { const double* res, double* x_new) {
double dres_s_dsdc[2]; double dres_s_dsdc[2];
@ -1582,7 +1582,7 @@ namespace
return true; return true;
} }
} }
bool solveNewtonStepC(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq, bool solveNewtonStepC(const double* xx, const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq,
const double* res, double* x_new) { const double* res, double* x_new) {
double dres_s_dsdc[2]; double dres_s_dsdc[2];

22
opm/polymer/TransportModelPolymer.hpp → opm/polymer/TransportSolverTwophasePolymer.hpp
View File

@ -17,11 +17,11 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef OPM_TRANSPORTMODELPOLYMER_HEADER_INCLUDED #ifndef OPM_TRANSPORTSOLVERTWOPHASEPOLYMER_HEADER_INCLUDED
#define OPM_TRANSPORTMODELPOLYMER_HEADER_INCLUDED #define OPM_TRANSPORTSOLVERTWOPHASEPOLYMER_HEADER_INCLUDED
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/core/transport/reorder/TransportModelInterface.hpp> #include <opm/core/transport/reorder/ReorderSolverInterface.hpp>
#include <opm/core/utility/linearInterpolation.hpp> #include <opm/core/utility/linearInterpolation.hpp>
#include <vector> #include <vector>
#include <list> #include <list>
@ -36,7 +36,7 @@ namespace Opm
/// Implements a reordering transport solver for incompressible two-phase flow /// Implements a reordering transport solver for incompressible two-phase flow
/// with polymer in the water phase. /// with polymer in the water phase.
/// \TODO Include permeability reduction effect. /// \TODO Include permeability reduction effect.
class TransportModelPolymer : public TransportModelInterface class TransportSolverTwophasePolymer : public ReorderSolverInterface
{ {
public: public:
@ -53,12 +53,12 @@ namespace Opm
/// (using gradient variant and bracketing as fallbacks). /// (using gradient variant and bracketing as fallbacks).
/// \param[in] tol Tolerance used in the solver. /// \param[in] tol Tolerance used in the solver.
/// \param[in] maxit Maximum number of non-linear iterations used. /// \param[in] maxit Maximum number of non-linear iterations used.
TransportModelPolymer(const UnstructuredGrid& grid, TransportSolverTwophasePolymer(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props, const IncompPropertiesInterface& props,
const PolymerProperties& polyprops, const PolymerProperties& polyprops,
const SingleCellMethod method, const SingleCellMethod method,
const double tol, const double tol,
const int maxit); const int maxit);
/// Set the preferred method, Bracketing or Newton. /// Set the preferred method, Bracketing or Newton.
void setPreferredMethod(SingleCellMethod method); void setPreferredMethod(SingleCellMethod method);
@ -191,4 +191,4 @@ namespace Opm
} // namespace Opm } // namespace Opm
#endif // OPM_TRANSPORTMODELPOLYMER_HEADER_INCLUDED #endif // OPM_TRANSPORTSOLVERTWOPHASEPOLYMER_HEADER_INCLUDED

6
opm/polymer/polymerUtilities.hpp
View File

@ -22,12 +22,12 @@
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/fluid/IncompPropertiesInterface.hpp> #include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp> #include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/PolymerState.hpp> #include <opm/polymer/PolymerState.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp> #include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/core/fluid/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/utility/SparseVector.hpp> #include <opm/core/utility/SparseVector.hpp>
#include <vector> #include <vector>