Replaced SimulatorState -> SimulationDataContainer

This commit is contained in:
Joakim Hove
2016-02-28 11:27:00 +01:00
parent 0e0e92aa71
commit 18c07d5d66
30 changed files with 569 additions and 383 deletions

View File

@@ -88,8 +88,8 @@ namespace Opm
PolymerBlackoilState& state,
WellState& well_state)
{
c_ = &state.concentration();
cmax_ = &state.maxconcentration();
c_ = &state.getCellData( state.CONCENTRATION );
cmax_ = &state.getCellData( state.CMAX );
CompressibleTpfa::solve(dt, state, well_state);
}

View File

@@ -20,6 +20,8 @@
#include <config.h>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
@@ -99,8 +101,8 @@ namespace Opm
PolymerState& state,
WellState& well_state)
{
c_ = &state.concentration();
cmax_ = &state.maxconcentration();
c_ = &state.getCellData( state.CONCENTRATION );
cmax_ = &state.getCellData( state.CMAX) ;
if (rock_comp_props_ != 0 && rock_comp_props_->isActive()) {
solveRockComp(dt, state, well_state);
} else {
@@ -116,7 +118,7 @@ namespace Opm
/// Compute per-solve dynamic properties.
void IncompTpfaPolymer::computePerSolveDynamicData(const double /*dt*/,
const SimulatorState& state,
const SimulationDataContainer& state,
const WellState& /*well_state*/)
{
// Computed here:

View File

@@ -37,6 +37,7 @@ namespace Opm
class LinearSolverInterface;
class PolymerState;
class WellState;
class SimulationDataContainer;
/// Encapsulating a tpfa pressure solver for the incompressible-fluid case with polymer.
/// Supports gravity, wells controlled by bhp or reservoir rates,
@@ -96,7 +97,7 @@ namespace Opm
private:
virtual void computePerSolveDynamicData(const double dt,
const SimulatorState& state,
const SimulationDataContainer& state,
const WellState& well_state);
private:
// ------ Data that will remain unmodified after construction. ------

View File

@@ -0,0 +1,42 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
namespace Opm
{
const std::string PolymerBlackoilState::CONCENTRATION = "CONCENTRATION";
const std::string PolymerBlackoilState::CMAX = "CMAX";
PolymerBlackoilState::PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases) :
BlackoilState( number_of_cells , number_of_faces, num_phases)
{
registerCellData(CONCENTRATION , 1 , 0 );
registerCellData(CMAX , 1 , 0 );
}
} // namespace Opm

View File

@@ -33,27 +33,10 @@ namespace Opm
class PolymerBlackoilState : public BlackoilState
{
public:
void init(const UnstructuredGrid& g, int num_phases)
{
this->init(g.number_of_cells, g.number_of_faces, num_phases);
}
static const std::string CONCENTRATION;
static const std::string CMAX;
void init(int number_of_cells, int number_of_faces, int num_phases)
{
BlackoilState::init(number_of_cells, number_of_faces, num_phases);
concentration_.resize(number_of_cells, 0.0);
cmax_.resize(number_of_cells, 0.0);
}
std::vector<double>& concentration() { return concentration_; }
std::vector<double>& maxconcentration() { return cmax_; }
const std::vector<double>& concentration() const { return concentration_; }
const std::vector<double>& maxconcentration() const { return cmax_; }
private:
std::vector<double> concentration_;
std::vector<double> cmax_;
PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases);
};
} // namespace Opm

View File

@@ -0,0 +1,42 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/polymer/PolymerState.hpp>
namespace Opm
{
const std::string PolymerState::CONCENTRATION = "CONCENTRATION";
const std::string PolymerState::CMAX = "CMAX";
PolymerState::PolymerState(int number_of_cells, int number_of_faces, int num_phases) :
SimulationDataContainer( number_of_cells , number_of_faces , num_phases )
{
registerCellData(CONCENTRATION , 1 , 0 );
registerCellData(CMAX , 1 , 0 );
}
} // namespace Opm

View File

@@ -20,8 +20,8 @@
#ifndef OPM_POLYMERSTATE_HEADER_INCLUDED
#define OPM_POLYMERSTATE_HEADER_INCLUDED
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/grid.h>
#include <vector>
@@ -29,34 +29,13 @@ namespace Opm
{
/// Simulator state for a two-phase simulator with polymer.
class PolymerState : public SimulatorState
class PolymerState : public SimulationDataContainer
{
public:
static const std::string CONCENTRATION;
static const std::string CMAX;
void init(int number_of_cells, int number_of_faces, int num_phases)
{
SimulatorState::init( number_of_cells , number_of_faces , num_phases );
registerCellData("CONCENTRATION" , 1 , 0 );
registerCellData("CMAX" , 1 , 0 );
}
const std::vector<double>& concentration() const {
return getCellData("CONCENTRATION");
}
const std::vector<double>& maxconcentration() const {
return getCellData("CMAX");
}
std::vector<double>& concentration() {
return getCellData("CONCENTRATION");
}
std::vector<double>& maxconcentration() {
return getCellData("CMAX");
}
PolymerState(int number_of_cells, int number_of_faces, int num_phases);
};
} // namespace Opm

