OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-16 20:24:48 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master' into PR/async-output

Browse Source
This commit is contained in:
Robert Kloefkorn 2016-03-30 10:45:40 +02:00
commit 9c67c0e135
31 changed files with 605 additions and 392 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@ -47,6 +47,9 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/VFPProdProperties.cpp
opm/autodiff/VFPInjProperties.cpp
opm/autodiff/WellMultiSegment.cpp
opm/autodiff/BlackoilSolventState.cpp
opm/polymer/PolymerState.cpp
opm/polymer/PolymerBlackoilState.cpp
opm/polymer/CompressibleTpfaPolymer.cpp
opm/polymer/IncompTpfaPolymer.cpp
opm/polymer/PolymerInflow.cpp

67
examples/sim_2p_incomp_ad.cpp
View File

@ -105,8 +105,9 @@ try
std::unique_ptr<GridManager> grid;
std::unique_ptr<IncompPropertiesInterface> props;
std::unique_ptr<RockCompressibility> rock_comp;
std::unique_ptr<TwophaseState> state;
EclipseStateConstPtr eclipseState;
TwophaseState state;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
@ -118,19 +119,25 @@ try
// Grid init
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
state.reset( new TwophaseState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid )));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(ug_grid, *props, param, gravity[2], *state);
} else {
initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
}
}
} else {
// Grid init.
@ -141,14 +148,20 @@ try
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
state.reset( new TwophaseState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid )));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(ug_grid, *props, param, gravity[2], *state);
}
}
// Warn if gravity but no density difference.
@ -170,7 +183,7 @@ try
// terms of total pore volume.
std::vector<double> porevol;
if (rock_comp->isActive()) {
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol);
} else {
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
}
@ -236,8 +249,8 @@ try
simtimer.init(param);
warnIfUnusedParams(param);
WellState well_state;
well_state.init(0, state);
rep = simulator.run(simtimer, state, well_state);
well_state.init(0, *state);
rep = simulator.run(simtimer, *state, well_state);
} else {
// With a deck, we may have more report steps etc.
WellState well_state;
@ -255,7 +268,7 @@ try
// properly transfer old well state to it every report step,
// since number of wells may change etc.
if (reportStepIdx == 0) {
well_state.init(wells.c_wells(), state);
well_state.init(wells.c_wells(), *state);
}
simtimer.setCurrentStepNum(reportStepIdx);
@ -273,7 +286,7 @@ try
if (reportStepIdx == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
if (output) {
epoch_rep.reportParam(epoch_os);
}

106
examples/sim_poly2p_comp_reorder.cpp
View File

@ -92,7 +92,7 @@ try
boost::scoped_ptr<RockCompressibility> rock_comp;
Opm::DeckConstPtr deck;
EclipseStateConstPtr eclipseState;
PolymerBlackoilState state;
std::unique_ptr<PolymerBlackoilState> state;
Opm::PolymerProperties poly_props;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
@ -106,23 +106,29 @@ try
// Grid init
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, ug_grid));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
state.reset( new PolymerBlackoilState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(ug_grid, *props, param, gravity[2], *state);
} else {
initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
}
initBlackoilSurfvol(ug_grid, *props, *state);
// Init polymer properties.
poly_props.readFromDeck(deck, eclipseState);
}
initBlackoilSurfvol(*grid->c_grid(), *props, state);
// Init polymer properties.
poly_props.readFromDeck(deck, eclipseState);
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
@ -132,31 +138,43 @@ try
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
initBlackoilSurfvol(*grid->c_grid(), *props, state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
double smin[2], smax[2];
props->satRange(1, &cell, smin, smax);
if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
state.concentration()[cell] = poly_init;
state.maxconcentration()[cell] = poly_init;
} else {
state.saturation()[2*cell + 0] = 0.;
state.saturation()[2*cell + 1] = 1.;
state.concentration()[cell] = 0.;
state.maxconcentration()[cell] = 0.;
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init.
props.reset(new BlackoilPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
state.reset( new PolymerBlackoilState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(ug_grid, *props, param, gravity[2], *state);
initBlackoilSurfvol(ug_grid, *props, *state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) {
double smin[2], smax[2];
auto& saturation = state->saturation();
auto& concentration = state->getCellData( state->CONCENTRATION );
auto& max_concentration = state->getCellData( state->CMAX );
props->satRange(1, &cell, smin, smax);
if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
concentration[cell] = poly_init;
max_concentration[cell] = poly_init;
} else {
saturation[2*cell + 0] = 0.;
saturation[2*cell + 1] = 1.;
concentration[cell] = 0.;
max_concentration[cell] = 0.;
}
}
}
}
// Init polymer properties.
// Setting defaults to provide a simple example case.
@ -240,8 +258,8 @@ try
simtimer.init(param);
warnIfUnusedParams(param);
WellState well_state;
well_state.init(0, state);
rep = simulator.run(simtimer, state, well_state);
well_state.init(0, *state);
rep = simulator.run(simtimer, *state, well_state);
} else {
// With a deck, we may have more epochs etc.
WellState well_state;
@ -284,7 +302,7 @@ try
// properly transfer old well state to it every report step,
// since number of wells may change etc.
if (reportStepIdx == 0) {
well_state.init(wells.c_wells(), state);
well_state.init(wells.c_wells(), *state);
}
// Create and run simulator.
@ -300,7 +318,7 @@ try
if (reportStepIdx == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
// Update total timing report and remember step number.
rep += epoch_rep;

101
examples/sim_poly2p_incomp_reorder.cpp
View File

@ -92,7 +92,7 @@ try
boost::scoped_ptr<IncompPropertiesInterface> props;
boost::scoped_ptr<RockCompressibility> rock_comp;
EclipseStateConstPtr eclipseState;
PolymerState state;
std::unique_ptr<PolymerState> state;
Opm::PolymerProperties poly_props;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
@ -106,22 +106,27 @@ try
// Grid init
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid ));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
state.reset( new PolymerState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(ug_grid, *props, param, gravity[2], *state);
} else {
initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
}
// Init polymer properties.
poly_props.readFromDeck(deck, eclipseState);
}
// Init polymer properties.
poly_props.readFromDeck(deck, eclipseState);
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
@ -131,28 +136,38 @@ try
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
double smin[2], smax[2];
props->satRange(1, &cell, smin, smax);
if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
state.concentration()[cell] = poly_init;
state.maxconcentration()[cell] = poly_init;
} else {
state.saturation()[2*cell + 0] = 0.;
state.saturation()[2*cell + 1] = 1.;
state.concentration()[cell] = 0.;
state.maxconcentration()[cell] = 0.;
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));;
state.reset( new PolymerState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(ug_grid, *props, param, gravity[2], *state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) {
double smin[2], smax[2];
auto& saturation = state->saturation();
auto& concentration = state->getCellData( state->CONCENTRATION );
auto& max_concentration = state->getCellData( state->CMAX );
props->satRange(1, &cell, smin, smax);
if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
concentration[cell] = poly_init;
max_concentration[cell] = poly_init;
} else {
saturation[2*cell + 0] = 0.;
saturation[2*cell + 1] = 1.;
concentration[cell] = 0.;
max_concentration[cell] = 0.;
}
}
}
}
@ -212,7 +227,7 @@ try
// terms of total pore volume.
std::vector<double> porevol;
if (rock_comp->isActive()) {
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol);
} else {
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
}
@ -276,8 +291,8 @@ try
simtimer.init(param);
warnIfUnusedParams(param);
WellState well_state;
well_state.init(0, state);
rep = simulator.run(simtimer, state, well_state);
well_state.init(0, *state);
rep = simulator.run(simtimer, *state, well_state);
} else {
// With a deck, we may have more epochs etc.
@ -321,7 +336,7 @@ try
// properly transfer old well state to it every report step,
// since number of wells may change etc.
if (reportStepIdx == 0) {
well_state.init(wells.c_wells(), state);
well_state.init(wells.c_wells(), *state);
}
// Create and run simulator.
@ -339,7 +354,7 @@ try
if (reportStepIdx == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
// Update total timing report and remember step number.
rep += epoch_rep;

12
examples/sim_poly_fi2p_comp_ad.cpp
View File

@ -115,7 +115,7 @@ try
std::shared_ptr<BlackoilPropertiesInterface> props;
std::shared_ptr<BlackoilPropsAdInterface> new_props;
std::shared_ptr<RockCompressibility> rock_comp;
PolymerBlackoilState state;
std::unique_ptr<PolymerBlackoilState> state;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
@ -187,6 +187,8 @@ try
Opm::UgGridHelpers::globalCell(cGrid),
Opm::UgGridHelpers::cartDims(cGrid),
param));
state.reset( new PolymerBlackoilState( Opm::UgGridHelpers::numCells(cGrid), Opm::UgGridHelpers::numFaces(cGrid), 2));
new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, materialLawManager, cGrid));
PolymerProperties polymer_props(deck, eclipseState);
PolymerPropsAd polymer_props_ad(polymer_props);
@ -199,10 +201,10 @@ try
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
initBlackoilSurfvol(*grid->c_grid(), *props, state);
initStateBasic(*grid->c_grid(), *props, param, gravity[2], *state);
initBlackoilSurfvol(*grid->c_grid(), *props, *state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], *state);
}
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
@ -254,7 +256,7 @@ try
deck,
*fis_solver,
grav);
fullReport= simulator.run(simtimer, state);
fullReport= simulator.run(simtimer, *state);
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.report(std::cout);

