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

Added IncompTpfaPolymer class.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-06-13 19:56:44 +02:00
parent 23580678dc
commit 7a40153e51
6 changed files with 316 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile.am
View File

@ -7,6 +7,7 @@ SUBDIRS = . examples
lib_LTLIBRARIES = libopmpolymer.la
libopmpolymer_la_SOURCES = \
opm/polymer/IncompTpfaPolymer.cpp \
opm/polymer/SimulatorPolymer.cpp \
opm/polymer/TransportModelPolymer.cpp \
opm/polymer/PolymerProperties.cpp \
@ -16,6 +17,7 @@ nobase_include_HEADERS = \
opm/polymer/GravityColumnSolverPolymer.hpp \
opm/polymer/GravityColumnSolverPolymer_impl.hpp \
opm/polymer/IncompPropertiesDefaultPolymer.hpp \
opm/polymer/IncompTpfaPolymer.hpp \
opm/polymer/PolymerProperties.hpp \
opm/polymer/PolymerState.hpp \
opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp \

10
examples/polymer_reorder.cpp
View File

@ -22,7 +22,7 @@
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
@ -575,10 +575,10 @@ main(int argc, char** argv)
Opm::LinearSolverFactory linsolver(param);
// Pressure solver.
const double *grav = use_gravity ? &gravity[0] : 0;
Opm::IncompTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
nl_pressure_residual_tolerance, nl_pressure_change_tolerance,
nl_pressure_maxiter,
grav, wells->c_wells(), src, bcs.c_bcs());
Opm::IncompTpfaPolymer psolver(*grid->c_grid(), *props, rock_comp.get(), polyprop, linsolver,
nl_pressure_residual_tolerance, nl_pressure_change_tolerance,
nl_pressure_maxiter,
grav, wells->c_wells(), src, bcs.c_bcs());
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);

168
opm/polymer/IncompTpfaPolymer.cpp Normal file
View File

@ -0,0 +1,168 @@
/*
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/>.
*/
#include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/mimetic/mimetic.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/polymer/PolymerState.hpp>
#include <opm/polymer/PolymerUtilities.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/newwells.h>
#include <iomanip>
#include <cmath>
#include <algorithm>
namespace Opm
{
/// Construct solver, possibly with rock compressibility.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] rock_comp_props Rock compressibility properties. May be null.
/// \param[in] linsolver Linear solver to use.
/// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller.
/// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller.
/// \param[in] maxiter Maximum acceptable number of iterations.
/// \param[in] gravity Gravity vector. If non-null, the array should
/// have D elements.
/// \param[in] wells The wells argument. Will be used in solution,
/// is ignored if NULL.
/// Note: this class observes the well object, and
/// makes the assumption that the well topology
/// and completions does not change during the
/// run. However, controls (only) are allowed
/// to change.
/// \param[in] src Source terms. May be empty().
/// \param[in] bcs Boundary conditions, treat as all noflow if null.
IncompTpfaPolymer::IncompTpfaPolymer(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
const PolymerProperties& poly_props,
LinearSolverInterface& linsolver,
const double residual_tol,
const double change_tol,
const int maxiter,
const double* gravity,
const Wells* wells,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs)
: IncompTpfa(grid, props, rock_comp_props, linsolver,
residual_tol, change_tol, maxiter,
gravity, wells, src, bcs),
poly_props_(poly_props),
c_(0),
cmax_(0)
{
}
/// Solve the pressure equation. If there is no pressure
/// dependency introduced by rock compressibility effects,
/// the equation is linear, and it is solved directly.
/// Otherwise, the nonlinear equations ares solved by a
/// Newton-Raphson scheme.
/// May throw an exception if the number of iterations
/// exceed maxiter (set in constructor).
void IncompTpfaPolymer::solve(const double dt,
PolymerState& state,
WellState& well_state)
{
c_ = &state.concentration();
cmax_ = &state.maxconcentration();
if (rock_comp_props_ != 0 && rock_comp_props_->isActive()) {
solveRockComp(dt, state.twophaseState(), well_state);
} else {
solveIncomp(dt, state.twophaseState(), well_state);
}
}
/// Compute per-solve dynamic properties.
void IncompTpfaPolymer::computePerSolveDynamicData(const double /*dt*/,
const TwophaseState& state,
const WellState& /*well_state*/)
{
// Computed here:
//
// std::vector<double> wdp_;
// std::vector<double> totmob_;
// std::vector<double> omega_;
// std::vector<double> trans_;
// std::vector<double> gpress_omegaweighted_;
// std::vector<double> initial_porevol_;
// ifs_tpfa_forces forces_;
// The only difference from IncompTpfa::computePerSolveDynamicData() is that
// we call the polymer-aware versions of the computeTotalMobility*() functions.
// wdp_
if (wells_) {
Opm::computeWDP(*wells_, grid_, state.saturation(), props_.density(),
gravity_ ? gravity_[2] : 0.0, true, wdp_);
}
// totmob_, omega_, gpress_omegaweighted_
if (gravity_) {
computeTotalMobilityOmega(props_, poly_props_, allcells_, state.saturation(), *c_, *cmax_,
totmob_, omega_);
mim_ip_density_update(grid_.number_of_cells, grid_.cell_facepos,
&omega_[0],
&gpress_[0], &gpress_omegaweighted_[0]);
} else {
computeTotalMobility(props_, poly_props_, allcells_, state.saturation(), *c_, *cmax_, totmob_);
}
// trans_
tpfa_eff_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &totmob_[0], &htrans_[0], &trans_[0]);
// initial_porevol_
if (rock_comp_props_ && rock_comp_props_->isActive()) {
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_);
}
// forces_
forces_.src = src_.empty() ? NULL : &src_[0];
forces_.bc = bcs_;
forces_.W = wells_;
forces_.totmob = &totmob_[0];
forces_.wdp = wdp_.empty() ? NULL : &wdp_[0];
}
} // namespace Opm

