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

DerivedGeology: pass the EclipseState object to its constructor

Browse Source 
this is required to implement pore volume and permeability multipliers
as discussed with [at]bska and [at]joakim-hove.

Note that this implies that the DerivedGeology class can't be
instantiated anymore if there is no EclipseState object. Thus all code
paths and tests that don't load a deck are removed by this patch. If
this is undesireable, there are two options: First, don't require
EclipseState for DerivedGeology which would imply to make the about 10
required multiplier functions part of the
BlackoilPropertiesAdInterface, or second, one can copy-and-paste the
DerivedGeology class as it was before this patch, derive from a newly
introduced DerivedGeologyInterface and pass DerivedGeologyInterface
objects to the simulator. IMHO, the second solution would be a bit
better but it would involve substantial overhead to implement and to
maintain it.

Anyway, in the mean time simulators cannot be instantiated without
decks.
This commit is contained in:
Andreas Lauser 2014-07-10 11:43:46 +02:00
parent 543d8d75b6
commit 63eaecf246
7 changed files with 60 additions and 373 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -75,8 +75,6 @@ list (APPEND EXAMPLE_SOURCE_FILES
examples/sim_2p_comp_ad.cpp examples/sim_2p_comp_ad.cpp
examples/sim_2p_incomp_ad.cpp examples/sim_2p_incomp_ad.cpp
examples/sim_simple.cpp examples/sim_simple.cpp
examples/test_impestpfa_ad.cpp
examples/test_implicit_ad.cpp
) )
# programs listed here will not only be compiled, but also marked for # programs listed here will not only be compiled, but also marked for

181
examples/sim_2p_comp_ad.cpp
View File

@ -86,13 +86,6 @@ try
parameter::ParameterGroup param(argc, argv, false); parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl; std::cout << "--------------- Reading parameters ---------------" << std::endl;
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
if (!use_deck) {
// This check should be removed when and if this simulator is verified and works without decks.
// The current code for the non-deck case fails for unknown reasons.
OPM_THROW(std::runtime_error, "This simulator cannot run without a deck with wells. Use deck_filename to specify deck.");
}
Opm::DeckConstPtr deck; Opm::DeckConstPtr deck;
std::unique_ptr<GridManager> grid; std::unique_ptr<GridManager> grid;
std::unique_ptr<BlackoilPropertiesInterface> props; std::unique_ptr<BlackoilPropertiesInterface> props;
@ -102,84 +95,33 @@ try
// bool check_well_controls = false; // bool check_well_controls = false;
// int max_well_control_iterations = 0; // int max_well_control_iterations = 0;
double gravity[3] = { 0.0 }; double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr parser(new Opm::Parser());
deck = parser->parseFile( deck_filename );
eclipseState.reset(new EclipseState(deck));
// Grid init std::string deck_filename = param.get<std::string>("deck_filename");
grid.reset(new GridManager(deck)); Opm::ParserPtr parser(new Opm::Parser());
// Rock and fluid init deck = parser->parseFile( deck_filename );
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param)); eclipseState.reset(new EclipseState(deck));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Grid init
// Rock compressibility. grid.reset(new GridManager(deck));
rock_comp.reset(new RockCompressibility(deck)); // Rock and fluid init
// Gravity. props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; // check_well_controls = param.getDefault("check_well_controls", false);
// Init state variables (saturation and pressure). // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
if (param.has("init_saturation")) { // Rock compressibility.
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); rock_comp.reset(new RockCompressibility(deck));
} else { // Gravity.
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
} // Init state variables (saturation and pressure).
initBlackoilSurfvol(*grid->c_grid(), *props, state); if (param.has("init_saturation")) {
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0);
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); 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);
} }
initBlackoilSurfvol(*grid->c_grid(), *props, state);
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0; const double *grav = use_gravity ? &gravity[0] : 0;
// Initialising wells if not from deck.
Wells* simple_wells = 0;
if (!use_deck) {
// Compute pore volumes, in order to enable specifying injection rate
// terms of total pore volume.
std::vector<double> porevol;
if (rock_comp->isActive()) {
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
} else {
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
}
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
*tot_porevol_init/unit::day;
simple_wells = create_wells(2, 2, 2);
const double inj_frac[2] = { 1.0, 0.0 };
const int inj_cell = 0;
const double WI = 1e-8; // This is a completely made-up number.
const double all_fluids[2] = { 1.0, 1.0 };
int ok = add_well(INJECTOR, 0.0, 1, inj_frac, &inj_cell, &WI, "Injector", simple_wells);
ok = ok && append_well_controls(SURFACE_RATE, 0.01*flow_per_sec, all_fluids, 0, simple_wells);
const int prod_cell = grid->c_grid()->number_of_cells - 1;
ok = ok && add_well(PRODUCER, 0.0, 1, NULL, &prod_cell, &WI, "Producer", simple_wells);
ok = ok && append_well_controls(SURFACE_RATE, -0.01*flow_per_sec, all_fluids, 1, simple_wells);
if (!ok) {
OPM_THROW(std::runtime_error, "Simple well init failed.");
}
well_controls_set_current( simple_wells->ctrls[0] , 0);
well_controls_set_current( simple_wells->ctrls[1] , 0);
}
// Linear solver. // Linear solver.
LinearSolverFactory linsolver(param); LinearSolverFactory linsolver(param);
@ -209,10 +151,30 @@ try
std::cout << "\n\n================ Starting main simulation loop ===============\n"; std::cout << "\n\n================ Starting main simulation loop ===============\n";
SimulatorReport rep; SimulatorReport rep;
Opm::DerivedGeology geology(*grid->c_grid(), *props); // With a deck, we may have more epochs etc.
if (!use_deck) { WellState well_state;
// Simple simulation without a deck. Opm::TimeMapConstPtr timeMap = eclipseState->getSchedule()->getTimeMap();
WellsManager wells(simple_wells); Opm::DerivedGeology geology(*grid->c_grid(), *props, eclipseState);
SimulatorTimer simtimer;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of report step.
std::cout << "\n\n-------------- Starting report step " << reportStepIdx << " --------------"
<< "\n (number of steps left: "
<< timeMap->numTimesteps() - reportStepIdx << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(eclipseState, reportStepIdx, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// 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);
}
simtimer.setCurrentStepNum(reportStepIdx);
// Create and run simulator.
SimulatorCompressibleAd simulator(param, SimulatorCompressibleAd simulator(param,
*grid->c_grid(), *grid->c_grid(),
geology, geology,
@ -221,54 +183,15 @@ try
wells, wells,
linsolver, linsolver,
grav); grav);
SimulatorTimer simtimer; if (reportStepIdx == 0) {
simtimer.init(param); warnIfUnusedParams(param);
warnIfUnusedParams(param);
WellState well_state;
well_state.init(simple_wells, state);
rep = simulator.run(simtimer, state, well_state);
} else {
// With a deck, we may have more epochs etc.
WellState well_state;
Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
SimulatorTimer simtimer;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of report step.
std::cout << "\n\n-------------- Starting report step " << reportStepIdx << " --------------"
<< "\n (number of steps left: "
<< timeMap->numTimesteps() - reportStepIdx << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(eclipseState, reportStepIdx, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// 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);
}
simtimer.setCurrentStepNum(reportStepIdx);
// Create and run simulator.
SimulatorCompressibleAd simulator(param,
*grid->c_grid(),
geology,
*props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
linsolver,
grav);
if (reportStepIdx == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
if (output) {
epoch_rep.reportParam(epoch_os);
}
// Update total timing report and remember step number.
rep += epoch_rep;
} }
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
if (output) {
epoch_rep.reportParam(epoch_os);
}
// Update total timing report and remember step number.
rep += epoch_rep;
} }
std::cout << "\n\n================ End of simulation ===============\n\n"; std::cout << "\n\n================ End of simulation ===============\n\n";