View File

@@ -268,7 +268,7 @@ namespace Opm
total_timer.start();
double init_surfvol[2] = { 0.0 };
double inplace_surfvol[2] = { 0.0 };
double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
double polymass_adsorbed = computePolymerAdsorbed(grid_, props_, poly_props_, state, rock_comp_props_);
double init_polymass = polymass + polymass_adsorbed;
double tot_injected[2] = { 0.0 };
@@ -301,7 +301,7 @@ namespace Opm
if (check_well_controls_) {
computeFractionalFlow(props_, poly_props_, allcells_,
state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
state.concentration(), state.maxconcentration(),
state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ) ,
fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
@@ -397,7 +397,7 @@ namespace Opm
state.pressure(), state.temperature(), &initial_porevol[0], &porevol[0],
&transport_src[0], &polymer_inflow_c[0], stepsize,
state.saturation(), state.surfacevol(),
state.concentration(), state.maxconcentration());
state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
double substep_polyinj = 0.0;
@@ -416,7 +416,7 @@ namespace Opm
if (gravity_ != 0 && use_segregation_split_) {
tsolver_.solveGravity(columns_, stepsize,
state.saturation(), state.surfacevol(),
state.concentration(), state.maxconcentration());
state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
}
}
transport_timer.stop();
@@ -426,7 +426,7 @@ namespace Opm
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
state, rock_comp_props_);
tot_injected[0] += injected[0];
@@ -539,8 +539,8 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
dm["cmax"] = &state.getCellData( state.CMAX ) ;
dm["surfvol"] = &state.surfacevol();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
@@ -556,8 +556,8 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
dm["cmax"] = &state.getCellData( state.CMAX ) ;
dm["surfvol"] = &state.surfacevol();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);

View File

@@ -292,8 +292,8 @@ namespace Opm
total_timer.start();
double init_satvol[2] = { 0.0 };
double satvol[2] = { 0.0 };
double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX ));
double init_polymass = polymass + polymass_adsorbed;
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
@@ -330,7 +330,7 @@ namespace Opm
// Solve pressure.
if (check_well_controls_) {
computeFractionalFlow(props_, poly_props_, allcells_,
state.saturation(), state.concentration(), state.maxconcentration(),
state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ),
fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
@@ -425,7 +425,7 @@ namespace Opm
injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
state.saturation(), state.concentration(), state.maxconcentration());
state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
Opm::computeInjectedProduced(props_, poly_props_,
state,
transport_src, polymer_inflow_c, stepsize,
@@ -438,7 +438,7 @@ namespace Opm
polyprod += substep_polyprod;
if (use_segregation_split_) {
tsolver_.solveGravity(columns_, &porevol[0], stepsize,
state.saturation(), state.concentration(), state.maxconcentration());
state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
}
}
transport_timer.stop();
@@ -448,8 +448,8 @@ namespace Opm
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX ));
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
@@ -559,8 +559,8 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
dm["cmax"] = &state.getCellData( state.CMAX ) ;
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
@@ -575,8 +575,8 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
dm["cmax"] = &state.getCellData( state.CMAX ) ;
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
@@ -611,8 +611,8 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
dm["cmax"] = &state.getCellData( state.CMAX ) ;
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;

View File