111
opm/polymer/IncompTpfaPolymer.hpp Normal file
View File

@ -0,0 +1,111 @@
/*
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/>.
*/
#ifndef OPM_INCOMPTPFAPOLYMER_HEADER_INCLUDED
#define OPM_INCOMPTPFAPOLYMER_HEADER_INCLUDED
#include <opm/core/pressure/IncompTpfa.hpp>
#include <vector>
struct UnstructuredGrid;
struct Wells;
struct FlowBoundaryConditions;
namespace Opm
{
class IncompPropertiesInterface;
class RockCompressibility;
class PolymerProperties;
class LinearSolverInterface;
class PolymerState;
class WellState;
/// Encapsulating a tpfa pressure solver for the incompressible-fluid case with polymer.
/// Supports gravity, wells controlled by bhp or reservoir rates,
/// boundary conditions and simple sources as driving forces.
/// Rock compressibility can be included, and necessary nonlinear
/// iterations are handled.
/// Below we use the shortcuts D for the number of dimensions, N
/// for the number of cells and F for the number of faces.
class IncompTpfaPolymer : public IncompTpfa
{
public:
/// Construct solver, possibly with rock compressibility.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] rock_comp_props Rock compressibility properties. May be null.
/// \param[in] poly_props Polymer properties.
/// \param[in] linsolver Linear solver to use.
/// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller.
/// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller.
/// \param[in] maxiter Maximum acceptable number of iterations.
/// \param[in] gravity Gravity vector. If non-null, the array should
/// have D elements.
/// \param[in] wells The wells argument. Will be used in solution,
/// is ignored if NULL.
/// Note: this class observes the well object, and
/// makes the assumption that the well topology
/// and completions does not change during the
/// run. However, controls (only) are allowed
/// to change.
/// \param[in] src Source terms. May be empty().
/// \param[in] bcs Boundary conditions, treat as all noflow if null.
IncompTpfaPolymer(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
const PolymerProperties& poly_props,
LinearSolverInterface& linsolver,
const double residual_tol,
const double change_tol,
const int maxiter,
const double* gravity,
const Wells* wells,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs);
/// Solve the pressure equation. If there is no pressure
/// dependency introduced by rock compressibility effects,
/// the equation is linear, and it is solved directly.
/// Otherwise, the nonlinear equations ares solved by a
/// Newton-Raphson scheme.
/// May throw an exception if the number of iterations
/// exceed maxiter (set in constructor).
void solve(const double dt,
PolymerState& state,
WellState& well_state);
private:
virtual void computePerSolveDynamicData(const double dt,
const TwophaseState& state,
const WellState& well_state);
private:
// ------ Data that will remain unmodified after construction. ------
const PolymerProperties& poly_props_;
// ------ Data that will be updated every solve() call. ------
const std::vector<double>* c_;
const std::vector<double>* cmax_;
};
} // namespace Opm
#endif // OPM_INCOMPTPFAPOLYMER_HEADER_INCLUDED

