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 #737 from atgeirr/sequential

Browse Source 
Experimental sequential solvers
This commit is contained in:
Atgeirr Flø Rasmussen 2016-06-27 11:38:13 +02:00 committed by GitHub
commit 80b8b5f12f
9 changed files with 1613 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -148,5 +148,6 @@ if (HAVE_OPM_DATA)
set_tests_properties(compare_restart_files PROPERTIES DEPENDS flow_SPE1CASE2_restart) # Compares the restart files from tests flow_SPE1CASE2_restart and flow_SPE1CASE2
add_test( NAME flow_sequential_SPE1 COMMAND flow_sequential ${OPM_DATA_ROOT}/spe1/SPE1CASE1.DATA )
endif()

7
CMakeLists_files.cmake
View File

@ -101,6 +101,7 @@ list (APPEND TEST_DATA_FILES
list (APPEND EXAMPLE_SOURCE_FILES
examples/find_zero.cpp
examples/flow.cpp
examples/flow_sequential.cpp
examples/flow_multisegment.cpp
examples/flow_solvent.cpp
examples/sim_2p_incomp.cpp
@ -121,6 +122,7 @@ list (APPEND PROGRAM_SOURCE_FILES
examples/sim_2p_incomp_ad.cpp
examples/sim_2p_comp_reorder.cpp
examples/flow.cpp
examples/flow_sequential.cpp
examples/flow_solvent.cpp
examples/opm_init_check.cpp
examples/sim_poly2p_comp_reorder.cpp
@ -143,21 +145,25 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/BlackoilModelBase_impl.hpp
opm/autodiff/BlackoilModelEnums.hpp
opm/autodiff/BlackoilModelParameters.hpp
opm/autodiff/BlackoilPressureModel.hpp
opm/autodiff/BlackoilPropsAdFromDeck.hpp
opm/autodiff/SolventPropsAdFromDeck.hpp
opm/autodiff/BlackoilPropsAdInterface.hpp
opm/autodiff/CPRPreconditioner.hpp
opm/autodiff/createGlobalCellArray.hpp
opm/autodiff/BlackoilSequentialModel.hpp
opm/autodiff/BlackoilSolventModel.hpp
opm/autodiff/BlackoilSolventModel_impl.hpp
opm/autodiff/BlackoilSolventState.hpp
opm/autodiff/BlackoilMultiSegmentModel.hpp
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
opm/autodiff/BlackoilTransportModel.hpp
opm/autodiff/fastSparseOperations.hpp
opm/autodiff/DuneMatrix.hpp
opm/autodiff/ExtractParallelGridInformationToISTL.hpp
opm/autodiff/FlowMain.hpp
opm/autodiff/FlowMainPolymer.hpp
opm/autodiff/FlowMainSequential.hpp
opm/autodiff/FlowMainSolvent.hpp
opm/autodiff/GeoProps.hpp
opm/autodiff/GridHelpers.hpp
@ -185,6 +191,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
opm/autodiff/SimulatorIncompTwophaseAd.hpp
opm/autodiff/SimulatorSequentialBlackoil.hpp
opm/autodiff/TransportSolverTwophaseAd.hpp
opm/autodiff/WellDensitySegmented.hpp
opm/autodiff/WellStateFullyImplicitBlackoil.hpp

39
examples/flow_sequential.cpp Normal file
View File

@ -0,0 +1,39 @@
/*
Copyright 2016 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/grid.h>
#include <opm/autodiff/SimulatorSequentialBlackoil.hpp>
#include <opm/autodiff/FlowMainSequential.hpp>
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
typedef UnstructuredGrid Grid;
typedef Opm::SimulatorSequentialBlackoil<Grid> Simulator;
Opm::FlowMainSequential<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}

330
opm/autodiff/BlackoilPressureModel.hpp Normal file
View File

@ -0,0 +1,330 @@
/*
Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil AS.
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_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
#define OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED
#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <algorithm>
namespace Opm {
/// A model implementation for the pressure equation in three-phase black oil.
///
/// The model is based on the normal black oil model.
/// It uses automatic differentiation via the class AutoDiffBlock
/// to simplify assembly of the jacobian matrix.
template<class Grid, class WellModel>
class BlackoilPressureModel : public BlackoilModelBase<Grid, WellModel, BlackoilPressureModel<Grid, WellModel> >
{
public:
typedef BlackoilModelBase<Grid, WellModel, BlackoilPressureModel<Grid, WellModel> > Base;
friend Base;
typedef typename Base::ReservoirState ReservoirState;
typedef typename Base::WellState WellState;
typedef typename Base::SolutionState SolutionState;
typedef typename Base::V V;
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
/// \param[in] param parameters
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells_arg well structure
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
BlackoilPressureModel(const typename Base::ModelParameters& param,
const Grid& grid,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const StandardWells& std_wells,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
eclState, has_disgas, has_vapoil, terminal_output),
state0_(3),
max_dp_rel_(std::numeric_limits<double>::infinity()),
scaling_{ ADB::null(), ADB::null(), ADB::null() }
{
}
/// Called once per timestep.
void prepareStep(const double dt,
const ReservoirState& reservoir_state,
const WellState& well_state)
{
asImpl().wellModel().setStoreWellPerforationFluxesFlag(true);
Base::prepareStep(dt, reservoir_state, well_state);
max_dp_rel_ = std::numeric_limits<double>::infinity();
state0_ = asImpl().variableState(reservoir_state, well_state);
asImpl().makeConstantState(state0_);
}
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
V solveJacobianSystem() const
{
// We make a reduced residual object which just
// contains the pressure residual.
// TODO: directly contruct that object in residual_ instead.
const int n1 = residual_.material_balance_eq[0].size();
const int n2 = residual_.well_flux_eq.size() + residual_.well_eq.size();
const int n_full = residual_.sizeNonLinear();
const auto& mb = residual_.material_balance_eq;
const auto& fe = residual_.well_flux_eq;
const auto& we = residual_.well_eq;
LinearisedBlackoilResidual pressure_res = {
{
// TODO: handle general 2-phase etc.
ADB::function(mb[0].value(), { mb[0].derivative()[0], mb[0].derivative()[3], mb[0].derivative()[4] })
},
ADB::function(fe.value(), { fe.derivative()[0], fe.derivative()[3], fe.derivative()[4] }),
ADB::function(we.value(), { we.derivative()[0], we.derivative()[3], we.derivative()[4] }),
residual_.matbalscale,
residual_.singlePrecision
};
assert(pressure_res.sizeNonLinear() == n1 + n2);
V dx_pressure = linsolver_.computeNewtonIncrement(pressure_res);
assert(dx_pressure.size() == n1 + n2);
V dx_full = V::Zero(n_full);
dx_full.topRows(n1) = dx_pressure.topRows(n1);
dx_full.bottomRows(n2) = dx_pressure.bottomRows(n2);
return dx_full;
}
using Base::numPhases;
using Base::numMaterials;
protected:
using Base::asImpl;
using Base::linsolver_;
using Base::residual_;
using Base::rq_;
using Base::grid_;
using Base::ops_;
using Base::has_vapoil_;
using Base::has_disgas_;
SolutionState state0_;
double max_dp_rel_ = std::numeric_limits<double>::infinity();
ADB scaling_[3] = { ADB::null(), ADB::null(), ADB::null() };
void assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly)
{
Base::assemble(reservoir_state, well_state, initial_assembly);
if (initial_assembly) {
}
// Compute pressure residual.
ADB pressure_residual = ADB::constant(V::Zero(residual_.material_balance_eq[0].size()));
for (int phase = 0; phase < numPhases(); ++phase) {
pressure_residual += residual_.material_balance_eq[phase] * scaling_[phase];
}
residual_.material_balance_eq[0] = pressure_residual; // HACK
// Compute total reservoir volume flux.
const int n = rq_[0].mflux.size();
V flux = V::Zero(n);
for (int phase = 0; phase < numPhases(); ++phase) {
UpwindSelector<double> upwind(grid_, ops_, rq_[phase].dh.value());
flux += rq_[phase].mflux.value() / upwind.select(rq_[phase].b.value());
}
// Storing the fluxes in the assemble() method is a bit of
// a hack, but alternatives either require a more
// significant redesign of the base class or are
// inefficient.
ReservoirState& s = const_cast<ReservoirState&>(reservoir_state);
s.faceflux().resize(n);
std::copy_n(flux.data(), n, s.faceflux().begin());
if (asImpl().localWellsActive()) {
const V& wflux = asImpl().wellModel().getStoredWellPerforationFluxes();
assert(int(well_state.perfRates().size()) == wflux.size());
std::copy_n(wflux.data(), wflux.size(), well_state.perfRates().begin());
}
}
SolutionState
variableState(const ReservoirState& x,
const WellState& xw) const
{
// As Base::variableState(), except making Sw and Xvar constants.
std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
std::vector<ADB> vars = ADB::variables(vars0);
const std::vector<int> indices = asImpl().variableStateIndices();
vars[indices[Sw]] = ADB::constant(vars[indices[Sw]].value());
vars[indices[Xvar]] = ADB::constant(vars[indices[Xvar]].value());
// OpmLog::debug("Injector pressure: " + std::to_string(vars[indices[Bhp]].value()[1]));
return asImpl().variableStateExtractVars(x, indices, vars);
}
void computeAccum(const SolutionState& state,
const int aix )
{
if (aix == 0) {
Base::computeAccum(state0_, aix);
} else {
Base::computeAccum(state, aix);
}
}
void assembleMassBalanceEq(const SolutionState& state)
{
Base::assembleMassBalanceEq(state);
// Compute the scaling factors.
// Scaling factors are:
// 1/bw, 1/bo - rs/bg, 1/bg - rv/bo
assert(numPhases() == 3);
assert(numMaterials() == 3);
V one = V::Constant(state.pressure.size(), 1.0);
scaling_[Water] = one / rq_[Water].b;
scaling_[Oil] = one / rq_[Oil].b - state.rs / rq_[Gas].b;
scaling_[Gas] = one / rq_[Gas].b - state.rv / rq_[Oil].b;
if (has_disgas_ && has_vapoil_) {
ADB r_factor = one / (one - state.rs * state.rv);
scaling_[Oil] = r_factor * scaling_[Oil];
scaling_[Gas] = r_factor * scaling_[Gas];
}
// @TODO: multiply the oil and gas scale by 1/(1-rs*rv)?
// OpmLog::debug("gas scaling: " + std::to_string(scaling_[2].value()[0]));
}
void updateState(const V& dx,
ReservoirState& reservoir_state,
WellState& well_state)
{
// Naively, rs and rv can get overwritten, so we
// avoid that by storing.
std::vector<double> rs_old = reservoir_state.gasoilratio();
std::vector<double> rv_old = reservoir_state.rv();
auto hs_old = reservoir_state.hydroCarbonState();
auto phasecond_old = Base::phaseCondition_;
auto isRs_old = Base::isRs_;
auto isRv_old = Base::isRv_;
auto isSg_old = Base::isSg_;
// Compute the pressure range.
const auto minmax_iters = std::minmax_element(reservoir_state.pressure().begin(),
reservoir_state.pressure().end());
const double range = *minmax_iters.second - *minmax_iters.first;
// Use the base class' updateState().
Base::updateState(dx, reservoir_state, well_state);
// Compute relative change.
max_dp_rel_ = dx.head(reservoir_state.pressure().size()).abs().maxCoeff() / range;
// Restore rs and rv, also various state flags.
reservoir_state.gasoilratio() = rs_old;
reservoir_state.rv() = rv_old;
reservoir_state.hydroCarbonState() = hs_old;
Base::phaseCondition_ = phasecond_old;
Base::isRs_ = isRs_old;
Base::isRv_ = isRv_old;
Base::isSg_ = isSg_old;
}
bool getConvergence(const double /* dt */, const int iteration)
{
const double tol_p = 1e-11;
const double resmax = residual_.material_balance_eq[0].value().abs().maxCoeff();
if (Base::terminalOutputEnabled()) {
// Only rank 0 does print to std::cout
if (iteration == 0) {
OpmLog::info("\nIter Res(p) Delta(p)\n");
}
std::ostringstream os;
os.precision(3);
os.setf(std::ios::scientific);
os << std::setw(4) << iteration;
os << std::setw(11) << resmax;
os << std::setw(11) << max_dp_rel_;
OpmLog::info(os.str());
}
return resmax < tol_p;
}
};
/// Providing types by template specialisation of ModelTraits for BlackoilPressureModel.
template <class Grid, class WellModel>
struct ModelTraits< BlackoilPressureModel<Grid, WellModel> >
{
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters;
typedef DefaultBlackoilSolutionState SolutionState;
};
} // namespace Opm
#endif // OPM_BLACKOILPRESSUREMODEL_HEADER_INCLUDED

