OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #954 from joakim-hove/set-component

Browse Source 
Changes in SimulatorState:
This commit is contained in:
Bård Skaflestad
2016-02-24 07:11:36 -06:00
parent 089229adf4 df6e6a2347
commit 893e6c851f
12 changed files with 165 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
opm/core/pressure/IncompTpfa.cpp
View File

@@ -28,7 +28,7 @@
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp>
@@ -155,7 +155,7 @@ namespace Opm
/// May throw an exception if the number of iterations
/// exceed maxiter (set in constructor).
void IncompTpfa::solve(const double dt,
TwophaseState& state,
SimulatorState& state,
WellState& well_state)
{
if (rock_comp_props_ != 0 && rock_comp_props_->isActive()) {
@@ -169,7 +169,7 @@ namespace Opm
// Solve with no rock compressibility (linear eqn).
void IncompTpfa::solveIncomp(const double dt,
TwophaseState& state,
SimulatorState& state,
WellState& well_state)
{
// Set up properties.
@@ -207,7 +207,7 @@ namespace Opm
// Solve with rock compressibility (nonlinear eqn).
void IncompTpfa::solveRockComp(const double dt,
TwophaseState& state,
SimulatorState& state,
WellState& well_state)
{
// This function is identical to CompressibleTpfa::solve().
@@ -321,7 +321,7 @@ namespace Opm
/// Compute per-solve dynamic properties.
void IncompTpfa::computePerSolveDynamicData(const double /*dt*/,
const TwophaseState& state,
const SimulatorState& state,
const WellState& /*well_state*/)
{
// Computed here:
@@ -369,7 +369,7 @@ namespace Opm
/// Compute per-iteration dynamic properties.
void IncompTpfa::computePerIterationDynamicData(const double /*dt*/,
const TwophaseState& state,
const SimulatorState& state,
const WellState& well_state)
{
// These are the variables that get computed by this function:
@@ -396,7 +396,7 @@ namespace Opm
/// Compute the residual in h_->b and Jacobian in h_->A.
void IncompTpfa::assemble(const double dt,
const TwophaseState& state,
const SimulatorState& state,
const WellState& /*well_state*/)
{
const double* pressures = wells_ ? &pressures_[0] : &state.pressure()[0];
@@ -462,7 +462,7 @@ namespace Opm
/// Compute the output.
void IncompTpfa::computeResults(TwophaseState& state,
void IncompTpfa::computeResults(SimulatorState& state,
WellState& well_state) const
{
// Make sure h_ contains the direct-solution matrix

18
opm/core/pressure/IncompTpfa.hpp
View File

@@ -20,6 +20,7 @@
#ifndef OPM_INCOMPTPFA_HEADER_INCLUDED
#define OPM_INCOMPTPFA_HEADER_INCLUDED
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <vector>
@@ -34,8 +35,9 @@ namespace Opm
class IncompPropertiesInterface;
class RockCompressibility;
class LinearSolverInterface;
class TwophaseState;
class WellState;
class SimulatoreState;
/// Encapsulating a tpfa pressure solver for the incompressible-fluid case.
/// Supports gravity, wells controlled by bhp or reservoir rates,
@@ -112,7 +114,7 @@ namespace Opm
/// May throw an exception if the number of iterations
/// exceed maxiter (set in constructor).
void solve(const double dt,
TwophaseState& state,
SimulatorState& state,
WellState& well_state);
@@ -122,28 +124,28 @@ namespace Opm
protected:
// Solve with no rock compressibility (linear eqn).
void solveIncomp(const double dt,
TwophaseState& state,
SimulatorState& state,
WellState& well_state);
// Solve with rock compressibility (nonlinear eqn).
void solveRockComp(const double dt,
TwophaseState& state,
SimulatorState& state,
WellState& well_state);
private:
// Helper functions.
void computeStaticData();
virtual void computePerSolveDynamicData(const double dt,
const TwophaseState& state,
const SimulatorState& state,
const WellState& well_state);
void computePerIterationDynamicData(const double dt,
const TwophaseState& state,
const SimulatorState& state,
const WellState& well_state);
void assemble(const double dt,
const TwophaseState& state,
const SimulatorState& state,
const WellState& well_state);
void solveIncrement();
double residualNorm() const;
double incrementNorm() const;
void computeResults(TwophaseState& state,
void computeResults(SimulatorState& state,
WellState& well_state) const;
protected:

10
opm/core/simulator/BlackoilState.cpp
View File

@@ -23,16 +23,6 @@ BlackoilState::init(const UnstructuredGrid& g, int num_phases)
{
init(g.number_of_cells, g.number_of_faces, num_phases);
}
/// Set the first saturation to either its min or max value in
/// the indicated cells. The second saturation value s2 is set
/// to (1.0 - s1) for each cell. Any further saturation values
/// are unchanged.
void
BlackoilState::setFirstSat(const std::vector<int>& cells,
const Opm::BlackoilPropertiesInterface& props,
ExtremalSat es) {
SimulatorState::setFirstSat(cells, props, es);
}
bool
BlackoilState::equals(const SimulatorState& other,

8
opm/core/simulator/BlackoilState.hpp
View File

@@ -39,14 +39,6 @@ namespace Opm
virtual void init(int number_of_cells, int number_of_faces, int num_phases);
/// Set the first saturation to either its min or max value in
/// the indicated cells. The second saturation value s2 is set
/// to (1.0 - s1) for each cell. Any further saturation values
/// are unchanged.
void setFirstSat(const std::vector<int>& cells,
const Opm::BlackoilPropertiesInterface& props,
ExtremalSat es);
virtual bool equals(const SimulatorState& other,
double epsilon = 1e-8) const;

77
opm/core/simulator/SimulatorState.cpp
View File

@@ -2,6 +2,7 @@
#include <opm/common/util/numeric/cmp.hpp>
#include <opm/core/simulator/SimulatorState.hpp>
#include <algorithm>
#include <cmath>
#include <cassert>
@@ -80,42 +81,50 @@ SimulatorState::registerFaceData( const std::string& name, const int components,
return pos ;
}
template <typename Props> void
SimulatorState::setFirstSat(const std::vector<int>& cells,
const Props& props,
ExtremalSat es)
{
if (cells.empty()) {
return;
}
int n = cells.size();
std::vector<double> smin(num_phases_*n);
std::vector<double> smax(num_phases_*n);
props.satRange(n, &cells[0], &smin[0], &smax[0]);
std::vector< double >& sat = saturation();
const double* svals = (es == MinSat) ? &smin[0] : &smax[0];
for (int ci = 0; ci < n; ++ci) {
const int cell = cells[ci];
sat[num_phases_*cell] = svals[num_phases_*ci];
sat[num_phases_*cell + 1] = 1.0 - sat[num_phases_*cell];
void SimulatorState::setCellDataComponent( const std::string& name , size_t component , const std::vector<int>& cells , const std::vector<double>& values) {
const auto iter = std::find( cellDataNames_.begin() , cellDataNames_.end() , name);
int id = iter - cellDataNames_.begin();
auto& data = cellData_[id];
if (component >= num_phases_)
throw std::invalid_argument("Invalid component");
if (cells.size() != values.size())
throw std::invalid_argument("size mismatch between cells and values");
/* This is currently quite broken; the setCellDataComponent
method assumes that the number of components in the field
we are currently focusing on has num_phases components in
total. This restriction should be lifted by allowing a per
field number of components.
*/
if (data.size() != num_phases_ * num_cells_)
throw std::invalid_argument("Can currently only be used on fields with num_components == num_phases (i.e. saturation...) ");
for (size_t i = 0; i < cells.size(); i++) {
if (cells[i] < num_cells_) {
auto field_index = cells[i] * num_phases_ + component;
auto value = values[i];
data[field_index] = value;
} else {
throw std::invalid_argument("Invalid cell number");
}
}
}
// template instantiations for all known (to this library) subclasses
// of SimulatorState that will call this method. notice that there are
// no empty angle brackets after "template" -- that would have been
// specialization instead
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
template void
SimulatorState::setFirstSat <IncompPropertiesInterface> (
const std::vector<int> &cells,
const IncompPropertiesInterface &props,
ExtremalSat es);
std::vector<double>& SimulatorState::getCellData( const std::string& name ) {
const auto iter = std::find( cellDataNames_.begin() , cellDataNames_.end() , name);
int id = iter - cellDataNames_.begin();
auto& data = cellData_[id];
return data;
}
const std::vector<double>& SimulatorState::getCellData( const std::string& name ) const {
const auto iter = std::find( cellDataNames_.begin() , cellDataNames_.end() , name);
int id = iter - cellDataNames_.begin();
const auto& data = cellData_[id];
return data;
}
template void
SimulatorState::setFirstSat <BlackoilPropertiesInterface> (
const std::vector<int> &cells,
const BlackoilPropertiesInterface &props,
ExtremalSat es);

22
opm/core/simulator/SimulatorState.hpp
View File

@@ -17,8 +17,6 @@ namespace Opm
virtual void init(int number_of_cells, int number_of_faces, int num_phases);
enum ExtremalSat { MinSat, MaxSat };
protected:
/// \brief pressure per cell.
static const int pressureId_ = 0;
@@ -32,19 +30,12 @@ namespace Opm
/// \brief The fluxes at the faces.
static const int faceFluxId_ = 1;
/**
* Initialize the first saturation to maximum value. This method
* should be considered deprecated. Avoid to use it!
*
* \tparam Props Fluid and rock properties that pertain to this
* kind of simulation. Currently, only Blackoil-
* and IncompPropertiesInterface are supported.
*/
template <typename Props>
void setFirstSat(const std::vector<int>& cells,
const Props& props,
ExtremalSat es);
public:
/// Will set the values of component nr @component in the
/// field @key. All the cells in @cells will be set to the
/// values in @values.
void setCellDataComponent( const std::string& key , size_t component , const std::vector<int>& cells , const std::vector<double>& values);
int numPhases() const { return num_phases_; }
int numCells () const { return num_cells_; }
int numFaces () const { return num_faces_; }
@@ -79,6 +70,9 @@ namespace Opm
size_t registerCellData( const std::string& name, const int components, const double initialValue = 0.0 );
size_t registerFaceData( const std::string& name, const int components, const double initialValue = 0.0 );
std::vector<double>& getCellData( const std::string& name );
const std::vector<double>& getCellData( const std::string& name ) const;
private:
int num_cells_;
int num_faces_;

4
opm/core/simulator/TwophaseState.hpp
View File

@@ -30,10 +30,6 @@ namespace Opm
class TwophaseState : public SimulatorState
{
public:
void setFirstSat(const std::vector<int>& cells,
const Opm::IncompPropertiesInterface& props,
ExtremalSat es);
virtual bool equals (const SimulatorState& other,
double epsilon = 1e-8) const;
};

7
opm/core/simulator/TwophaseState_impl.hpp
View File

@@ -4,13 +4,6 @@
namespace Opm {
inline void
TwophaseState::setFirstSat(const std::vector<int>& cells,
const Opm::IncompPropertiesInterface& props,
ExtremalSat es) {
SimulatorState::setFirstSat(cells, props, es);
}
inline bool
TwophaseState::equals (const SimulatorState& other,
double epsilon) const {

13
opm/core/simulator/initState.hpp
View File

@@ -22,6 +22,7 @@
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/core/simulator/SimulatorState.hpp>
struct UnstructuredGrid;
namespace Opm
@@ -35,6 +36,18 @@ namespace Opm
///
/// Functions for initializing a reservoir state.
/// Will initialize the first and second component of the
/// SATURATION field in all the cells in the set @cells. The
/// @props object will be queried, and depending on the value
/// @satType either the minimum or the maximum saturation is
/// applied to thee first component in the SATURATION field.
/// For the second component (1 - first_sat) is used.
enum ExtremalSat { MinSat, MaxSat };
template <class Props>
static void initSaturation(const std::vector<int>& cells , const Props& props , SimulatorState& state , ExtremalSat satType);
/// Initialize a two-phase state from parameters.
/// The following parameters are accepted (defaults):
/// - convection_testcase (false) -- Water in the 'left' part of the grid.

83
opm/core/simulator/initState_impl.hpp
View File

@@ -40,6 +40,38 @@
namespace Opm
{
template <class Props>
static void initSaturation(const std::vector<int>& cells , const Props& props , SimulatorState& state , ExtremalSat satType) {
const int num_phases = state.numPhases();
std::vector<double> min_sat(num_phases * cells.size());
std::vector<double> max_sat(num_phases * cells.size());
props.satRange(cells.size() ,cells.data() , min_sat.data() , max_sat.data());
{
std::vector<double> second_sat(cells.size());
std::vector<double> first_sat(cells.size());
for (size_t index = 0; index < cells.size(); index++) {
if (satType == MinSat) {
first_sat[index] = min_sat[ num_phases * index];
second_sat[index] = 1 - min_sat[ num_phases * index];
} else {
first_sat[index] = max_sat[ num_phases * index];
second_sat[index] = 1 - max_sat[ num_phases * index];
}
}
state.setCellDataComponent( "SATURATION" , 0 , cells , first_sat );
state.setCellDataComponent( "SATURATION" , 1 , cells , second_sat );
}
}
// Initialize saturations so that there is water below woc,
// and oil above.
namespace
{
#ifdef __clang__
@@ -75,8 +107,41 @@ namespace Opm
#pragma clang diagnostic pop
#endif /* __clang__ */
enum WaterInit { WaterBelow, WaterAbove };
/// Will initialize the first and second component of the
/// SATURATION field in all the cells in the set @cells. The
/// @props object will be queried, and depending on the value
/// @satType either the minimum or the maximum saturation is
/// applied to thee first component in the SATURATION field.
/// For the second component (1 - first_sat) is used.
template <class Props, class State>
static void initSaturation(const std::vector<int>& cells , const Props& props , State& state , ExtremalSat satType) {
std::vector<double> min_sat(state.numPhases() * cells.size());
std::vector<double> max_sat(state.numPhases() * cells.size());
props.satRange(cells.size() ,cells.data() , min_sat.data() , max_sat.data());
{
std::vector<double> first_sat(cells.size());
std::vector<double> second_sat(cells.size());
for (size_t index=0; index < cells.size(); index++) {
if (satType == MinSat) {
first_sat[index] = min_sat[index * state.numPhases()];
second_sat[index] = 1 - min_sat[index * state.numPhases()];
} else {
first_sat[index] = max_sat[index * state.numPhases()];
second_sat[index] = 1 - max_sat[index * state.numPhases()];
}
}
state.setCellDataComponent( "SATURATION" , 0 , cells , first_sat );
state.setCellDataComponent( "SATURATION" , 1 , cells , second_sat );
}
}
// Initialize saturations so that there is water below woc,
// and oil above.
// If invert is true, water is instead above, oil below.
@@ -106,9 +171,9 @@ namespace Opm
cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions,
woc, oil, water);
}
// Set saturations.
state.setFirstSat(oil, props, State::MinSat);
state.setFirstSat(water, props, State::MaxSat);
initSaturation( oil , props , state , MinSat );
initSaturation( water , props , state , MaxSat );
}
@@ -380,7 +445,10 @@ namespace Opm
for (int i = 0; i < num_cells; ++i) {
all_cells[i] = i;
}
state.setFirstSat(all_cells, props, State::MinSat);
initSaturation( all_cells , props , state , MinSat );
const bool convection_testcase = param.getDefault("convection_testcase", false);
const bool segregation_testcase = param.getDefault("segregation_testcase", false);
if (convection_testcase) {
@@ -394,7 +462,8 @@ namespace Opm
left_cells.push_back(cell);
}
}
state.setFirstSat(left_cells, props, State::MaxSat);
initSaturation( left_cells , props , state , MaxSat );
const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
std::fill(state.pressure().begin(), state.pressure().end(), init_p);
} else if (segregation_testcase) {
@@ -501,7 +570,7 @@ namespace Opm
for (int i = 0; i < num_cells; ++i) {
all_cells[i] = i;
}
state.setFirstSat(all_cells, props, State::MinSat);
initSaturation(all_cells , props , state , MinSat);
const bool convection_testcase = param.getDefault("convection_testcase", false);
if (convection_testcase) {
// Initialise water saturation to max in the 'left' part.
@@ -514,7 +583,7 @@ namespace Opm
left_cells.push_back(cell);
}
}
state.setFirstSat(left_cells, props, State::MaxSat);
initSaturation(left_cells , props , state , MaxSat );
const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
std::fill(state.pressure().begin(), state.pressure().end(), init_p);
} else if (param.has("water_oil_contact")) {

4
tutorials/tutorial3.cpp
View File

@@ -37,6 +37,7 @@
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
@@ -267,7 +268,8 @@ try
/// \internal [two-phase state]
TwophaseState state;
state.init(grid.number_of_cells , grid.number_of_faces, 2);
state.setFirstSat(allcells, props, TwophaseState::MinSat);
initSaturation( allcells , props , state , MinSat );
/// \internal [two-phase state]
/// \endinternal

5
tutorials/tutorial4.cpp
View File

@@ -37,6 +37,7 @@
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
@@ -214,7 +215,9 @@ try
/// \internal[two-phase state]
TwophaseState state;
state.init(grid.number_of_cells , grid.number_of_faces, 2);
state.setFirstSat(allcells, props, TwophaseState::MinSat);
initSaturation( allcells , props , state , MinSat );
/// \internal[two-phase state]
/// \endinternal