@@ -125,8 +125,10 @@ namespace Opm {
WellState& well_state)
{
Base::prepareStep(dt, reservoir_state, well_state);
auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX );
// Initial max concentration of this time step from PolymerBlackoilState.
cmax_ = Eigen::Map<const V>(reservoir_state.maxconcentration().data(), Opm::AutoDiffGrid::numCells(grid_));
cmax_ = Eigen::Map<const V>(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_));
}
@@ -168,9 +170,10 @@ namespace Opm {
// Initial polymer concentration.
if (has_polymer_) {
assert (not x.concentration().empty());
const int nc = x.concentration().size();
const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
const auto& concentration = x.getCellData( x.CONCENTRATION );
assert (not concentration.empty());
const int nc = concentration.size();
const V c = Eigen::Map<const V>(concentration.data() , nc);
// Concentration belongs after other reservoir vars but before well vars.
auto concentration_pos = vars0.begin() + fluid_.numPhases();
assert(concentration_pos == vars0.end() - 2);
@@ -255,12 +258,14 @@ namespace Opm {
template <class Grid>
void BlackoilPolymerModel<Grid>::computeCmax(ReservoirState& state)
{
const int nc = AutoDiffGrid::numCells(grid_);
V tmp = V::Zero(nc);
for (int i = 0; i < nc; ++i) {
tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
}
std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
auto& max_concentration = state.getCellData( state.CMAX );
const auto& concentration = state.getCellData( state.CONCENTRATION );
std::transform( max_concentration.begin() ,
max_concentration.end() ,
concentration.begin() ,
max_concentration.begin() ,
[](double c_max , double c) { return std::max( c_max , c ); });
}
@@ -397,10 +402,13 @@ namespace Opm {
// Call base version.
Base::updateState(modified_dx, reservoir_state, well_state);
// Update concentration.
const V c_old = Eigen::Map<const V>(&reservoir_state.concentration()[0], nc, 1);
const V c = (c_old - dc).max(zero);
std::copy(&c[0], &c[0] + nc, reservoir_state.concentration().begin());
{
auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION );
// Update concentration.
const V c_old = Eigen::Map<const V>(concentration.data(), nc, 1);
const V c = (c_old - dc).max(zero);
std::copy(&c[0], &c[0] + nc, concentration.begin());
}
} else {
// Just forward call to base version.
Base::updateState(dx, reservoir_state, well_state);

View File

@@ -206,7 +206,7 @@ namespace {
const std::vector<double>& polymer_inflow = xw.polymerInflow();
// Initial max concentration of this time step from PolymerBlackoilState.
cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_));
cmax_ = Eigen::Map<V>(&x.getCellData( x.CMAX )[0], Opm::AutoDiffGrid::numCells(grid_));
const SolutionState state = constantState(x, xw);
computeAccum(state, 0);
@@ -366,9 +366,18 @@ namespace {
state.saturation[1] = ADB::constant(so, bpat);
// Concentration
assert(not x.concentration().empty());
const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
state.concentration = ADB::constant(c);
{
auto& concentration = x.getCellData( x.CONCENTRATION );
assert(concentration.empty());
const V c = Eigen::Map<const V>(concentration.data(), nc);
// Do not understand:
//concentration = ADB::constant(c);
// Old code based on concentraton() method had the statement:
//
// state.concentration = ADB::constant(c)
//
// This looks like it was a method assignment - how did it even compile?
}
// Well rates.
assert (not xw.wellRates().empty());
@@ -413,9 +422,12 @@ namespace {
vars0.push_back(sw0);
// Initial concentration.
assert (not x.concentration().empty());
const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
vars0.push_back(c);
{
auto& concentration = x.getCellData( x.CONCENTRATION );
assert (not concentration.empty());
const V c = Eigen::Map<const V>(concentration.data() , nc);
vars0.push_back(c);
}
// Initial well rates.
assert (not xw.wellRates().empty());
@@ -511,11 +523,14 @@ namespace {
{
const int nc = grid_.number_of_cells;
V tmp = V::Zero(nc);
const auto& concentration = state.getCellData( state.CONCENTRATION );
auto& cmax = state.getCellData( state.CMAX );
for (int i = 0; i < nc; ++i) {
tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
tmp[i] = std::max(cmax[i], concentration[i]);
}
std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
std::copy(&tmp[0], &tmp[0] + nc, cmax.begin());
}
@@ -764,11 +779,14 @@ namespace {
// Concentration updates.
// const double dcmax = 0.3 * polymer_props_ad_.cMax();
// std::cout << "\n the max concentration: " << dcmax / 0.3 << std::endl;
const V c_old = Eigen::Map<const V>(&state.concentration()[0], nc, 1);
// const V dc_limited = sign(dc) * dc.abs().min(dcmax);
// const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02());
const V c = (c_old - dc).max(zero);
std::copy(&c[0], &c[0] + nc, state.concentration().begin());
{
auto& concentration = state.getCellData( state.CONCENTRATION );
const V c_old = Eigen::Map<const V>(concentration.data() , nc, 1);
// const V dc_limited = sign(dc) * dc.abs().min(dcmax);
// const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02());
const V c = (c_old - dc).max(zero);
std::copy(&c[0], &c[0] + nc, concentration.begin());
}
// Qs update.
// Since we need to update the wellrates, that are ordered by wells,

View File

@@ -211,8 +211,8 @@ namespace Opm
OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells.");
}
const std::vector<double>& s = state.saturation();
const std::vector<double>& c = state.concentration();
const std::vector<double>& cmax = state.maxconcentration();
const std::vector<double>& c = state.getCellData( state.CONCENTRATION );
const std::vector<double>& cmax = state.getCellData( state.CMAX );
std::fill(injected, injected + np, 0.0);
std::fill(produced, produced + np, 0.0);
polyinj = 0.0;
@@ -282,8 +282,8 @@ namespace Opm
const std::vector<double>& temp = state.temperature();
const std::vector<double>& s = state.saturation();
const std::vector<double>& z = state.surfacevol();
const std::vector<double>& c = state.concentration();
const std::vector<double>& cmax = state.maxconcentration();
const std::vector<double>& c = state.getCellData( state.CONCENTRATION );
const std::vector<double>& cmax = state.getCellData( state.CMAX );
std::fill(injected, injected + np, 0.0);
std::fill(produced, produced + np, 0.0);
polyinj = 0.0;
@@ -406,7 +406,7 @@ namespace Opm
porosity.assign(props.porosity(), props.porosity() + num_cells);
}
double abs_mass = 0.0;
const std::vector<double>& cmax = state.maxconcentration();
const std::vector<double>& cmax = state.getCellData( state.CMAX );
for (int cell = 0; cell < num_cells; ++cell) {
double c_ads;
polyprops.simpleAdsorption(cmax[cell], c_ads);