267
opm/autodiff/BlackoilSequentialModel.hpp Normal file
View File

@ -0,0 +1,267 @@
/*
Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil AS.
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_BLACKOILSEQUENTIALMODEL_HEADER_INCLUDED
#define OPM_BLACKOILSEQUENTIALMODEL_HEADER_INCLUDED
#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/autodiff/BlackoilPressureModel.hpp>
#include <opm/autodiff/BlackoilTransportModel.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
namespace Opm {
struct BlackoilSequentialModelParameters : public BlackoilModelParameters
{
bool iterate_to_fully_implicit;
explicit BlackoilSequentialModelParameters( const parameter::ParameterGroup& param )
: BlackoilModelParameters(param),
iterate_to_fully_implicit(param.getDefault("iterate_to_fully_implicit", false))
{
}
};
/// A sequential splitting model implementation for three-phase black oil.
template<class Grid, class WellModel>
class BlackoilSequentialModel
{
public:
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilSequentialModelParameters ModelParameters;
typedef DefaultBlackoilSolutionState SolutionState;
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
/// \param[in] param parameters
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] vfp_properties Vertical flow performance tables
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
BlackoilSequentialModel(const ModelParameters& param,
const Grid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const WellModel well_model,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output)
: pressure_model_(new PressureModel(param, grid, fluid, geo, rock_comp_props, well_model,
linsolver, eclState, has_disgas, has_vapoil, terminal_output)),
transport_model_(new TransportModel(param, grid, fluid, geo, rock_comp_props, well_model,
linsolver, eclState, has_disgas, has_vapoil, terminal_output)),
pressure_solver_(typename PressureSolver::SolverParameters(), std::move(pressure_model_)),
transport_solver_(typename TransportSolver::SolverParameters(), std::move(transport_model_)),
initial_reservoir_state_(0, 0, 0), // will be overwritten
iterate_to_fully_implicit_(param.iterate_to_fully_implicit)
{
}
/// Called once before each time step.
/// \param[in] dt time step size
/// \param[in] reservoir_state reservoir state variables
/// \param[in] well_state well state variables
void prepareStep(const double /* dt */,
const ReservoirState& reservoir_state,
const WellState& well_state)
{
initial_reservoir_state_ = reservoir_state;
initial_well_state_ = well_state;
}
/// Called once per nonlinear iteration.
/// This model will first solve the pressure model to convergence, then the
/// transport model.
/// \param[in] iteration should be 0 for the first call of a new timestep
/// \param[in] dt time step size
/// \param[in] nonlinear_solver nonlinear solver used (for oscillation/relaxation control)
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
template <class NonlinearSolverType>
IterationReport nonlinearIteration(const int iteration,
const double dt,
NonlinearSolverType& /* nonlinear_solver */,
ReservoirState& reservoir_state,
WellState& well_state)
{
if (!iterate_to_fully_implicit_) {
// Do a single pressure solve, followed by a single transport solve.
if (terminalOutputEnabled()) {
OpmLog::info("Using sequential model.");
}
// Pressure solve.
if (terminalOutputEnabled()) {
OpmLog::info("Solving the pressure equation.");
}
ReservoirState initial_state = reservoir_state;
const int pressure_liniter = pressure_solver_.step(dt, reservoir_state, well_state);
if (pressure_liniter == -1) {
OPM_THROW(std::runtime_error, "Pressure solver failed to converge.");
}
// Transport solve.
if (terminalOutputEnabled()) {
OpmLog::info("Solving the transport equations.");
}
const int transport_liniter = transport_solver_.step(dt, initial_state, well_state, reservoir_state, well_state);
if (transport_liniter == -1) {
OPM_THROW(std::runtime_error, "Transport solver failed to converge.");
}
// Report and return.
return IterationReport { false, true, pressure_liniter + transport_liniter };
} else {
// Iterate to fully implicit solution.
// This call is just for a single iteration (one pressure and one transport solve),
// we return a 'false' converged status if more are needed
if (terminalOutputEnabled()) {
OpmLog::info("Using sequential model in iterative mode, outer iteration " + std::to_string(iteration));
}
// Pressure solve.
if (terminalOutputEnabled()) {
OpmLog::info("Solving the pressure equation.");
}
const int pressure_liniter = pressure_solver_.step(dt, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
if (pressure_liniter == -1) {
OPM_THROW(std::runtime_error, "Pressure solver failed to converge.");
}
// Transport solve.
if (terminalOutputEnabled()) {
OpmLog::info("Solving the transport equations.");
}
const int transport_liniter = transport_solver_.step(dt, initial_reservoir_state_, initial_well_state_, reservoir_state, well_state);
if (transport_liniter == -1) {
OPM_THROW(std::runtime_error, "Transport solver failed to converge.");
}
// Report and return.
const bool converged = iteration >= 3; // TODO: replace this with a proper convergence check
return IterationReport { false, converged, pressure_liniter + transport_liniter };
}
}
/// Called once after each time step.
/// In this class, this function does nothing.
/// \param[in] dt time step size
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void afterStep(const double /* dt */,
ReservoirState& /* reservoir_state */,
WellState& /* well_state */)
{
}
/// \brief Set threshold pressures that prevent or reduce flow.
/// This prevents flow across faces if the potential
/// difference is less than the threshold. If the potential
/// difference is greater, the threshold value is subtracted
/// before calculating flow. This is treated symmetrically, so
/// flow is prevented or reduced in both directions equally.
/// \param[in] threshold_pressures_by_face array of size equal to the number of faces
/// of the grid passed in the constructor.
void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face)
{
pressure_solver_.model().setThresholdPressures(threshold_pressures_by_face);
transport_solver_.model().setThresholdPressures(threshold_pressures_by_face);
}
/// Return true if output to cout is wanted.
bool terminalOutputEnabled() const
{
return pressure_solver_.model().terminalOutputEnabled();
}
/// Return the relative change in variables relevant to this model.
double relativeChange(const SimulationDataContainer& previous,
const SimulationDataContainer& current ) const
{
// TODO: this is a quick stopgap implementation, and should be evaluated more carefully.
return std::max(pressure_solver_.model().relativeChange(previous, current),
transport_solver_.model().relativeChange(previous, current));
}
protected:
typedef BlackoilPressureModel<Grid, WellModel> PressureModel;
typedef BlackoilTransportModel<Grid, WellModel> TransportModel;
typedef NonlinearSolver<PressureModel> PressureSolver;
typedef NonlinearSolver<TransportModel> TransportSolver;
std::unique_ptr<PressureModel> pressure_model_;
std::unique_ptr<TransportModel> transport_model_;
PressureSolver pressure_solver_;
TransportSolver transport_solver_;
ReservoirState initial_reservoir_state_;
WellState initial_well_state_;
bool iterate_to_fully_implicit_;
};
} // namespace Opm
#endif // OPM_BLACKOILSEQUENTIALMODEL_HEADER_INCLUDED