17
opm/autodiff/BackupRestore.hpp
View File

@ -21,8 +21,9 @@
#define OPM_BACKUPRESTORE_HEADER_INCLUDED
#include <iostream>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@ -159,7 +160,7 @@ namespace Opm {
BlackoilStateId = 2,
WellStateFullyImplicitBackoilId = 3 };
inline int objectId( const SimulatorState& /* state */) {
inline int objectId( const SimulationDataContainer& /* state */) {
return SimulatorStateId;
}
inline int objectId( const WellState& /* state */) {
@ -184,9 +185,9 @@ namespace Opm {
}
}
// SimulatorState
// SimulationDataContainer
inline
std::ostream& operator << (std::ostream& out, const SimulatorState& state )
std::ostream& operator << (std::ostream& out, const SimulationDataContainer& state )
{
// write id of object to stream
writeValue( out, objectId( state ) );
@ -205,7 +206,7 @@ namespace Opm {
}
inline
std::istream& operator >> (std::istream& in, SimulatorState& state )
std::istream& operator >> (std::istream& in, SimulationDataContainer& state )
{
// check id of stored object
checkObjectId( in, state );
@ -213,7 +214,7 @@ namespace Opm {
int numPhases = 0;
readValue( in, numPhases );
if( numPhases != state.numPhases() )
if( numPhases != (int) state.numPhases() )
OPM_THROW(std::logic_error,"num phases wrong");
// read variables
@ -234,7 +235,7 @@ namespace Opm {
writeValue( out, objectId( state ) );
// backup simulator state
const SimulatorState& simstate = static_cast< const SimulatorState& > (state);
const SimulationDataContainer& simstate = static_cast< const SimulationDataContainer& > (state);
out << simstate;
// backup additional blackoil state variables
@ -252,7 +253,7 @@ namespace Opm {
checkObjectId( in, state );
// restore simulator state
SimulatorState& simstate = static_cast< SimulatorState& > (state);
SimulationDataContainer& simstate = static_cast< SimulationDataContainer& > (state);
in >> simstate;
// restore additional blackoil state variables

4
opm/autodiff/BlackoilModelBase.hpp
View File

@ -45,7 +45,7 @@ namespace Opm {
class RockCompressibility;
class NewtonIterationBlackoilInterface;
class VFPProperties;
class SimulationDataContainer;
/// Struct for containing iteration variables.
struct DefaultBlackoilSolutionState
@ -207,7 +207,7 @@ namespace Opm {
/// \brief compute the relative change between to simulation states
// \return || u^n+1 - u^n || / || u^n+1 ||
double relativeChange( const SimulatorState& previous, const SimulatorState& current ) const;
double relativeChange( const SimulationDataContainer& previous, const SimulationDataContainer& current ) const;
/// The size (number of unknowns) of the nonlinear system of equations.
int sizeNonLinear() const;

5
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -48,6 +48,7 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <cassert>
#include <cmath>
#include <iostream>
@ -2520,8 +2521,8 @@ namespace detail {
template <class Grid, class Implementation>
double
BlackoilModelBase<Grid, Implementation>::
relativeChange(const SimulatorState& previous,
const SimulatorState& current ) const
relativeChange(const SimulationDataContainer& previous,
const SimulationDataContainer& current ) const
{
std::vector< double > p0 ( previous.pressure() );
std::vector< double > sat0( previous.saturation() );

14
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -137,11 +137,14 @@ namespace Opm {
// Initial polymer concentration.
if (has_solvent_) {
assert (not x.solvent_saturation().empty());
const int nc = x.solvent_saturation().size();
const V ss = Eigen::Map<const V>(&x.solvent_saturation()[0], nc);
const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL );
const int nc = solvent_saturation.size();
const V ss = Eigen::Map<const V>(solvent_saturation.data() , nc);
// This is some insanely detailed flickety flackety code;
// Solvent belongs after other reservoir vars but before well vars.
auto solvent_pos = vars0.begin() + fluid_.numPhases();
assert (not solvent_saturation.empty());
assert(solvent_pos == vars0.end() - 2);
vars0.insert(solvent_pos, ss);
}
@ -548,9 +551,10 @@ namespace Opm {
Base::updateState(modified_dx, reservoir_state, well_state);
// Update solvent.
const V ss_old = Eigen::Map<const V>(&reservoir_state.solvent_saturation()[0], nc, 1);
auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
const V ss = (ss_old - dss).max(zero);
std::copy(&ss[0], &ss[0] + nc, reservoir_state.solvent_saturation().begin());
std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
// adjust oil saturation
const Opm::PhaseUsage& pu = fluid_.phaseUsage();

34
opm/autodiff/BlackoilSolventState.cpp Normal file
View File

@ -0,0 +1,34 @@
/*
Copyright 2015 IRIS AS, 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/autodiff/BlackoilSolventState.hpp>
namespace Opm
{
const std::string BlackoilSolventState::SSOL = "SSOL";
BlackoilSolventState::BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases)
: BlackoilState( number_of_cells , number_of_faces , number_of_phases)
{
registerCellData( SSOL , 1 );
};
} // namespace Opm

23
opm/autodiff/BlackoilSolventState.hpp
View File

@ -20,9 +20,9 @@
#ifndef OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
#define OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/grid.h>
#include <vector>
#include <string>
#include <opm/core/simulator/BlackoilState.hpp>
namespace Opm
@ -33,22 +33,9 @@ namespace Opm
class BlackoilSolventState : 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 SSOL;
void init(int number_of_cells, int number_of_faces, int num_phases)
{
BlackoilState::init(number_of_cells, number_of_faces, num_phases);
solventId_ = SimulatorState::registerCellData( "SSOL", 1 );
}
std::vector<double>& solvent_saturation() { return cellData()[ solventId_ ]; }
const std::vector<double>& solvent_saturation() const { return cellData()[ solventId_ ]; }
private:
int solventId_;
BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases);
};
} // namespace Opm

49
opm/autodiff/FlowMain.hpp
View File

@ -183,7 +183,8 @@ namespace Opm
bool use_local_perm_ = true;
std::unique_ptr<DerivedGeology> geoprops_;
// setupState()
ReservoirState state_;
std::unique_ptr<ReservoirState> state_;
std::vector<double> threshold_pressures_;
// distributeData()
boost::any parallel_information_;
@ -441,8 +442,13 @@ namespace Opm
Opm::UgGridHelpers::cartDims(grid),
param_);
// Init state variables (saturation and pressure).
if (param_.has("init_saturation")) {
state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::numFaces(grid),
props.numPhases() ));
initStateBasic(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
@ -451,26 +457,35 @@ namespace Opm
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, param_, gravity_[2], state_);
props, param_, gravity_[2], *state_);
initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state_);
initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_);
enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
const int numPhases = props.numPhases();
const int numCells = Opm::UgGridHelpers::numCells(grid);
// Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState.
auto& gor = state_->getCellData( BlackoilState::GASOILRATIO );
const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL );
for (int c = 0; c < numCells; ++c) {
state_.gasoilratio()[c] = state_.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
/ state_.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
// Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field.
gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]];
}
}
} else if (deck_->hasKeyword("EQUIL") && props.numPhases() == 3) {
state_.init(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::numFaces(grid),
props.numPhases());
initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], state_);
state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
// Which state class are we really using - what a f... mess?
state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::numFaces(grid),
props.numPhases()));
initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], *state_);
//state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
} else {
state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::numFaces(grid),
props.numPhases()));
initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::numFaces(grid),
@ -478,12 +493,12 @@ namespace Opm
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, deck_, gravity_[2], state_);
props, deck_, gravity_[2], *state_);
}
// Threshold pressures.
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), state_, props, gravity_[2]);
computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), *state_, props, gravity_[2]);
threshold_pressures_ = thresholdPressures(deck_, eclipse_state_, grid, maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipse_state_, geoprops_->nnc(), maxDp);
threshold_pressures_.insert(threshold_pressures_.end(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
@ -493,9 +508,9 @@ namespace Opm
const int numCells = Opm::UgGridHelpers::numCells(grid);
std::vector<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> pc = state_.saturation();
props.capPress(numCells, state_.saturation().data(), cells.data(), pc.data(), nullptr);
fluidprops_->setSwatInitScaling(state_.saturation(), pc);
std::vector<double> pc = state_->saturation();
props.capPress(numCells, state_->saturation().data(), cells.data(), pc.data(), nullptr);
fluidprops_->setSwatInitScaling(state_->saturation(), pc);
}
}
@ -518,7 +533,7 @@ namespace Opm
// and initilialize new properties and states for it.
if (must_distribute_) {
distributeGridAndData(grid_init_->grid(), deck_, eclipse_state_,
state_, *fluidprops_, *geoprops_,
*state_, *fluidprops_, *geoprops_,
material_law_manager_, threshold_pressures_,
parallel_information_, use_local_perm_);
}
@ -617,7 +632,7 @@ namespace Opm
<< std::flush;
}
SimulatorReport fullReport = simulator_->run(simtimer, state_);
SimulatorReport fullReport = simulator_->run(simtimer, *state_);
if (output_cout_) {
std::cout << "\n\n================ End of simulation ===============\n\n";

95
opm/autodiff/ParallelDebugOutput.hpp
View File

@ -19,8 +19,10 @@
#ifndef OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
#define OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/core/grid.h>
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/wells/WellsManager.hpp>
@ -35,6 +37,7 @@ namespace Opm
{
class ParallelDebugOutputInterface
{
protected:
ParallelDebugOutputInterface () {}
@ -42,11 +45,11 @@ namespace Opm
virtual ~ParallelDebugOutputInterface() {}
// gather solution to rank 0 for EclipseWriter
virtual bool collectToIORank( const SimulatorState& localReservoirState,
virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
const WellState& localWellState,
const int reportStep ) = 0;
virtual const SimulatorState& globalReservoirState() const = 0 ;
virtual const SimulationDataContainer& globalReservoirState() const = 0 ;
virtual const WellState& globalWellState() const = 0 ;
virtual bool isIORank() const = 0;
virtual bool isParallel() const = 0;
@ -60,7 +63,7 @@ namespace Opm
protected:
const GridImpl& grid_;
const SimulatorState* globalState_;
const SimulationDataContainer* globalState_;
const WellState* wellState_;
public:
@ -71,7 +74,7 @@ namespace Opm
: grid_( grid ) {}
// gather solution to rank 0 for EclipseWriter
virtual bool collectToIORank( const SimulatorState& localReservoirState,
virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
const WellState& localWellState,
const int /* reportStep */)
{
@ -80,7 +83,7 @@ namespace Opm
return true ;
}
virtual const SimulatorState& globalReservoirState() const { return *globalState_; }
virtual const SimulationDataContainer& globalReservoirState() const { return *globalState_; }
virtual const WellState& globalWellState() const { return *wellState_; }
virtual bool isIORank () const { return true; }
virtual bool isParallel () const { return false; }
@ -243,7 +246,7 @@ namespace Opm
Dune::CpGrid& globalGrid = *grid_;
// initialize global state with correct sizes
globalReservoirState_.init( globalGrid.numCells(), globalGrid.numFaces(), numPhases );
globalReservoirState_.reset( new SimulationDataContainer( globalGrid.numCells(), globalGrid.numFaces(), numPhases ));
// copy global cartesian index
globalIndex_ = globalGrid.globalCell();
@ -306,18 +309,18 @@ namespace Opm
}
}
class PackUnPackSimulatorState : public P2PCommunicatorType::DataHandleInterface
class PackUnPackSimulationDataContainer : public P2PCommunicatorType::DataHandleInterface
{
const SimulatorState& localState_;
SimulatorState& globalState_;
const SimulationDataContainer& localState_;
SimulationDataContainer& globalState_;
const WellState& localWellState_;
WellState& globalWellState_;
const IndexMapType& localIndexMap_;
const IndexMapStorageType& indexMaps_;
public:
PackUnPackSimulatorState( const SimulatorState& localState,
SimulatorState& globalState,
PackUnPackSimulationDataContainer( const SimulationDataContainer& localState,
SimulationDataContainer& globalState,
const WellState& localWellState,
WellState& globalWellState,
const IndexMapType& localIndexMap,
@ -333,12 +336,11 @@ namespace Opm
if( isIORank )
{
// add missing data to global state
for( size_t i=globalState_.cellData().size();
i<localState.cellData().size(); ++i )
{
const size_t components = localState.cellData()[ i ].size() / localState.numCells();
assert( components * localState.numCells() == localState.cellData()[ i ].size() );
globalState_.registerCellData( localState.cellDataNames()[ i ], components );
for (const auto& pair : localState.cellData()) {
const std::string& key = pair.first;
if (!globalState_.hasCellData( key )) {
globalState_.registerCellData( key , localState.numCellDataComponents( key ));
}
}
MessageBufferType buffer;
@ -357,13 +359,10 @@ namespace Opm
}
// write all cell data registered in local state
const size_t numCells = localState_.numCells();
const size_t numCellData = localState_.cellData().size();
for( size_t d=0; d<numCellData; ++d )
{
const std::vector< double >& data = localState_.cellData()[ d ];
const size_t stride = data.size() / numCells ;
assert( numCells * stride == data.size() );
for (const auto& pair : localState_.cellData()) {
const std::string& key = pair.first;
const auto& data = pair.second;
const size_t stride = localState_.numCellDataComponents( key );
for( size_t i=0; i<stride; ++i )
{
@ -378,22 +377,20 @@ namespace Opm
void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
{
// read all cell data registered in local state
const size_t numCells = globalState_.numCells();
const size_t numCellData = globalState_.cellData().size();
for( size_t d=0; d<numCellData; ++d )
{
std::vector< double >& data = globalState_.cellData()[ d ];
const size_t stride = data.size() / numCells ;
assert( numCells * stride == data.size() );
// write all cell data registered in local state
for (auto& pair : globalState_.cellData()) {
const std::string& key = pair.first;
auto& data = pair.second;
const size_t stride = globalState_.numCellDataComponents( key );
for( size_t i=0; i<stride; ++i )
{
// write all data from local state to buffer
read( buffer, indexMap, data, i, stride );
//write all data from local state to buffer
read( buffer, indexMap, data, i, stride );
}
}
// read well data from buffer
readWells( buffer );
}
@ -516,7 +513,7 @@ namespace Opm
};
// gather solution to rank 0 for EclipseWriter
bool collectToIORank( const SimulatorState& localReservoirState,
bool collectToIORank( const SimulationDataContainer& localReservoirState,
const WellState& localWellState,
const int reportStep )
{
@ -536,10 +533,10 @@ namespace Opm
false);
const Wells* wells = wells_manager.c_wells();
globalWellState_.init(wells, globalReservoirState_, globalWellState_ );
globalWellState_.init(wells, *globalReservoirState_, globalWellState_ );
}
PackUnPackSimulatorState packUnpack( localReservoirState, globalReservoirState_,
PackUnPackSimulationDataContainer packUnpack( localReservoirState, *globalReservoirState_,
localWellState, globalWellState_,
localIndexMap_, indexMaps_, isIORank() );
@ -552,7 +549,7 @@ namespace Opm
return isIORank();
}
const SimulatorState& globalReservoirState() const { return globalReservoirState_; }
const SimulationDataContainer& globalReservoirState() const { return *globalReservoirState_; }
const WellState& globalWellState() const { return globalWellState_; }
bool isIORank() const
@ -573,18 +570,18 @@ namespace Opm
}
protected:
std::unique_ptr< Dune::CpGrid > grid_;
Opm::EclipseStateConstPtr eclipseState_;
const double* permeability_;
P2PCommunicatorType toIORankComm_;
IndexMapType globalIndex_;
IndexMapType localIndexMap_;
IndexMapStorageType indexMaps_;
SimulatorState globalReservoirState_;
std::unique_ptr< Dune::CpGrid > grid_;
Opm::EclipseStateConstPtr eclipseState_;
const double* permeability_;
P2PCommunicatorType toIORankComm_;
IndexMapType globalIndex_;
IndexMapType localIndexMap_;
IndexMapStorageType indexMaps_;
std::unique_ptr<SimulationDataContainer> globalReservoirState_;
// this needs to be revised
WellStateFullyImplicitBlackoil globalWellState_;
WellStateFullyImplicitBlackoil globalWellState_;
// true if we are on I/O rank
const bool isIORank_;
const bool isIORank_;
};
#endif // #if HAVE_DUNE_CORNERPOINT

16
opm/autodiff/RedistributeDataHandles.hpp
View File

@ -232,7 +232,7 @@ public:
{
assert( T::codimension == 0);
for ( int i=0; i<sendState_.numPhases(); ++i )
for ( size_t i=0; i<sendState_.numPhases(); ++i )
{
buffer.write(sendState_.surfacevol()[e.index()*sendState_.numPhases()+i]);
}
@ -240,7 +240,7 @@ public:
buffer.write(sendState_.rv()[e.index()]);
buffer.write(sendState_.pressure()[e.index()]);
buffer.write(sendState_.temperature()[e.index()]);
for ( int i=0; i<sendState_.numPhases(); ++i )
for ( size_t i=0; i<sendState_.numPhases(); ++i )
{
buffer.write(sendState_.saturation()[e.index()*sendState_.numPhases()+i]);
}
@ -257,11 +257,11 @@ public:
void scatter(B& buffer, const T& e, std::size_t size_arg)
{
assert( T::codimension == 0);
assert( int(size_arg) == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
assert( size_arg == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
static_cast<void>(size_arg);
double val;
for ( int i=0; i<recvState_.numPhases(); ++i )
for ( size_t i=0; i<recvState_.numPhases(); ++i )
{
buffer.read(val);
recvState_.surfacevol()[e.index()*sendState_.numPhases()+i]=val;
@ -274,7 +274,7 @@ public:
recvState_.pressure()[e.index()]=val;
buffer.read(val);
recvState_.temperature()[e.index()]=val;
for ( int i=0; i<recvState_.numPhases(); ++i )
for ( size_t i=0; i<recvState_.numPhases(); ++i )
{
buffer.read(val);
recvState_.saturation()[e.index()*sendState_.numPhases()+i]=val;
@ -411,8 +411,6 @@ void distributeGridAndData( Dune::CpGrid& grid,
boost::any& parallelInformation,
const bool useLocalPerm)
{
BlackoilState distributed_state;
Dune::CpGrid global_grid ( grid );
global_grid.switchToGlobalView();
@ -441,8 +439,8 @@ void distributeGridAndData( Dune::CpGrid& grid,
BlackoilPropsAdFromDeck distributed_props(properties,
distributed_material_law_manager,
grid.numCells());
distributed_state.init(grid.numCells(), grid.numFaces(), state.numPhases());
// init does not resize surfacevol. Do it manually.
BlackoilState distributed_state(grid.numCells(), grid.numFaces(), state.numPhases());
// construction does not resize surfacevol. Do it manually.
distributed_state.surfacevol().resize(grid.numCells()*state.numPhases(),
std::numeric_limits<double>::max());
BlackoilStateDataHandle state_handle(global_grid, grid,

53
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
View File

@ -23,6 +23,11 @@
#include "SimulatorFullyImplicitBlackoilOutput.hpp"
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/DataMap.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/common/ErrorMacros.hpp>
@ -32,9 +37,6 @@
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BackupRestore.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
#include <sstream>
#include <iomanip>
#include <fstream>
@ -52,7 +54,7 @@ namespace Opm
void outputStateVtk(const UnstructuredGrid& grid,
const SimulatorState& state,
const SimulationDataContainer& state,
const int step,
const std::string& output_dir)
{
@ -89,16 +91,15 @@ namespace Opm
void outputStateMatlab(const UnstructuredGrid& grid,
const Opm::SimulatorState& state,
const Opm::SimulationDataContainer& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
for( unsigned int i=0; i<state.cellDataNames().size(); ++ i )
{
const std::string& name = state.cellDataNames()[ i ];
for (const auto& pair : state.cellData()) {
const std::string& name = pair.first;
std::string key;
if( name == "SURFACEVOL" ) {
key = "surfvolume";
@ -113,7 +114,7 @@ namespace Opm
continue;
}
// set data to datmap
dm[ key ] = &state.cellData()[ i ];
dm[ key ] = &pair.second;
}
std::vector<double> cell_velocity;
@ -206,7 +207,7 @@ namespace Opm
#ifdef HAVE_DUNE_CORNERPOINT
void outputStateVtk(const Dune::CpGrid& grid,
const Opm::SimulatorState& state,
const Opm::SimulationDataContainer& state,
const int step,
const std::string& output_dir)
{
@ -296,7 +297,7 @@ namespace Opm
void
BlackoilOutputWriter::
writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& localState,
const SimulationDataContainer& localState,
const WellState& localWellState,
bool substep)
{
@ -312,7 +313,7 @@ namespace Opm
isIORank = parallelOutput_->collectToIORank( localState, localWellState, timer.reportStepNum() );
}
const SimulatorState& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
// serial output is only done on I/O rank
@ -379,7 +380,35 @@ namespace Opm
}
catch ( const std::bad_cast& e )
{
// store report step
lastBackupReportStep_ = reportStep;
// write resport step number
backupfile_.write( (const char *) &reportStep, sizeof(int) );
/*
try {
const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state);
backupfile_ << boState;
const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState);
backupfile_ << boWellState;
}
catch ( const std::bad_cast& e )
{
}
*/
/*
const WellStateFullyImplicitBlackoil* boWellState =
dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState);
if( boWellState ) {
backupfile_ << (*boWellState);
}
else
OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed");
*/
backupfile_ << std::flush;
}
/*

28
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@ -20,7 +20,6 @@
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/DataMap.hpp>
#include <opm/common/ErrorMacros.hpp>
@ -54,14 +53,17 @@
namespace Opm
{
class SimulationDataContainer;
class BlackoilState;
void outputStateVtk(const UnstructuredGrid& grid,
const Opm::SimulatorState& state,
const Opm::SimulationDataContainer& state,
const int step,
const std::string& output_dir);
void outputStateMatlab(const UnstructuredGrid& grid,
const Opm::SimulatorState& state,
const Opm::SimulationDataContainer& state,
const int step,
const std::string& output_dir);
@ -70,23 +72,23 @@ namespace Opm
const std::string& output_dir);
#ifdef HAVE_DUNE_CORNERPOINT
void outputStateVtk(const Dune::CpGrid& grid,
const Opm::SimulatorState& state,
const Opm::SimulationDataContainer& state,
const int step,
const std::string& output_dir);
#endif
template<class Grid>
void outputStateMatlab(const Grid& grid,
const Opm::SimulatorState& state,
const Opm::SimulationDataContainer& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
for( unsigned int i=0; i<state.cellDataNames().size(); ++ i )
for (const auto& pair : state.cellData())
{
const std::string& name = state.cellDataNames()[ i ];
const std::string& name = pair.first;
std::string key;
if( name == "SURFACEVOL" ) {
key = "surfvolume";
@ -101,7 +103,7 @@ namespace Opm
continue;
}
// set data to datmap
dm[ key ] = &state.cellData()[ i ];
dm[ key ] = &pair.second;
}
std::vector<double> cell_velocity;
@ -155,7 +157,7 @@ namespace Opm
/** \copydoc Opm::OutputWriter::writeTimeStep */
void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& state,
const SimulationDataContainer& state,
const WellState&,
bool /*substep*/ = false)
{
@ -185,7 +187,7 @@ namespace Opm
/** \copydoc Opm::OutputWriter::writeTimeStep */
void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& reservoirState,
const SimulationDataContainer& reservoirState,
const WellState& wellState,
bool /*substep*/ = false)
{
@ -215,7 +217,7 @@ namespace Opm
/** \copydoc Opm::OutputWriter::writeTimeStep */
void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& reservoirState,
const SimulationDataContainer& reservoirState,
const Opm::WellState& wellState,
bool substep = false);
@ -242,7 +244,7 @@ namespace Opm
void initFromRestartFile(const PhaseUsage& phaseusage,
const double* permeability,
const Grid& grid,
SimulatorState& simulatorstate,
SimulationDataContainer& simulatorstate,
WellStateFullyImplicitBlackoil& wellstate);
bool isRestart() const;
@ -325,7 +327,7 @@ namespace Opm
initFromRestartFile( const PhaseUsage& phaseusage,
const double* permeability,
const Grid& grid,
SimulatorState& simulatorstate,
SimulationDataContainer& simulatorstate,
WellStateFullyImplicitBlackoil& wellstate)
{
WellsManager wellsmanager(eclipseState_,

1
opm/autodiff/TransportSolverTwophaseAd.cpp
View File

@ -22,6 +22,7 @@
#include <opm/autodiff/TransportSolverTwophaseAd.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/common/ErrorMacros.hpp>

4
opm/polymer/CompressibleTpfaPolymer.cpp
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);
}

8
opm/polymer/IncompTpfaPolymer.cpp
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:

3
opm/polymer/IncompTpfaPolymer.hpp
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. ------

42
opm/polymer/PolymerBlackoilState.cpp Normal file
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

23
opm/polymer/PolymerBlackoilState.hpp
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

42
opm/polymer/PolymerState.cpp Normal file
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

31
opm/polymer/PolymerState.hpp
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

18
opm/polymer/SimulatorCompressiblePolymer.cpp
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);

26
opm/polymer/SimulatorPolymer.cpp
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;

36
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
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);

46
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
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,

10
opm/polymer/polymerUtilities.cpp
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);

4
tests/test_rateconverter.cpp
View File

@ -32,6 +32,7 @@
#include <boost/test/unit_test.hpp>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/utility/Units.hpp>
@ -107,8 +108,7 @@ BOOST_FIXTURE_TEST_CASE(ThreePhase, TestFixture<SetupSimple>)
Region reg{ 0 };
RCvrt cvrt(ad_props, reg);
Opm::BlackoilState x;
x.init(*grid.c_grid(), 3);
Opm::BlackoilState x( Opm::UgGridHelpers::numCells( *grid.c_grid()) , Opm::UgGridHelpers::numFaces( *grid.c_grid()) , 3);
cvrt.defineState(x);

76
tests/test_singlecellsolves.cpp
View File

@ -71,7 +71,7 @@ try
boost::scoped_ptr<GridManager> grid;
boost::scoped_ptr<IncompPropertiesInterface> props;
PolymerState state;
std::unique_ptr<PolymerState> state;
Opm::PolymerProperties poly_props;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
@ -81,23 +81,35 @@ try
// Grid init.
grid.reset(new GridManager(2, 1, 1, 1.0, 1.0, 1.0));
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, 0.0, state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
double smin[2], smax[2];
props->satRange(1, &cell, smin, smax);
if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
state.concentration()[cell] = poly_init;
state.maxconcentration()[cell] = poly_init;
} else {
state.saturation()[2*cell + 0] = 0.;
state.saturation()[2*cell + 1] = 1.;
state.concentration()[cell] = 0.;
state.maxconcentration()[cell] = 0.;
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
state.reset( new PolymerState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, 0.0, *state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
double smin[2], smax[2];
props->satRange(1, &cell, smin, smax);
auto& saturation = state->saturation();
auto& concentration = state->getCellData( state->CONCENTRATION );
auto& max_concentration = state->getCellData( state->CMAX );
props->satRange(1, &cell, smin, smax);
if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
concentration[cell] = poly_init;
max_concentration[cell] = poly_init;
} else {
saturation[2*cell + 0] = 0.;
saturation[2*cell + 1] = 1.;
concentration[cell] = 0.;
max_concentration[cell] = 0.;
}
}
}
}
@ -210,7 +222,7 @@ try
if (face01 == -1) {
OPM_THROW(std::runtime_error, "Could not find face adjacent to cells [0 1]");
}
state.faceflux()[face01] = src[0];
state->faceflux()[face01] = src[0];
for (int sats = 0; sats < num_sats; ++sats) {
const double s = double(sats)/double(num_sats - 1);
const double ff = s; // Simplified a lot...
@ -220,20 +232,26 @@ try
// std::cout << "(s, c) = (" << s << ", " << c << ")\n";
transport_src[0] = src[0]*ff;
// Resetting the state for next run.
state.saturation()[0] = 0.0;
state.saturation()[1] = 0.0;
state.concentration()[0] = 0.0;
state.concentration()[1] = 0.0;
state.maxconcentration()[0] = 0.0;
state.maxconcentration()[1] = 0.0;
reorder_model.solve(&state.faceflux()[0],
auto& saturation = state->saturation();
auto& concentration = state->getCellData( state->CONCENTRATION );
auto& max_concentration = state->getCellData( state->CMAX );
saturation[0] = 0.0;
saturation[1] = 0.0;
concentration[0] = 0.0;
concentration[1] = 0.0;
max_concentration[0] = 0.0;
max_concentration[1] = 0.0;
reorder_model.solve(&state->faceflux()[0],
&porevol[0],
&transport_src[0],
&polymer_inflow_c[0],
dt,
state.saturation(),
state.concentration(),
state.maxconcentration());
saturation,
concentration,
max_concentration);
#ifdef PROFILING
// Extract residual counts.