2
examples/sim_fibo_ad.cpp
View File

@ -193,7 +193,7 @@ try
// initialize variables // initialize variables
simtimer.init(timeMap); simtimer.init(timeMap);
Opm::DerivedGeology geology(*grid->c_grid(), *new_props); Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState);
SimulatorReport fullReport; SimulatorReport fullReport;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) { for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {

2
examples/sim_fibo_ad_cp.cpp
View File

@ -236,7 +236,7 @@ try
// initialize variables // initialize variables
simtimer.init(timeMap); simtimer.init(timeMap);
Opm::DerivedGeology geology(*grid, *new_props); Opm::DerivedGeology geology(*grid, *new_props, eclipseState);
SimulatorReport fullReport; SimulatorReport fullReport;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) { for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {

110
examples/test_impestpfa_ad.cpp
View File

@ -1,110 +0,0 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013 Statoil ASA.
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 <config.h>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/ImpesTPFAAD.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/wells.h>
#include <algorithm>
#include <iostream>
int
main(int argc, char* argv[])
try
{
const Opm::parameter::ParameterGroup param(argc, argv, false);
const Opm::GridManager gm(5, 5);
const UnstructuredGrid* g = gm.c_grid();
const int nc = g->number_of_cells;
const Opm::BlackoilPropertiesBasic oldprops(param, 2, nc);
const Opm::BlackoilPropsAd props(oldprops);
Wells* wells = create_wells(2, 2, 5);
const double inj_frac[] = { 1.0, 0.0 };
const double prod_frac[] = { 0.0, 0.0 };
const int num_inj = 3;
const int inj_cells[num_inj] = { 0, 1, 2 };
const int num_prod = 2;
const int prod_cells[num_prod] = { 20, 21 };
const double WI[3] = { 1e-12, 1e-12, 1e-12 };
bool ok = add_well(INJECTOR, 0.0, num_inj, inj_frac, inj_cells, WI, "Inj", wells);
ok = ok && add_well(PRODUCER, 0.0, num_prod, prod_frac, prod_cells, WI, "Prod", wells);
ok = ok && append_well_controls(BHP, 500.0*Opm::unit::barsa, 0, 0, wells);
// ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
double oildistr[2] = { 0.0, 1.0 };
ok = ok && append_well_controls(SURFACE_RATE, 1e-3, oildistr, 1, wells);
if (!ok) {
OPM_THROW(std::runtime_error, "Something went wrong with well init.");
}
set_current_control(0, 0, wells);
set_current_control(1, 0, wells);
double grav[] = { /*1.0*/ 0.0, 0.0 };
Opm::DerivedGeology geo(*g, props, grav);
Opm::LinearSolverFactory linsolver(param);
Opm::ImpesTPFAAD ps(*g, props, geo, *wells, linsolver);
Opm::BlackoilState state;
initStateBasic(*g, oldprops, param, 0.0, state);
initBlackoilSurfvol(*g, oldprops, state);
Opm::WellState well_state;
well_state.init(wells, state);
ps.solve(1.0, state, well_state);
std::cout << "Cell pressure:" << std::endl;
std::copy(state.pressure().begin(), state.pressure().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
std::cout << "Face flux:" << std::endl;
std::copy(state.faceflux().begin(), state.faceflux().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
std::cout << "Well bhp pressure:" << std::endl;
std::copy(well_state.bhp().begin(), well_state.bhp().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
return 0;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

127
examples/test_implicit_ad.cpp
View File

@ -1,127 +0,0 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013 Statoil ASA.
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 <config.h>
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <memory>
#include <algorithm>
#include <iostream>
namespace {
std::shared_ptr<Wells>
createWellConfig()
{
std::shared_ptr<Wells> wells(create_wells(2, 2, 2),
destroy_wells);
const double inj_frac[] = { 1.0, 0.0 };
const double prod_frac[] = { 0.0, 0.0 };
const int num_inj = 1;
const int inj_cells[num_inj] = { 0 };
const int num_prod = 1;
const int prod_cells[num_prod] = { 19 };
const double WI[3] = { 1e-12, 1e-12, 1e-12 };
bool ok = add_well(INJECTOR, 0.0, num_inj, inj_frac, inj_cells, WI, "Inj", wells.get());
ok = ok && add_well(PRODUCER, 0.0, num_prod, prod_frac, prod_cells, WI, "Prod", wells.get());
ok = ok && append_well_controls(BHP, 500.0*Opm::unit::barsa, 0, 0, wells.get());
// ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
double oildistr[2] = { 0.0, 1.0 };
ok = ok && append_well_controls(SURFACE_RATE, 1e-3, oildistr, 1, wells.get());
if (!ok) {
OPM_THROW(std::runtime_error, "Something went wrong with well init.");
}
set_current_control(0, 0, wells.get());
set_current_control(1, 0, wells.get());
return wells;
}
template <class Ostream, typename T, class A>
Ostream&
operator<<(Ostream& os, const std::vector<T,A>& v)
{
std::copy(v.begin(), v.end(), std::ostream_iterator<T>(os, " "));
return os;
}
}
int
main(int argc, char* argv[])
try
{
const Opm::parameter::ParameterGroup param(argc, argv, false);
const Opm::GridManager gm(20, 1);
const UnstructuredGrid* g = gm.c_grid();
const int nc = g->number_of_cells;
const Opm::BlackoilPropertiesBasic props0(param, 2, nc);
const Opm::BlackoilPropsAd props(props0);
std::shared_ptr<Wells> wells = createWellConfig();
double grav[] = { 0.0, 0.0 };
Opm::DerivedGeology geo(*g, props, grav);
Opm::NewtonIterationBlackoilSimple fis_solver(param);
Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(param, *g, props, geo, 0, *wells, fis_solver, /*hasDisgas*/ true, /*hasVapoil=*/false);
Opm::BlackoilState state;
initStateBasic(*g, props0, param, 0.0, state);
initBlackoilSurfvol(*g, props0, state);
Opm::WellStateFullyImplicitBlackoil well_state;
well_state.init(wells.get(), state);
solver.step(1.0, state, well_state);
std::cout << state.pressure() << '\n'
<< well_state.bhp() << '\n';
return 0;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

9
opm/autodiff/GeoProps.hpp
View File

@ -25,6 +25,8 @@
//#include <opm/core/pressure/tpfa/trans_tpfa.h> //#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/tpfa/TransTpfa.hpp> #include <opm/core/pressure/tpfa/TransTpfa.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include "disable_warning_pragmas.h" #include "disable_warning_pragmas.h"
#include <Eigen/Eigen> #include <Eigen/Eigen>
@ -49,9 +51,10 @@ namespace Opm
/// Construct contained derived geological properties /// Construct contained derived geological properties
/// from grid and property information. /// from grid and property information.
template <class Props, class Grid> template <class Props, class Grid>
DerivedGeology(const Grid& grid, DerivedGeology(const Grid& grid,
const Props& props , const Props& props ,
const double* grav = 0) Opm::EclipseStateConstPtr eclState,
const double* grav = 0)
: pvol_ (Opm::AutoDiffGrid::numCells(grid)) : pvol_ (Opm::AutoDiffGrid::numCells(grid))
, trans_(Opm::AutoDiffGrid::numFaces(grid)) , trans_(Opm::AutoDiffGrid::numFaces(grid))
, gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1)) , gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))