638
opm/autodiff/BlackoilTransportModel.hpp Normal file
View File

@ -0,0 +1,638 @@
/*
Copyright 2015, 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil AS.
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_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED
#define OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED
#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
namespace Opm {
/// A model implementation for the transport equation in three-phase black oil.
template<class Grid, class WellModel>
class BlackoilTransportModel : public BlackoilModelBase<Grid, WellModel, BlackoilTransportModel<Grid, WellModel> >
{
public:
typedef BlackoilModelBase<Grid, WellModel, BlackoilTransportModel<Grid, WellModel> > Base;
friend Base;
typedef typename Base::ReservoirState ReservoirState;
typedef typename Base::WellState WellState;
typedef typename Base::SolutionState SolutionState;
typedef typename Base::V V;
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
/// \param[in] param parameters
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells_arg well structure
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
BlackoilTransportModel(const typename Base::ModelParameters& param,
const Grid& grid,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const StandardWells& std_wells,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
eclState, has_disgas, has_vapoil, terminal_output)
{
}
void prepareStep(const double dt,
const ReservoirState& reservoir_state,
const WellState& well_state)
{
Base::prepareStep(dt, reservoir_state, well_state);
Base::param_.solve_welleq_initially_ = false;
state0_ = variableState(reservoir_state, well_state);
asImpl().makeConstantState(state0_);
}
void assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly)
{
using namespace Opm::AutoDiffGrid;
// If we have VFP tables, we need the well connection
// pressures for the "simple" hydrostatic correction
// between well depth and vfp table depth.
if (isVFPActive()) {
SolutionState state = asImpl().variableState(reservoir_state, well_state);
SolutionState state0 = state;
asImpl().makeConstantState(state0);
asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
}
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
asImpl().wellModel().updateWellControls(terminal_output_, well_state);
// Create the primary variables.
SolutionState state = asImpl().variableState(reservoir_state, well_state);
if (initial_assembly) {
is_first_iter_ = true;
// Create the (constant, derivativeless) initial state.
SolutionState state0 = state;
asImpl().makeConstantState(state0);
// Compute initial accumulation contributions
// and well connection pressures.
asImpl().computeAccum(state0, 0);
asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
} else {
is_first_iter_ = false;
}
// -------- Mass balance equations --------
asImpl().assembleMassBalanceEq(state);
// -------- Well equations ----------
if ( ! wellsActive() ) {
return;
}
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
asImpl().wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
}
V aliveWells;
std::vector<ADB> cq_s;
// @afr changed
// asImpl().wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
// end of changed
asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
if (param_.compute_well_potentials_) {
SolutionState state0 = state;
asImpl().makeConstantState(state0);
asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
}
}
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
V solveJacobianSystem() const
{
const int n_transport = residual_.material_balance_eq[1].size();
const int n_full = residual_.sizeNonLinear();
const auto& mb = residual_.material_balance_eq;
LinearisedBlackoilResidual transport_res = {
{
// TODO: handle general 2-phase etc.
ADB::function(mb[1].value(), { mb[1].derivative()[1], mb[1].derivative()[2] }),
ADB::function(mb[2].value(), { mb[2].derivative()[1], mb[2].derivative()[2] })
},
ADB::null(),
ADB::null(),
residual_.matbalscale,
residual_.singlePrecision
};
assert(transport_res.sizeNonLinear() == 2*n_transport);
V dx_transport = linsolver_.computeNewtonIncrement(transport_res);
assert(dx_transport.size() == 2*n_transport);
V dx_full = V::Zero(n_full);
for (int i = 0; i < 2*n_transport; ++i) {
dx_full(n_transport + i) = dx_transport(i);
}
return dx_full;
}
using Base::numPhases;
using Base::numMaterials;
protected:
using Base::asImpl;
using Base::materialName;
using Base::convergenceReduction;
using Base::maxResidualAllowed;
using Base::linsolver_;
using Base::residual_;
using Base::rq_;
using Base::geo_;
using Base::ops_;
using Base::grid_;
using Base::use_threshold_pressure_;
using Base::canph_;
using Base::active_;
using Base::pvdt_;
using Base::fluid_;
using Base::param_;
using Base::terminal_output_;
using Base::isVFPActive;
using Base::phaseCondition;
using Base::vfp_properties_;
using Base::wellsActive;
V total_flux_; // HACK, should be part of a revised (transport-specific) SolutionState.
V total_wellperf_flux_;
DataBlock comp_wellperf_flux_;
SolutionState state0_ = SolutionState(3);
bool is_first_iter_ = false;
Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind_flags_;
SolutionState
variableState(const ReservoirState& x,
const WellState& xw) const
{
// HACK
const_cast<V&>(total_flux_)
= Eigen::Map<const V>(x.faceflux().data(), x.faceflux().size());
const_cast<V&>(total_wellperf_flux_)
= Eigen::Map<const V>(xw.perfRates().data(), xw.perfRates().size());
const_cast<DataBlock&>(comp_wellperf_flux_)
= Eigen::Map<const DataBlock>(xw.perfPhaseRates().data(), xw.perfRates().size(), numPhases());
assert(numPhases() * xw.perfRates().size() == xw.perfPhaseRates().size());
// As Base::variableState(), except making Pressure, Qs and Bhp constants.
std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
std::vector<ADB> vars = ADB::variables(vars0);
const std::vector<int> indices = asImpl().variableStateIndices();
vars[indices[Pressure]] = ADB::constant(vars[indices[Pressure]].value());
vars[indices[Qs]] = ADB::constant(vars[indices[Qs]].value());
vars[indices[Bhp]] = ADB::constant(vars[indices[Bhp]].value());
return asImpl().variableStateExtractVars(x, indices, vars);
}
void computeAccum(const SolutionState& state,
const int aix )
{
if (aix == 0) {
// The pressure passed in state is from after
// the pressure solver, but we need to use the original
// b factors etc. to get the initial accumulation term
// correct.
Base::computeAccum(state0_, aix);
} else {
Base::computeAccum(state, aix);
}
}
void assembleMassBalanceEq(const SolutionState& state)
{
// Compute b_p and the accumulation term b_p*s_p for each phase,
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
// These quantities are stored in rq_[phase].accum[1].
// The corresponding accumulation terms from the start of
// the timestep (b^0_p*s^0_p etc.) were already computed
// on the initial call to assemble() and stored in rq_[phase].accum[0].
asImpl().computeAccum(state, 1);
// Set up the common parts of the mass balance equations
// for each active phase.
const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
const V trans_nnc = ops_.nnc_trans;
V trans_all(transi.size() + trans_nnc.size());
trans_all << transi, trans_nnc;
const ADB tr_mult = asImpl().transMult(state.pressure);
const V gdz = geo_.gravity()[2] * (ops_.grad * geo_.z().matrix());
if (is_first_iter_) {
upwind_flags_.resize(gdz.size(), numPhases());
}
// Compute mobilities and heads
const std::vector<PhasePresence>& cond = asImpl().phaseCondition();
const std::vector<ADB> kr = asImpl().computeRelPerm(state);
#pragma omp parallel for schedule(static)
for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
// Compute and store mobilities.
const int canonical_phase_idx = canph_[ phase_idx ];
const ADB& phase_pressure = state.canonical_phase_pressures[canonical_phase_idx];
const ADB mu = asImpl().fluidViscosity(canonical_phase_idx, phase_pressure, state.temperature, state.rs, state.rv, cond);
// Note that the pressure-dependent transmissibility multipliers are considered
// part of the mobility here.
rq_[ phase_idx ].mob = tr_mult * kr[phase_idx] / mu;
// Compute head differentials. Gravity potential is done using the face average as in eclipse and MRST.
const ADB rho = asImpl().fluidDensity(canonical_phase_idx, rq_[phase_idx].b, state.rs, state.rv);
const ADB rhoavg = ops_.caver * rho;
rq_[ phase_idx ].dh = ops_.grad * phase_pressure - rhoavg * gdz;
if (is_first_iter_) {
upwind_flags_.col(phase_idx) = -rq_[phase_idx].dh.value();
}
if (use_threshold_pressure_) {
asImpl().applyThresholdPressures(rq_[ phase_idx ].dh);
}
}
// Extract saturation-dependent part of head differences.
const ADB gradp = ops_.grad * state.pressure;
std::vector<ADB> dh_sat(numPhases(), ADB::null());
for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
dh_sat[phase_idx] = gradp - rq_[phase_idx].dh;
}
// Find upstream directions for each phase.
upwind_flags_ = multiPhaseUpwind(dh_sat, trans_all);
// Compute (upstream) phase and total mobilities for connections.
// Also get upstream b, rs, and rv values to avoid recreating the UpwindSelector.
std::vector<ADB> mob(numPhases(), ADB::null());
std::vector<ADB> b(numPhases(), ADB::null());
ADB rs = ADB::null();
ADB rv = ADB::null();
ADB tot_mob = ADB::constant(V::Zero(gdz.size()));
for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
UpwindSelector<double> upwind(grid_, ops_, upwind_flags_.col(phase_idx));
mob[phase_idx] = upwind.select(rq_[phase_idx].mob);
tot_mob += mob[phase_idx];
b[phase_idx] = upwind.select(rq_[phase_idx].b);
if (canph_[phase_idx] == Oil) {
rs = upwind.select(state.rs);
}
if (canph_[phase_idx] == Gas) {
rv = upwind.select(state.rv);
}
}
// Compute phase fluxes.
for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
ADB gflux = ADB::constant(V::Zero(gdz.size()));
for (int other_phase = 0; other_phase < numPhases(); ++other_phase) {
if (phase_idx != other_phase) {
gflux += mob[other_phase] * (dh_sat[phase_idx] - dh_sat[other_phase]);
}
}
rq_[phase_idx].mflux = b[phase_idx] * (mob[phase_idx] / tot_mob) * (total_flux_ + trans_all * gflux);
}
#pragma omp parallel for schedule(static)
for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
// const int canonical_phase_idx = canph_[ phase_idx ];
// const ADB& phase_pressure = state.canonical_phase_pressures[canonical_phase_idx];
// asImpl().computeMassFlux(phase_idx, trans_all, kr[canonical_phase_idx], phase_pressure, state);
// Material balance equation for this phase.
residual_.material_balance_eq[ phase_idx ] =
pvdt_ * (rq_[phase_idx].accum[1] - rq_[phase_idx].accum[0])
+ ops_.div*rq_[phase_idx].mflux;
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
// Add the extra (flux) terms to the mass balance equations
// From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
// The extra terms in the accumulation part of the equation are already handled.
if (active_[ Oil ] && active_[ Gas ]) {
const int po = fluid_.phaseUsage().phase_pos[ Oil ];
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
residual_.material_balance_eq[ pg ] += ops_.div * (rs * rq_[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv * rq_[pg].mflux);
}
if (param_.update_equations_scaling_) {
asImpl().updateEquationsScaling();
}
}
Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> multiPhaseUpwind(const std::vector<ADB>& head_diff,
const V& transmissibility)
{
// Based on the paper "Upstream Differencing for Multiphase Flow in Reservoir Simulation",
// by Yann Brenier and Jérôme Jaffré,
// SIAM J. Numer. Anal., 28(3), 685–696.
// DOI:10.1137/0728036
// Using the data members:
// total_flux_
// rq_[].mob
// Notation based on paper cited above.
const int num_connections = head_diff[0].size();
Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind(num_connections, numPhases());
using ValueAndIndex = std::pair<double, int>;
const int num_phases = numPhases();
std::vector<ValueAndIndex> g(num_phases);
std::vector<double> theta(num_phases);
for (int conn = 0; conn < num_connections; ++conn) {
const double q = total_flux_[conn];
const double t = transmissibility[conn];
const int a = ops_.connection_cells(conn, 0); // first cell of connection
const int b = ops_.connection_cells(conn, 1); // second cell of connection
// Get and sort the g values (also called "weights" in the paper) for this connection.
for (int phase_idx = 0; phase_idx < num_phases; ++phase_idx) {
g[phase_idx] = ValueAndIndex(head_diff[phase_idx].value()[conn], phase_idx);
}
std::sort(g.begin(), g.end());
// Compute theta and r.
// Paper notation: subscript l -> ell (for read/searchability)
// Note that since we index phases from 0, r is one less than in the paper.
int r = -1;
for (int ell = 0; ell < num_phases; ++ell) {
theta[ell] = q;
for (int j = 0; j < num_phases; ++j) {
if (j < ell) {
theta[ell] += t * (g[ell].first - g[j].first) * rq_[g[j].second].mob.value()[b];
}
if (j > ell) {
theta[ell] += t * (g[ell].first - g[j].first) * rq_[g[j].second].mob.value()[a];
}
}
if (theta[ell] <= 0.0) {
r = ell;
} else {
break; // r is correct, no need to continue
}
}
for (int ell = 0; ell < num_phases; ++ell) {
const int phase_idx = g[ell].second;
upwind(conn, phase_idx) = ell > r ? 1.0 : -1.0;
}
}
return upwind;
}
void computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& /* aliveWells */,
std::vector<ADB>& cq_s) const
{
// Note that use of this function replaces using the well models'
// function of the same name.
if( ! asImpl().localWellsActive() ) return ;
const int np = asImpl().wells().number_of_phases;
const int nw = asImpl().wells().number_of_wells;
const int nperf = asImpl().wells().well_connpos[nw];
const Opm::PhaseUsage& pu = asImpl().fluid_.phaseUsage();
// Compute total mobilities for perforations.
ADB totmob_perfcells = ADB::constant(V::Zero(nperf));
for (int phase = 0; phase < numPhases(); ++phase) {
totmob_perfcells += mob_perfcells[phase];
}
// Compute fractional flow.
std::vector<ADB> frac_flow(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
frac_flow[phase] = mob_perfcells[phase] / totmob_perfcells;
}
// Identify injecting and producing perforations.
V is_inj = V::Zero(nperf);
V is_prod = V::Zero(nperf);
for (int c = 0; c < nperf; ++c){
if (total_wellperf_flux_[c] > 0.0) {
is_inj[c] = 1;
} else {
is_prod[c] = 1;
}
}
// Compute fluxes for producing perforations.
std::vector<ADB> cq_s_prod(3, ADB::null());
for (int phase = 0; phase < np; ++phase) {
// For producers, we use the total reservoir flux from the pressure solver.
cq_s_prod[phase] = b_perfcells[phase] * frac_flow[phase] * total_wellperf_flux_;
}
if (asImpl().has_disgas_ || asImpl().has_vapoil_) {
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
const ADB cq_s_prod_oil = cq_s_prod[oilpos];
const ADB cq_s_prod_gas = cq_s_prod[gaspos];
cq_s_prod[gaspos] += subset(state.rs, Base::well_model_.wellOps().well_cells) * cq_s_prod_oil;
cq_s_prod[oilpos] += subset(state.rv, Base::well_model_.wellOps().well_cells) * cq_s_prod_gas;
}
// Compute well perforation surface volume fluxes.
cq_s.resize(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase];
// For injectors, we use the component fluxes computed by the pressure solver.
const V cq_s_inj = comp_wellperf_flux_.col(pos);
cq_s[phase] = is_prod * cq_s_prod[phase] + is_inj * cq_s_inj;
}
}
bool getConvergence(const double dt, const int iteration)
{
const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = param_.tolerance_cnv_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = asImpl().numPhases();
const int nm = asImpl().numMaterials();
assert(int(rq_.size()) == nm);
const V& pv = geo_.poreVolume();
std::vector<double> R_sum(nm);
std::vector<double> B_avg(nm);
std::vector<double> maxCoeff(nm);
std::vector<double> maxNormWell(np);
Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
Eigen::Array<typename V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
for ( int idx = 0; idx < nm; ++idx )
{
const ADB& tempB = rq_[idx].b;
B.col(idx) = 1./tempB.value();
R.col(idx) = residual_.material_balance_eq[idx].value();
tempV.col(idx) = R.col(idx).abs()/pv;
}
const double pvSum = convergenceReduction(B, tempV, R,
R_sum, maxCoeff, B_avg, maxNormWell,
nc);
std::vector<double> CNV(nm);
std::vector<double> mass_balance_residual(nm);
std::vector<double> well_flux_residual(np);
bool converged_MB = true;
bool converged_CNV = true;
// Finish computation
for ( int idx = 1; idx < nm; ++idx ) {
CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
assert(nm >= np);
}
const bool converged = converged_MB && converged_CNV;
for (int idx = 0; idx < nm; ++idx) {
if (std::isnan(mass_balance_residual[idx])
|| std::isnan(CNV[idx])
|| (idx < np && std::isnan(well_flux_residual[idx]))) {
OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx));
}
if (mass_balance_residual[idx] > maxResidualAllowed()
|| CNV[idx] > maxResidualAllowed()
|| (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx));
}
}
if ( terminal_output_ ) {
// Only rank 0 does print to std::cout
std::ostringstream os;
if (iteration == 0) {
os << "\nIter";
for (int idx = 1; idx < nm; ++idx) {
os << " MB(" << materialName(idx).substr(0, 3) << ") ";
}
for (int idx = 1; idx < nm; ++idx) {
os << " CNV(" << materialName(idx).substr(0, 1) << ") ";
}
os << '\n';
}
os.precision(3);
os.setf(std::ios::scientific);
os << std::setw(4) << iteration;
for (int idx = 1; idx < nm; ++idx) {
os << std::setw(11) << mass_balance_residual[idx];
}
for (int idx = 1; idx < nm; ++idx) {
os << std::setw(11) << CNV[idx];
}
OpmLog::info(os.str());
}
return converged;
}
};
/// Providing types by template specialisation of ModelTraits for BlackoilTransportModel.
template <class Grid, class WellModel>
struct ModelTraits< BlackoilTransportModel<Grid, WellModel> >
{
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters;
typedef DefaultBlackoilSolutionState SolutionState;
};
} // namespace Opm
#endif // OPM_BLACKOILTRANSPORTMODEL_HEADER_INCLUDED

