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

Use new linear solver interface in fully implicit solver.

Browse Source 
This includes:
 - Using the class FullyImplicitBlackoilResidual instead of
   in-class definition for the residual object.
 - Changing residual field name mass_balance to material_balance_eq.
 - Letting the simulator and solver classes accept a
   FullyImplicitSystemSolverInterface instead of a LinearSolverInterface.
 - In sim_fibo_ad and test_implicit_ad, instantiate class
   FullyImplicitSystemSolverSimple, replicating existing behaviour.
This commit is contained in:
Atgeirr Flø Rasmussen 2014-02-10 10:47:21 +01:00
parent 8850868f50
commit a14cff7834
6 changed files with 30 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
examples/sim_fibo_ad.cpp
View File

@ -41,6 +41,7 @@
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/FullyImplicitSystemSolverSimple.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
@ -137,6 +138,7 @@ try
// Linear solver.
LinearSolverFactory linsolver(param);
FullyImplicitSystemSolverSimple fis_solver(linsolver);
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
@ -211,7 +213,7 @@ try
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
linsolver,
fis_solver,
grav,
outputWriter);
if (epoch == 0) {

4
examples/test_implicit_ad.cpp
View File

@ -23,6 +23,7 @@
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/FullyImplicitSystemSolverSimple.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
@ -103,8 +104,9 @@ try
Opm::DerivedGeology geo(*g, props, grav);
Opm::LinearSolverFactory linsolver(param);
Opm::FullyImplicitSystemSolverSimple fis_solver(linsolver);
Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, linsolver);
Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, fis_solver);
Opm::BlackoilState state;
initStateBasic(*g, props0, param, 0.0, state);

57
opm/autodiff/FullyImplicitBlackoilSolver.cpp
View File