33
opm/polymer/PolymerState.hpp
View File

@ -29,25 +29,46 @@ namespace Opm
{
/// Simulator state for a two-phase simulator with polymer.
/// Inheritance is not used for runtime polymorphism,
/// just for simplicity of implementation.
class PolymerState : public TwophaseState
class PolymerState
{
public:
void init(const UnstructuredGrid& g, int num_phases)
{
TwophaseState::init(g, num_phases);
state2p_.init(g, num_phases);
concentration_.resize(g.number_of_cells, 0.0);
cmax_.resize(g.number_of_cells, 0.0);
}
std::vector<double>& concentration() { return concentration_; }
enum ExtremalSat { MinSat = TwophaseState::MinSat, MaxSat = TwophaseState::MaxSat };
void setFirstSat(const std::vector<int>& cells,
const Opm::IncompPropertiesInterface& props,
ExtremalSat es)
{
// A better solution for embedding TwophaseState::ExtremalSat could perhaps
// be found, to avoid the cast.
state2p_.setFirstSat(cells, props, static_cast<TwophaseState::ExtremalSat>(es));
}
std::vector<double>& pressure () { return state2p_.pressure(); }
std::vector<double>& facepressure() { return state2p_.facepressure(); }
std::vector<double>& faceflux () { return state2p_.faceflux(); }
std::vector<double>& saturation () { return state2p_.saturation(); }
std::vector<double>& concentration() { return concentration_; }
std::vector<double>& maxconcentration() { return cmax_; }
const std::vector<double>& concentration() const { return concentration_; }
const std::vector<double>& pressure () const { return state2p_.pressure(); }
const std::vector<double>& facepressure() const { return state2p_.facepressure(); }
const std::vector<double>& faceflux () const { return state2p_.faceflux(); }
const std::vector<double>& saturation () const { return state2p_.saturation(); }
const std::vector<double>& concentration() const { return concentration_; }
const std::vector<double>& maxconcentration() const { return cmax_; }
TwophaseState& twophaseState() { return state2p_; }
const TwophaseState& twophaseState() const { return state2p_; }
private:
TwophaseState state2p_;
std::vector<double> concentration_;
std::vector<double> cmax_;
};

6
opm/polymer/SimulatorPolymer.cpp
View File

@ -26,7 +26,7 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/core/grid.h>
#include <opm/core/newwells.h>
@ -114,7 +114,7 @@ namespace Opm
const LinearSolverInterface& linsolver_;
const double* gravity_;
// Solvers
IncompTpfa psolver_;
IncompTpfaPolymer psolver_;
TransportModelPolymer tsolver_;
// Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_;
@ -175,7 +175,7 @@ namespace Opm
bcs_(bcs),
linsolver_(linsolver),
gravity_(gravity),
psolver_(grid, props, rock_comp_props, linsolver,
psolver_(grid, props, rock_comp_props, poly_props, linsolver,
param.getDefault("nl_pressure_residual_tolerance", 1e-8),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),