141
opm/autodiff/FlowMainSequential.hpp Normal file
View File

@ -0,0 +1,141 @@
/*
Copyright 2016 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_FLOWMAINSEQUENTIAL_HEADER_INCLUDED
#define OPM_FLOWMAINSEQUENTIAL_HEADER_INCLUDED
#include <opm/autodiff/FlowMain.hpp>
namespace Opm
{
// The FlowMainSequential class is for a black-oil simulator using the sequential models.
template <class Grid, class Simulator>
class FlowMainSequential : public FlowMainBase<FlowMainSequential<Grid, Simulator>, Grid, Simulator>
{
protected:
using Base = FlowMainBase<FlowMainSequential<Grid, Simulator>, Grid, Simulator>;
using Base::eclipse_state_;
using Base::param_;
using Base::fis_solver_;
using Base::parallel_information_;
friend Base;
// ------------ Methods ------------
// Print startup message if on output rank.
void printStartupMessage()
{
if (Base::output_cout_) {
const std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow-Sequential (version " << version << ")"
<< std::string(17 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow-Sequential is a simulator for fully implicit three-phase, *\n";
std::cout << "* black-oil flow, and is part of OPM. *\n";
std::cout << "* For more information see http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
}
// Setup linear solver.
// Writes to:
// fis_solver_
// param_ (conditionally)
// The CPR solver cannot be used with the sequential model.
// Also, the interleaved solver requires the full sparsity pattern option.
void setupLinearSolver()
{
const std::string cprSolver = "cpr";
const std::string interleavedSolver = "interleaved";
const std::string directSolver = "direct";
std::string flowDefaultSolver = interleavedSolver;
if (!param_.has("solver_approach")) {
if (eclipse_state_->getSimulationConfig()->useCPR()) {
flowDefaultSolver = cprSolver;
}
}
const std::string solver_approach = param_.getDefault("solver_approach", flowDefaultSolver);
if (solver_approach == cprSolver) {
OPM_THROW( std::runtime_error , "CPR solver is not ready for use with sequential simulator.");
} else if (solver_approach == interleavedSolver) {
if (!param_.has("require_full_sparsity_pattern")) {
param_.insertParameter("require_full_sparsity_pattern", "true");
}
fis_solver_.reset(new NewtonIterationBlackoilInterleaved(param_, parallel_information_));
} else if (solver_approach == directSolver) {
fis_solver_.reset(new NewtonIterationBlackoilSimple(param_, parallel_information_));
} else {
OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized.");
}
}
// Create simulator instance.
// Writes to:
// simulator_
void createSimulator()
{
// We must override the min_iter argument unless it was already supplied, to avoid requiring iteration.
if (!param_.has("min_iter")) {
param_.insertParameter("min_iter", "0");
}
// Create the simulator instance.
Base::simulator_.reset(new Simulator(Base::param_,
Base::grid_init_->grid(),
*Base::geoprops_,
*Base::fluidprops_,
Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr,
*Base::fis_solver_,
Base::gravity_.data(),
Base::deck_->hasKeyword("DISGAS"),
Base::deck_->hasKeyword("VAPOIL"),
Base::eclipse_state_,
*Base::output_writer_,
Base::threshold_pressures_));
}
};
} // namespace Opm
#endif // OPM_FLOWMAINSEQUENTIAL_HEADER_INCLUDED

116
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@ -567,6 +567,116 @@ namespace Opm
}
};
std::pair<NewtonIterationBlackoilInterleaved::SolutionVector, Dune::InverseOperatorResult>
computePressureIncrement(const LinearisedBlackoilResidual& residual)
{
typedef LinearisedBlackoilResidual::ADB ADB;
typedef ADB::V V;
// Build the vector of equations (should be just a single material balance equation
// in which the pressure equation is stored).
const int np = residual.material_balance_eq.size();
assert(np == 1);
std::vector<ADB> eqs;
eqs.reserve(np + 2);
for (int phase = 0; phase < np; ++phase) {
eqs.push_back(residual.material_balance_eq[phase]);
}
// Check if wells are present.
const bool hasWells = residual.well_flux_eq.size() > 0 ;
std::vector<ADB> elim_eqs;
if (hasWells) {
// Eliminate the well-related unknowns, and corresponding equations.
eqs.push_back(residual.well_flux_eq);
eqs.push_back(residual.well_eq);
elim_eqs.reserve(2);
elim_eqs.push_back(eqs[np]);
eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns.
elim_eqs.push_back(eqs[np]);
eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns.
assert(int(eqs.size()) == np);
}
// Solve the linearised oil equation.
Eigen::SparseMatrix<double, Eigen::RowMajor> eigenA = eqs[0].derivative()[0].getSparse();
DuneMatrix opA(eigenA);
const int size = eqs[0].size();
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > Vector1;
Vector1 x;
x.resize(size);
x = 0.0;
Vector1 b;
b.resize(size);
b = 0.0;
std::copy_n(eqs[0].value().data(), size, b.begin());
// Solve with AMG solver.
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1> > Mat;
typedef Dune::MatrixAdapter<Mat, Vector1, Vector1> Operator;
Operator sOpA(opA);
typedef Dune::Amg::SequentialInformation ParallelInformation;
typedef Dune::SeqILU0<Mat,Vector1,Vector1> EllipticPreconditioner;
typedef EllipticPreconditioner Smoother;
typedef Dune::Amg::AMG<Operator, Vector1, Smoother, ParallelInformation> AMG;
typedef Dune::Amg::FirstDiagonal CouplingMetric;
typedef Dune::Amg::SymmetricCriterion<Mat, CouplingMetric> CritBase;
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
// TODO: revise choice of parameters
const int coarsenTarget = 1200;
Criterion criterion(15, coarsenTarget);
criterion.setDebugLevel(0); // no debug information, 1 for printing hierarchy information
criterion.setDefaultValuesIsotropic(2);
criterion.setNoPostSmoothSteps(1);
criterion.setNoPreSmoothSteps(1);
// for DUNE 2.2 we also need to pass the smoother args
typedef typename AMG::Smoother Smoother;
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
SmootherArgs smootherArgs;
smootherArgs.iterations = 1;
smootherArgs.relaxationFactor = 1.0;
AMG precond(sOpA, criterion, smootherArgs);
const int verbosity = 0;
const int maxit = 30;
const double tolerance = 1e-5;
// Construct linear solver.
Dune::BiCGSTABSolver<Vector1> linsolve(sOpA, precond, tolerance, maxit, verbosity);
// Solve system.
Dune::InverseOperatorResult result;
linsolve.apply(x, b, result);
// Check for failure of linear solver.
if (!result.converged) {
OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver in computePressureIncrement().");
}
// Copy solver output to dx.
NewtonIterationBlackoilInterleaved::SolutionVector dx(size);
for (int i = 0; i < size; ++i) {
dx(i) = x[i];
}
if (hasWells) {
// Compute full solution using the eliminated equations.
// Recovery in inverse order of elimination.
dx = recoverVariable(elim_eqs[1], dx, np);
dx = recoverVariable(elim_eqs[0], dx, np);
}
return std::make_pair(dx, result);
}
} // end namespace detail
@ -575,6 +685,12 @@ namespace Opm
{
// get np and call appropriate template method
const int np = residual.material_balance_eq.size();
if (np == 1) {
auto result = detail::computePressureIncrement(residual);
iterations_ = result.second.iterations;
return result.first;
}
const NewtonIterationBlackoilInterface& newtonIncrement = residual.singlePrecision ?
detail::NewtonIncrement< maxNumberEquations_, float > :: get( newtonIncrementSinglePrecision_, parameters_, parallelInformation_, np ) :
detail::NewtonIncrement< maxNumberEquations_, double > :: get( newtonIncrementDoublePrecision_, parameters_, parallelInformation_, np );

74
opm/autodiff/SimulatorSequentialBlackoil.hpp Normal file
View File

@ -0,0 +1,74 @@
/*
Copyright 2013, 2015, 2016 SINTEF ICT, Applied Mathematics.
Copyright 2015 Andreas Lauser
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_SIMULATORSEQUENTIALBLACKOIL_HEADER_INCLUDED
#define OPM_SIMULATORSEQUENTIALBLACKOIL_HEADER_INCLUDED
#include <opm/autodiff/SimulatorBase.hpp>
#include <opm/autodiff/NonlinearSolver.hpp>
#include <opm/autodiff/BlackoilSequentialModel.hpp>
#include <opm/autodiff/StandardWells.hpp>
namespace Opm {
template <class GridT>
class SimulatorSequentialBlackoil;
class StandardWells;
template <class GridT>
struct SimulatorTraits<SimulatorSequentialBlackoil<GridT> >
{
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilState ReservoirState;
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef BlackoilSequentialModel<Grid, StandardWells> Model;
typedef NonlinearSolver<Model> Solver;
typedef StandardWells WellModel;
};
/// a simulator for the blackoil model
template <class GridT>
class SimulatorSequentialBlackoil
: public SimulatorBase<SimulatorSequentialBlackoil<GridT> >
{
typedef SimulatorBase<SimulatorSequentialBlackoil<GridT> > Base;
public:
// forward the constructor to the base class
SimulatorSequentialBlackoil(const parameter::ParameterGroup& param,
const typename Base::Grid& grid,
DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool disgas,
const bool vapoil,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
const std::vector<double>& threshold_pressures_by_face)
: Base(param, grid, geo, props, rock_comp_props, linsolver, gravity, disgas, vapoil,
eclipse_state, output_writer, threshold_pressures_by_face)
{}
};
} // namespace Opm
#endif // OPM_SIMULATORSEQUENTIALBLACKOIL_HEADER_INCLUDED