@ -194,7 +194,7 @@ namespace {
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const Wells& wells,
const LinearSolverInterface& linsolver)
const FullyImplicitSystemSolverInterface& linsolver)
: grid_ (grid)
, fluid_ (fluid)
, geo_ (geo)
@ -643,13 +643,13 @@ namespace {
// std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl;
// std::cout << rq_[phase].mflux;
residual_.mass_balance[ phaseIdx ] =
residual_.material_balance_eq[ phaseIdx ] =
pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ ops_.div*rq_[phaseIdx].mflux;
// DUMP(ops_.div*rq_[phase].mflux);
// DUMP(residual_.mass_balance[phase]);
// DUMP(residual_.material_balance_eq[phase]);
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
@ -663,16 +663,16 @@ namespace {
rq_[po].head.value());
const ADB rs_face = upwindOil.select(state.rs);
residual_.mass_balance[ Gas ] += ops_.div * (rs_face * rq_[po].mflux);
residual_.material_balance_eq[ Gas ] += ops_.div * (rs_face * rq_[po].mflux);
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
const UpwindSelector<double> upwindGas(grid_, ops_,
rq_[pg].head.value());
const ADB rv_face = upwindGas.select(state.rv);
residual_.mass_balance[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux);
residual_.material_balance_eq[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux);
// DUMP(residual_.mass_balance[ Gas ]);
// DUMP(residual_.material_balance_eq[ Gas ]);
}
@ -772,7 +772,7 @@ namespace {
// const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc);
// DUMP(well_contribs[phase]);
residual_.mass_balance[phase] += well_contribs[phase];
residual_.material_balance_eq[phase] += well_contribs[phase];
}
if (active_[Gas] && active_[Oil]) {
const int oilpos = pu.phase_pos[Oil];
@ -782,8 +782,8 @@ namespace {
well_rates_all += superset(wops_.p2w * (well_perf_rates[oilpos]*rs_perf), Span(nw, 1, gaspos*nw), nw*np);
well_rates_all += superset(wops_.p2w * (well_perf_rates[gaspos]*rv_perf), Span(nw, 1, oilpos*nw), nw*np);
// DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs);
residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs;
residual_.mass_balance[oilpos] += well_contribs[gaspos]*state.rv;
residual_.material_balance_eq[gaspos] += well_contribs[oilpos]*state.rs;
residual_.material_balance_eq[oilpos] += well_contribs[gaspos]*state.rv;
}
// Set the well flux equation
@ -826,40 +826,7 @@ namespace {
V FullyImplicitBlackoilSolver::solveJacobianSystem() const
{
const int np = fluid_.numPhases();
ADB mass_res = residual_.mass_balance[0];
for (int phase = 1; phase < np; ++phase) {
mass_res = vertcat(mass_res, residual_.mass_balance[phase]);
}
const ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq);
const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));
// DUMP(total_residual);
const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual.derivative()[0];
V dx(V::Zero(total_residual.size()));
Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(matr.rows(), matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
total_residual.value().data(), dx.data());
/*
std::ofstream filestream("matrix.out");
filestream << matr;
filestream.close();
std::ofstream filestream2("sol.out");
filestream2 << dx;
filestream2.close();
std::ofstream filestream3("r.out");
filestream3 << total_residual.value();
filestream3.close(); */
if (!rep.converged) {
OPM_THROW(std::runtime_error,
"FullyImplicitBlackoilSolver::solveJacobianSystem(): "
"Linear solver convergence failure.");
}
return dx;
return linsolver_.linearSolve(residual_);
}
@ -1268,8 +1235,8 @@ namespace {
{
double r = 0;
for (std::vector<ADB>::const_iterator
b = residual_.mass_balance.begin(),
e = residual_.mass_balance.end();
b = residual_.material_balance_eq.begin(),
e = residual_.material_balance_eq.end();
b != e; ++b)
{
r = std::max(r, (*b).value().matrix().norm());

48
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@ -23,6 +23,8 @@
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/FullyImplicitBlackoilResidual.hpp>
#include <opm/autodiff/FullyImplicitSystemSolverInterface.hpp>
struct UnstructuredGrid;
struct Wells;
@ -31,7 +33,7 @@ namespace Opm {
class DerivedGeology;
class RockCompressibility;
class LinearSolverInterface;
class FullyImplicitSystemSolverInterface;
class BlackoilState;
class WellState;
@ -62,7 +64,7 @@ namespace Opm {
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const Wells& wells,
const LinearSolverInterface& linsolver);
const FullyImplicitSystemSolverInterface& linsolver);
/// Take a single forward step, modifiying
/// state.pressure()
@ -78,45 +80,9 @@ namespace Opm {
BlackoilState& state ,
WellState& wstate);
/// Typedef used throughout the solver, in the public section
/// since it is used in the definition of Residual.
typedef AutoDiffBlock<double> ADB;
/// Residual structure of the fully implicit solver.
/// All equations are given as AD types, with multiple
/// jacobian blocks corresponding to the primary unknowns. The
/// primary unknowns are for a three-phase simulation, in order:
/// p (pressure)
/// sw (water saturation)
/// xvar (gas saturation, gas-oil ratio or oil-gas ratio)
/// qs (well outflows by well and phase)
/// bhp (bottom hole pressures)
/// In the above, the xvar variable will have a different
/// meaning from cell to cell, corresponding to the state in
/// that cell (saturated, undersaturated oil or undersaturated
/// gas). In a two-phase simulation, either sw or xvar is not
/// used, depending on which face is missing.
struct Residual {
/// The mass_balance vector has one element for each
/// active phase, each of which has size equal to the
/// number of cells. Each mass balance equation is given
/// in terms of surface volumes.
std::vector<ADB> mass_balance;
/// The well_flux_eq has size equal to the number of wells
/// times the number of phases. It contains the well flow
/// equations, relating the total well flows to
/// bottom-hole pressures and reservoir conditions.
ADB well_flux_eq;
/// The well_eq has size equal to the number of wells. It
/// contains the well control equations, that is for each
/// well either a rate specification or bottom hole
/// pressure specification.
ADB well_eq;
};
private:
// Types and enums
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
@ -159,7 +125,7 @@ namespace Opm {
const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
const Wells& wells_;
const LinearSolverInterface& linsolver_;
const FullyImplicitSystemSolverInterface& linsolver_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
// Size = # active faces. Maps active -> canonical phase indices.
@ -172,7 +138,7 @@ namespace Opm {
std::vector<ReservoirResidualQuant> rq_;
std::vector<PhasePresence> phaseCondition_;
Residual residual_;
FullyImplicitBlackoilResidual residual_;
// Private methods.
SolutionState

6
opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
View File

@ -71,7 +71,7 @@ namespace Opm
const BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
FullyImplicitSystemSolverInterface& linsolver,
const double* gravity,
EclipseWriter &writer);
@ -113,7 +113,7 @@ namespace Opm
const BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
FullyImplicitSystemSolverInterface& linsolver,
const double* gravity,
EclipseWriter &eclipseWriter)
@ -261,7 +261,7 @@ namespace Opm
const BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
FullyImplicitSystemSolverInterface& linsolver,
const double* gravity,
EclipseWriter &eclipseWriter)
: grid_(grid),

4
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -34,7 +34,7 @@ namespace Opm
class EclipseWriter;
class RockCompressibility;
class WellsManager;
class LinearSolverInterface;
class FullyImplicitSystemSolverInterface;
class SimulatorTimer;
class BlackoilState;
class WellState;
@ -71,7 +71,7 @@ namespace Opm
const BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
FullyImplicitSystemSolverInterface& linsolver,
const double* gravity,
EclipseWriter &writer);