mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-01-13 18:22:02 -06:00
Add sequential models for black oil.
This commit adds sequential solvers, including a simulator variant using them (flow_sequential.cpp) with an integration test (running SPE1, same as for fully implicit). The sequential code is capable of running several (but not all) test cases without tuning or special parameters, but reducing ds_max a bit (from default 0.2 to say 0.1) helps with transport solver convergence. The Norne model runs fine (esp. with a little tuning). A parameter iterate_to_fully_implicit (defaults to false) is available, when set the simulator will iterate with alternating pressure and transport solves towards the fully implicit solution. Although that takes a lot extra time it serves as a correctness check. Performance is not competitive with fully implicit at this point: essentially both the pressure and transport models inherit the fully implicit model and do a lot of double (or triple) work. The point has been to establish a proof of concept and baseline for further experiments, without disturbing the base model too much (or at all, if possible). Changes to existing code has been minimized by merging most such changes as smaller PRs already, the only remaining such change is to NewtonIterationBlackoilInterleaved. Admittedly, that code (to solve the pressure system with AMG) is not ideal because it duplicates similar code in CPRPreconditioner.hpp and is not parallel. I propose to address this later by refactoring the "solve elliptic system" code from CPRPreconditioner into a separate class that can be used also from here
This commit is contained in:
parent
82d0c81c13
commit
7489f15ee9
@ -146,5 +146,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()
|
||||
|
@ -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
39
examples/flow_sequential.cpp
Normal 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
330
opm/autodiff/BlackoilPressureModel.hpp
Normal 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
267
opm/autodiff/BlackoilSequentialModel.hpp
Normal 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
638
opm/autodiff/BlackoilTransportModel.hpp
Normal 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
141
opm/autodiff/FlowMainSequential.hpp
Normal 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
|
@ -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
74
opm/autodiff/SimulatorSequentialBlackoil.hpp
Normal 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
|
Loading…
Reference in New Issue
Block a user