mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
remove historical and un-used files.
This commit is contained in:
Liu Ming
parent
c0a61c9655
commit
a23e4ca63b
@ -1,262 +0,0 @@
|
||||
/*
|
||||
*/
|
||||
|
||||
#include <opm/core/pressure/FlowBCManager.hpp>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/wells.h>
|
||||
#include <opm/core/wells/WellsManager.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/simulator/initState.hpp>
|
||||
#include <opm/core/simulator/SimulatorReport.hpp>
|
||||
#include <opm/core/simulator/SimulatorTimer.hpp>
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
#include <opm/core/props/IncompPropertiesBasic.hpp>
|
||||
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
|
||||
#include <opm/core/props/rock/RockCompressibility.hpp>
|
||||
|
||||
#include <opm/core/linalg/LinearSolverFactory.hpp>
|
||||
|
||||
#include <opm/core/simulator/TwophaseState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
#include <opm/autodiff/polymer/SimulatorFullyImplicitTwophase.hpp>
|
||||
#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
|
||||
#include <opm/autodiff/polymer/IncompPropsAdBasic.hpp>
|
||||
#include <opm/autodiff/polymer/IncompPropsAdFromDeck.hpp>
|
||||
|
||||
#include <boost/scoped_ptr.hpp>
|
||||
#include <boost/filesystem.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <numeric>
|
||||
|
||||
|
||||
namespace
|
||||
{
|
||||
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
|
||||
{
|
||||
if (param.anyUnused()) {
|
||||
std::cout << "-------------------- Unused parameters: --------------------\n";
|
||||
param.displayUsage();
|
||||
std::cout << "----------------------------------------------------------------" << std::endl;
|
||||
}
|
||||
}
|
||||
} // anon namespace
|
||||
|
||||
|
||||
|
||||
// ----------------- Main program -----------------
|
||||
int
|
||||
main(int argc, char** argv)
|
||||
try
|
||||
{
|
||||
using namespace Opm;
|
||||
|
||||
std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n";
|
||||
parameter::ParameterGroup param(argc, argv, false);
|
||||
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");
|
||||
boost::scoped_ptr<EclipseGridParser> deck;
|
||||
boost::scoped_ptr<GridManager> grid;
|
||||
boost::scoped_ptr<IncompPropsAdInterface> props;
|
||||
TwophaseState state;
|
||||
double gravity[3] = { 0.0 };
|
||||
if (use_deck) {
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
deck.reset(new EclipseGridParser(deck_filename));
|
||||
// Grid init
|
||||
grid.reset(new GridManager(*deck));
|
||||
// Rock and fluid init
|
||||
props.reset(new IncompPropsAdFromDeck(*deck, *grid->c_grid()));
|
||||
// Gravity.
|
||||
gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
|
||||
// Init state variables (saturation and pressure).
|
||||
int num_cells = grid->c_grid()->number_of_cells;
|
||||
if (param.has("init_saturation")) {
|
||||
//initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
|
||||
const double init_saturation = param.get<double>("init_saturation");
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
state.saturation()[2*c] = init_saturation;
|
||||
state.saturation()[2*c+1] = 1. - init_saturation;
|
||||
}
|
||||
} else {
|
||||
if (deck->hasField("PRESSURE")) {
|
||||
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
|
||||
std::vector<double>& s = state.saturation();
|
||||
std::vector<double>& p = state.pressure();
|
||||
const std::vector<double>& p_deck = deck->getFloatingPointValue("PRESSURE");
|
||||
// water-oil or water-gas: we require SWAT
|
||||
if (!deck->hasField("SWAT")) {
|
||||
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
|
||||
}
|
||||
const std::vector<double>& sw_deck = deck->getFloatingPointValue("SWAT");
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
int c_deck = (grid->c_grid()->global_cell == NULL) ? c : grid->c_grid()->global_cell[c];
|
||||
s[2*c] = sw_deck[c_deck];
|
||||
s[2*c + 1] = 1.0 - sw_deck[c_deck];
|
||||
p[c] = p_deck[c_deck];
|
||||
}
|
||||
}
|
||||
}
|
||||
} 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 IncompPropsAdBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
|
||||
// Rock compressibility.
|
||||
// Gravity.
|
||||
gravity[2] = param.getDefault("gravity", 0.0);
|
||||
int num_cells = grid->c_grid()->number_of_cells;
|
||||
}
|
||||
|
||||
// Warn if gravity but no density difference.
|
||||
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
|
||||
const double *grav = use_gravity ? &gravity[0] : 0;
|
||||
|
||||
// Initialising src
|
||||
std::vector<double> src(num_cells, 0.0);
|
||||
if (use_deck) {
|
||||
// Do nothing, wells will be the driving force, not source terms.
|
||||
if (deck->hasField("SRC")) {
|
||||
const std::vector<double>& src_deck = deck->getFloatingPointValue("SRC");
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
int c_deck = (grid->c_grid()->global_cell == NULL) ? c : grid->c_grid()->global_cell[c];
|
||||
src[c] = src_deck[c_deck];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Compute pore volumes, in order to enable specifying injection rate
|
||||
// terms of total pore volume.
|
||||
std::vector<double> porevol;
|
||||
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;
|
||||
src[0] = flow_per_sec;
|
||||
src[num_cells - 1] = -flow_per_sec;
|
||||
}
|
||||
|
||||
in: num_cells = grid->c_grid()->number_of_cells;
|
||||
|
||||
// Linear solver.
|
||||
LinearSolverFactory linsolver(param);
|
||||
|
||||
// Write parameters used for later reference.
|
||||
bool output = param.getDefault("output", true);
|
||||
std::ofstream epoch_os;
|
||||
std::string output_dir;
|
||||
if (output) {
|
||||
output_dir =
|
||||
param.getDefault("output_dir", std::string("output"));
|
||||
boost::filesystem::path fpath(output_dir);
|
||||
try {
|
||||
create_directories(fpath);
|
||||
}
|
||||
catch (...) {
|
||||
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
|
||||
}
|
||||
std::string filename = output_dir + "/epoch_timing.param";
|
||||
epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
|
||||
// open file to clean it. The file is appended to in SimulatorTwophase
|
||||
filename = output_dir + "/step_timing.param";
|
||||
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
|
||||
step_os.close();
|
||||
param.writeParam(output_dir + "/simulation.param");
|
||||
}
|
||||
|
||||
|
||||
std::cout << "\n\n================ Starting main simulation loop ===============\n"
|
||||
<< " (number of epochs: "
|
||||
<< (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush;
|
||||
|
||||
SimulatorReport rep;
|
||||
if (!use_deck) {
|
||||
// Simple simulation without a deck.
|
||||
SimulatorFullyImplicitTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
linsolver,
|
||||
src);
|
||||
SimulatorTimer simtimer;
|
||||
simtimer.init(param);
|
||||
warnIfUnusedParams(param);
|
||||
rep = simulator.run(simtimer, state, src);
|
||||
} else {
|
||||
// With a deck, we may have more epochs etc.
|
||||
int step = 0;
|
||||
SimulatorTimer simtimer;
|
||||
// Use timer for last epoch to obtain total time.
|
||||
deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
|
||||
simtimer.init(*deck);
|
||||
const double total_time = simtimer.totalTime();
|
||||
for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
|
||||
// Set epoch index.
|
||||
deck->setCurrentEpoch(epoch);
|
||||
|
||||
// Update the timer.
|
||||
if (deck->hasField("TSTEP")) {
|
||||
simtimer.init(*deck);
|
||||
} else {
|
||||
if (epoch != 0) {
|
||||
OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch);
|
||||
}
|
||||
simtimer.init(param);
|
||||
}
|
||||
simtimer.setCurrentStepNum(step);
|
||||
simtimer.setTotalTime(total_time);
|
||||
|
||||
// Report on start of epoch.
|
||||
std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------"
|
||||
<< "\n (number of steps: "
|
||||
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
|
||||
|
||||
|
||||
// Create and run simulator.
|
||||
SimulatorFullyImplicitTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
linsolver,
|
||||
src);
|
||||
if (epoch == 0) {
|
||||
warnIfUnusedParams(param);
|
||||
}
|
||||
SimulatorReport epoch_rep = simulator.run(simtimer, state, src);
|
||||
if (output) {
|
||||
epoch_rep.reportParam(epoch_os);
|
||||
}
|
||||
// Update total timing report and remember step number.
|
||||
rep += epoch_rep;
|
||||
step = simtimer.currentStepNum();
|
||||
}
|
||||
}
|
||||
|
||||
std::cout << "\n\n================ End of simulation ===============\n\n";
|
||||
rep.report(std::cout);
|
||||
|
||||
if (output) {
|
||||
std::string filename = output_dir + "/walltime.param";
|
||||
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
|
||||
rep.reportParam(tot_os);
|
||||
}
|
||||
|
||||
}
|
||||
catch (const std::exception &e) {
|
||||
std::cerr << "Program threw an exception: " << e.what() << "\n";
|
||||
throw;
|
||||
}
|
||||
|
||||
@ -1,115 +0,0 @@
|
||||
#include "config.h"
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
#include <cassert>
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
|
||||
#include <opm/core/wells.h>
|
||||
#include <opm/core/wells/WellsManager.hpp>
|
||||
#include <opm/core/io/vtk/writeVtkData.hpp>
|
||||
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
|
||||
#include <opm/core/pressure/FlowBCManager.hpp>
|
||||
#include <opm/core/props/IncompPropertiesBasic.hpp>
|
||||
#include <opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
|
||||
|
||||
|
||||
#include <opm/core/simulator/TwophaseState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
int main (int argc, char** argv)
|
||||
try
|
||||
{
|
||||
|
||||
using namespace Opm;
|
||||
parameter::ParameterGroup param(argc, argv, false);
|
||||
|
||||
bool use_deck = param.has("deck_filename");
|
||||
if (!use_deck) {
|
||||
OPM_THROW(std::runtime_error, "FullyImplicitTwoPhaseSolver cannot run without deckfile.");
|
||||
}
|
||||
double gravity[3] = { 0.0 };
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
EclipseGridParser deck = EclipseGridParser(deck_filename);
|
||||
int nx = param.getDefault("nx", 30);
|
||||
int ny = param.getDefault("ny", 1);
|
||||
int nz = 1;
|
||||
double dx = 10./nx;
|
||||
double dy = 1.0;
|
||||
double dz = 1.0;
|
||||
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
|
||||
const UnstructuredGrid& grid = *grid_manager.c_grid();
|
||||
int num_cells = grid.number_of_cells;
|
||||
int num_phases = 2;
|
||||
using namespace Opm::unit;
|
||||
using namespace Opm::prefix;
|
||||
std::vector<double> density(num_phases, 1000.0);
|
||||
std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
|
||||
viscosity[0] = 0.5 * centi * Poise;
|
||||
viscosity[1] = 5 * centi * Poise;
|
||||
double porosity = 0.35;
|
||||
double permeability = 10.0*milli*darcy;
|
||||
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
|
||||
IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity,
|
||||
porosity, permeability, grid.dimensions, num_cells);
|
||||
/*
|
||||
std::vector<double> src(num_cells, 0.0);
|
||||
src[0] = 1. / day;
|
||||
src[num_cells-1] = -1. / day;
|
||||
*/
|
||||
|
||||
FlowBCManager bcs;
|
||||
LinearSolverUmfpack linsolver;
|
||||
TwophaseState state;
|
||||
state.init(grid, 2);
|
||||
WellState well_state;
|
||||
WellsManager wells(deck, grid, props.permeability());
|
||||
well_state.init(wells.c_wells(), state);
|
||||
FullyImplicitTwoPhaseSolver solver(grid, props, *wells.c_wells(), linsolver, gravity);
|
||||
std::vector<double> porevol;
|
||||
Opm::computePorevolume(grid, props.porosity(), porevol);
|
||||
const double dt = param.getDefault("dt", 10.) * day;
|
||||
const int num_time_steps = param.getDefault("nsteps", 10);
|
||||
std::vector<int> allcells(num_cells);
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
allcells[cell] = cell;
|
||||
}
|
||||
std::vector<double> src; // empty src term.
|
||||
gravity[2] = param.getDefault("gravity", 0.0);
|
||||
//initial sat
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
state.saturation()[2*c] = 0.2;
|
||||
state.saturation()[2*c+1] = 0.8;
|
||||
}
|
||||
std::vector<double> p(num_cells, 100*Opm::unit::barsa);
|
||||
state.pressure() = p;
|
||||
std::ostringstream vtkfilename;
|
||||
vtkfilename.str("");
|
||||
vtkfilename << "sim_2p_fincomp_" << std::setw(3) << std::setfill('0') << 0 << ".vtu";
|
||||
std::ofstream vtkfile(vtkfilename.str().c_str());
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
Opm::writeVtkData(grid, dm, vtkfile);
|
||||
for (int i = 0; i < num_time_steps; ++i) {
|
||||
solver.step(dt, state, src, well_state);
|
||||
vtkfilename.str("");
|
||||
vtkfilename << "sim_2p_fincomp_" << std::setw(3) << std::setfill('0') << i + 1 << ".vtu";
|
||||
std::ofstream vtkfile(vtkfilename.str().c_str());
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
Opm::writeVtkData(grid, dm, vtkfile);
|
||||
}
|
||||
}
|
||||
catch (const std::exception &e) {
|
||||
std::cerr << "Program threw an exception: " << e.what() << "\n";
|
||||
throw;
|
||||
}
|
||||
@ -1,242 +0,0 @@
|
||||
/*
|
||||
Copyright 2013 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/core/pressure/FlowBCManager.hpp>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/wells.h>
|
||||
#include <opm/core/wells/WellsManager.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/simulator/initState.hpp>
|
||||
#include <opm/core/simulator/SimulatorReport.hpp>
|
||||
#include <opm/core/simulator/SimulatorTimer.hpp>
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
#include <opm/core/io/eclipse/EclipseWriter.hpp>
|
||||
#include <opm/core/props/IncompPropertiesBasic.hpp>
|
||||
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
|
||||
#include <opm/core/props/rock/RockCompressibility.hpp>
|
||||
|
||||
#include <opm/core/linalg/LinearSolverFactory.hpp>
|
||||
|
||||
#include <opm/core/simulator/TwophaseState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
|
||||
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp>
|
||||
#include <opm/core/utility/share_obj.hpp>
|
||||
|
||||
#include <boost/scoped_ptr.hpp>
|
||||
#include <boost/filesystem.hpp>
|
||||
#include <boost/algorithm/string.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <numeric>
|
||||
|
||||
|
||||
namespace
|
||||
{
|
||||
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
|
||||
{
|
||||
if (param.anyUnused()) {
|
||||
std::cout << "-------------------- Unused parameters: --------------------\n";
|
||||
param.displayUsage();
|
||||
std::cout << "----------------------------------------------------------------" << std::endl;
|
||||
}
|
||||
}
|
||||
} // anon namespace
|
||||
|
||||
|
||||
|
||||
// ----------------- Main program -----------------
|
||||
int
|
||||
main(int argc, char** argv)
|
||||
try
|
||||
{
|
||||
using namespace Opm;
|
||||
|
||||
std::cout << "\n================ Test program for fully implicit three-phase black-oil flow ===============\n\n";
|
||||
parameter::ParameterGroup param(argc, argv, false);
|
||||
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) {
|
||||
OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
|
||||
"Specify the deck with deck_filename=deckname.data (for example).");
|
||||
}
|
||||
boost::scoped_ptr<EclipseGridParser> deck;
|
||||
boost::scoped_ptr<GridManager> grid;
|
||||
boost::scoped_ptr<IncompPropertiesInterface> props;
|
||||
boost::scoped_ptr<IncompPropsAdInterface> new_props;
|
||||
TwophaseState state;
|
||||
// bool check_well_controls = false;
|
||||
// int max_well_control_iterations = 0;
|
||||
double gravity[3] = { 0.0 };
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
deck.reset(new EclipseGridParser(deck_filename));
|
||||
// Grid init
|
||||
grid.reset(new GridManager(*deck));
|
||||
|
||||
// use the capitalized part of the deck's filename between the
|
||||
// last '/' and the last '.' character as base name.
|
||||
std::string baseName = deck_filename;
|
||||
auto charPos = baseName.rfind('/');
|
||||
if (charPos != std::string::npos)
|
||||
baseName = baseName.substr(charPos + 1);
|
||||
charPos = baseName.rfind('.');
|
||||
if (charPos != std::string::npos)
|
||||
baseName = baseName.substr(0, charPos);
|
||||
baseName = boost::to_upper_copy(baseName);
|
||||
|
||||
// Rock and fluid init
|
||||
props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
|
||||
new_props.reset(new IncompPropsAdFromDeck(*deck, *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.
|
||||
// Gravity.
|
||||
gravity[2] = deck->hasField("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);
|
||||
}
|
||||
|
||||
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
|
||||
const double *grav = use_gravity ? &gravity[0] : 0;
|
||||
|
||||
// Linear solver.
|
||||
LinearSolverFactory linsolver(param);
|
||||
|
||||
// Write parameters used for later reference.
|
||||
bool output = param.getDefault("output", true);
|
||||
std::ofstream epoch_os;
|
||||
std::string output_dir;
|
||||
if (output) {
|
||||
output_dir =
|
||||
param.getDefault("output_dir", std::string("output"));
|
||||
boost::filesystem::path fpath(output_dir);
|
||||
try {
|
||||
create_directories(fpath);
|
||||
}
|
||||
catch (...) {
|
||||
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
|
||||
}
|
||||
std::string filename = output_dir + "/epoch_timing.param";
|
||||
epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
|
||||
// open file to clean it. The file is appended to in SimulatorTwophase
|
||||
filename = output_dir + "/step_timing.param";
|
||||
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
|
||||
step_os.close();
|
||||
param.writeParam(output_dir + "/simulation.param");
|
||||
}
|
||||
|
||||
|
||||
std::cout << "\n\n================ Starting main simulation loop ===============\n"
|
||||
<< " (number of epochs: "
|
||||
<< (deck->numberOfEpochs()) << ")\n\n" << std::flush;
|
||||
|
||||
SimulatorReport rep;
|
||||
// With a deck, we may have more epochs etc.
|
||||
WellState well_state;
|
||||
int step = 0;
|
||||
SimulatorTimer simtimer;
|
||||
// Use timer for last epoch to obtain total time.
|
||||
deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
|
||||
simtimer.init(*deck);
|
||||
const double total_time = simtimer.totalTime();
|
||||
for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
|
||||
// Set epoch index.
|
||||
deck->setCurrentEpoch(epoch);
|
||||
|
||||
// Update the timer.
|
||||
if (deck->hasField("TSTEP")) {
|
||||
simtimer.init(*deck);
|
||||
} else {
|
||||
if (epoch != 0) {
|
||||
OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch);
|
||||
}
|
||||
simtimer.init(param);
|
||||
}
|
||||
simtimer.setCurrentStepNum(step);
|
||||
simtimer.setTotalTime(total_time);
|
||||
|
||||
// Report on start of epoch.
|
||||
std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------"
|
||||
<< "\n (number of steps: "
|
||||
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
|
||||
|
||||
// Create new wells, well_state
|
||||
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
|
||||
// @@@ HACK: we should really make a new well state and
|
||||
// properly transfer old well state to it every epoch,
|
||||
// since number of wells may change etc.
|
||||
if (epoch == 0) {
|
||||
well_state.init(wells.c_wells(), state);
|
||||
}
|
||||
|
||||
// Create and run simulator.
|
||||
std::vector<double> src(grid->c_grid()->number_of_cells, 0.0);
|
||||
src[0] = 10. / Opm::unit::day;
|
||||
src[grid->c_grid()->number_of_cells-1] = -10. / Opm::unit::day;
|
||||
SimulatorFullyImplicitTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*new_props,
|
||||
wells,
|
||||
linsolver,
|
||||
src,
|
||||
grav);
|
||||
if (epoch == 0) {
|
||||
warnIfUnusedParams(param);
|
||||
}
|
||||
SimulatorReport epoch_rep = simulator.run(simtimer, state, src, well_state);
|
||||
if (output) {
|
||||
epoch_rep.reportParam(epoch_os);
|
||||
}
|
||||
// Update total timing report and remember step number.
|
||||
rep += epoch_rep;
|
||||
step = simtimer.currentStepNum();
|
||||
}
|
||||
|
||||
std::cout << "\n\n================ End of simulation ===============\n\n";
|
||||
rep.report(std::cout);
|
||||
|
||||
if (output) {
|
||||
std::string filename = output_dir + "/walltime.param";
|
||||
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
|
||||
rep.reportParam(tot_os);
|
||||
}
|
||||
|
||||
}
|
||||
catch (const std::exception &e) {
|
||||
std::cerr << "Program threw an exception: " << e.what() << "\n";
|
||||
throw;
|
||||
}
|
||||
|
||||
@ -1,157 +0,0 @@
|
||||
#include "config.h"
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
#include <cassert>
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/io/vtk/writeVtkData.hpp>
|
||||
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
|
||||
#include <opm/core/pressure/FlowBCManager.hpp>
|
||||
#include <opm/core/props/IncompPropertiesBasic.hpp>
|
||||
#include <opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp>
|
||||
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
|
||||
#include <opm/polymer/PolymerState.hpp>
|
||||
#include <opm/polymer/PolymerInflow.hpp>
|
||||
#include <opm/polymer/PolymerProperties.hpp>
|
||||
|
||||
#include <opm/core/simulator/TwophaseState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
|
||||
int main (int argc, char** argv)
|
||||
try
|
||||
{
|
||||
using namespace Opm;
|
||||
parameter::ParameterGroup param(argc, argv, false);
|
||||
bool use_poly_deck = param.has("deck_filename");
|
||||
if (!use_poly_deck) {
|
||||
OPM_THROW(std::runtime_error, "Polymer Properties must be read from deck_filename\n");
|
||||
}
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
EclipseGridParser deck = EclipseGridParser(deck_filename);
|
||||
int nx = param.getDefault("nx", 30);
|
||||
int ny = param.getDefault("ny", 30);
|
||||
int nz = 1;
|
||||
double dx = 10./nx;
|
||||
double dy = 1.0;
|
||||
double dz = 1.0;
|
||||
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
|
||||
const UnstructuredGrid& grid = *grid_manager.c_grid();
|
||||
int num_cells = grid.number_of_cells;
|
||||
int num_phases = 2;
|
||||
using namespace Opm::unit;
|
||||
using namespace Opm::prefix;
|
||||
std::vector<double> density(num_phases, 1000.0);
|
||||
std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
|
||||
viscosity[0] = 0.5 * centi * Poise;
|
||||
viscosity[1] = 5 * centi * Poise;
|
||||
double porosity = 0.35;
|
||||
double permeability = 10.0*milli*darcy;
|
||||
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
|
||||
IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity,
|
||||
porosity, permeability, grid.dimensions, num_cells);
|
||||
|
||||
// Init polymer properties.
|
||||
// Setting defaults to provide a simple example case.
|
||||
PolymerProperties polymer_props(deck);
|
||||
#if 0
|
||||
if (use_poly_deck) {
|
||||
} else {
|
||||
double c_max = param.getDefault("c_max_limit", 5.0);
|
||||
double mix_param = param.getDefault("mix_param", 1.0);
|
||||
double rock_density = param.getDefault("rock_density", 1000.0);
|
||||
double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15);
|
||||
double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability
|
||||
double c_max_ads = param.getDefault("c_max_ads", 1.);
|
||||
int ads_index = param.getDefault<int>("ads_index", Opm::PolymerProperties::NoDesorption);
|
||||
std::vector<double> c_vals_visc(2, -1e100);
|
||||
c_vals_visc[0] = 0.0;
|
||||
c_vals_visc[1] = 7.0;
|
||||
std::vector<double> visc_mult_vals(2, -1e100);
|
||||
visc_mult_vals[0] = 1.0;
|
||||
// poly_props.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
|
||||
visc_mult_vals[1] = 20.0;
|
||||
std::vector<double> c_vals_ads(3, -1e100);
|
||||
c_vals_ads[0] = 0.0;
|
||||
c_vals_ads[1] = 2.0;
|
||||
c_vals_ads[2] = 8.0;
|
||||
std::vector<double> ads_vals(3, -1e100);
|
||||
ads_vals[0] = 0.0;
|
||||
ads_vals[1] = 0.0015;
|
||||
ads_vals[2] = 0.0025;
|
||||
PolymerProperties polymer_props;
|
||||
polymer_props.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads,
|
||||
static_cast<Opm::PolymerProperties::AdsorptionBehaviour>(ads_index),
|
||||
c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals);
|
||||
}
|
||||
#endif
|
||||
PolymerPropsAd polymer_props_ad(polymer_props);
|
||||
std::vector<double> omega;
|
||||
std::vector<double> src(num_cells, 0.0);
|
||||
std::vector<double> src_polymer(num_cells);
|
||||
src[0] = param.getDefault("insrc", 1.) / day;
|
||||
src[num_cells-1] = -param.getDefault("insrc", 1.) / day;
|
||||
|
||||
PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
|
||||
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
|
||||
param.getDefault("poly_amount", polymer_props.cMax()));
|
||||
FlowBCManager bcs;
|
||||
LinearSolverUmfpack linsolver;
|
||||
FullyImplicitTwophasePolymerSolver solver(grid, props,polymer_props_ad, linsolver);
|
||||
std::vector<double> porevol;
|
||||
Opm::computePorevolume(grid, props.porosity(), porevol);
|
||||
const double dt = param.getDefault("dt", 10.) * day;
|
||||
const int num_time_steps = param.getDefault("nsteps", 10);
|
||||
std::vector<int> allcells(num_cells);
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
allcells[cell] = cell;
|
||||
}
|
||||
PolymerState state;
|
||||
state.init(grid, 2);
|
||||
//initial sat
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
state.saturation()[2*c] = 0.2;
|
||||
state.saturation()[2*c+1] = 0.8;
|
||||
}
|
||||
std::vector<double> p(num_cells, 100*Opm::unit::barsa);
|
||||
state.pressure() = p;
|
||||
|
||||
std::vector<double> c(num_cells, 0.0);
|
||||
state.concentration() = c;
|
||||
std::ostringstream vtkfilename;
|
||||
double currentime = 0;
|
||||
|
||||
// Write the initial state.
|
||||
vtkfilename.str("");
|
||||
vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << 0<< ".vtu";
|
||||
std::ofstream vtkfile(vtkfilename.str().c_str());
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
dm["concentration"] = &state.concentration();
|
||||
Opm::writeVtkData(grid, dm, vtkfile);
|
||||
for (int i = 0; i < num_time_steps; ++i) {
|
||||
currentime += dt;
|
||||
polymer_inflow.getInflowValues(currentime, currentime+dt, src_polymer);
|
||||
solver.step(dt, state, src, src_polymer);
|
||||
vtkfilename.str("");
|
||||
vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << i + 1<< ".vtu";
|
||||
std::ofstream vtkfile(vtkfilename.str().c_str());
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
dm["concentration"] = &state.concentration();
|
||||
Opm::writeVtkData(grid, dm, vtkfile);
|
||||
}
|
||||
}
|
||||
catch (const std::exception &e) {
|
||||
std::cerr << "Program threw an exception: " << e.what() << "\n";
|
||||
throw;
|
||||
}
|
||||
@ -1,274 +0,0 @@
|
||||
/*
|
||||
Copyright 2013 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/core/pressure/FlowBCManager.hpp>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/grid/GridManager.hpp>
|
||||
#include <opm/core/wells.h>
|
||||
#include <opm/core/wells/WellsManager.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/simulator/initState.hpp>
|
||||
#include <opm/core/simulator/SimulatorReport.hpp>
|
||||
#include <opm/core/simulator/SimulatorTimer.hpp>
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
#include <opm/core/io/eclipse/EclipseWriter.hpp>
|
||||
#include <opm/core/props/IncompPropertiesBasic.hpp>
|
||||
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
|
||||
#include <opm/core/props/rock/RockCompressibility.hpp>
|
||||
|
||||
#include <opm/core/linalg/LinearSolverFactory.hpp>
|
||||
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
|
||||
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
|
||||
#include <opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp>
|
||||
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
|
||||
#include <opm/polymer/PolymerProperties.hpp>
|
||||
#include <opm/polymer/PolymerInflow.hpp>
|
||||
#include <opm/polymer/PolymerState.hpp>
|
||||
#include <opm/core/utility/share_obj.hpp>
|
||||
|
||||
#include <boost/scoped_ptr.hpp>
|
||||
#include <boost/filesystem.hpp>
|
||||
#include <boost/algorithm/string.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <numeric>
|
||||
|
||||
|
||||
namespace
|
||||
{
|
||||
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
|
||||
{
|
||||
if (param.anyUnused()) {
|
||||
std::cout << "-------------------- Unused parameters: --------------------\n";
|
||||
param.displayUsage();
|
||||
std::cout << "----------------------------------------------------------------" << std::endl;
|
||||
}
|
||||
}
|
||||
} // anon namespace
|
||||
|
||||
|
||||
|
||||
// ----------------- Main program -----------------
|
||||
int
|
||||
main(int argc, char** argv)
|
||||
try
|
||||
{
|
||||
using namespace Opm;
|
||||
|
||||
std::cout << "\n================ Test program for fully implicit three-phase black-oil flow ===============\n\n";
|
||||
parameter::ParameterGroup param(argc, argv, false);
|
||||
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) {
|
||||
OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
|
||||
"Specify the deck with deck_filename=deckname.data (for example).");
|
||||
}
|
||||
boost::scoped_ptr<EclipseGridParser> deck;
|
||||
boost::scoped_ptr<GridManager> grid;
|
||||
boost::scoped_ptr<IncompPropertiesInterface> props;
|
||||
boost::scoped_ptr<IncompPropsAdInterface> new_props;
|
||||
// boost::scoped_ptr<PolymerPropsAd> polymer_props;
|
||||
PolymerState state;
|
||||
// bool check_well_controls = false;
|
||||
// int max_well_control_iterations = 0;
|
||||
double gravity[3] = { 0.0 };
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
deck.reset(new EclipseGridParser(deck_filename));
|
||||
// Grid init
|
||||
grid.reset(new GridManager(*deck));
|
||||
|
||||
// use the capitalized part of the deck's filename between the
|
||||
// last '/' and the last '.' character as base name.
|
||||
std::string baseName = deck_filename;
|
||||
auto charPos = baseName.rfind('/');
|
||||
if (charPos != std::string::npos)
|
||||
baseName = baseName.substr(charPos + 1);
|
||||
charPos = baseName.rfind('.');
|
||||
if (charPos != std::string::npos)
|
||||
baseName = baseName.substr(0, charPos);
|
||||
baseName = boost::to_upper_copy(baseName);
|
||||
|
||||
// Rock and fluid init
|
||||
props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
|
||||
new_props.reset(new IncompPropsAdFromDeck(*deck, *grid->c_grid()));
|
||||
PolymerProperties polymer_props(*deck);
|
||||
PolymerPropsAd polymer_props_ad(polymer_props);
|
||||
// polymer_props.reset(new PolymerPropsAd(*deck, *grid->c_grid()));
|
||||
// check_well_controls = param.getDefault("check_well_controls", false);
|
||||
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
|
||||
// Gravity.
|
||||
gravity[2] = deck->hasField("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);
|
||||
}
|
||||
|
||||
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
|
||||
const double* grav = use_gravity ? &gravity[0] : 0;
|
||||
|
||||
// Linear solver.
|
||||
LinearSolverFactory linsolver(param);
|
||||
|
||||
// Write parameters used for later reference.
|
||||
bool output = param.getDefault("output", true);
|
||||
std::ofstream epoch_os;
|
||||
std::string output_dir;
|
||||
if (output) {
|
||||
output_dir =
|
||||
param.getDefault("output_dir", std::string("output"));
|
||||
boost::filesystem::path fpath(output_dir);
|
||||
try {
|
||||
create_directories(fpath);
|
||||
}
|
||||
catch (...) {
|
||||
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
|
||||
}
|
||||
std::string filename = output_dir + "/epoch_timing.param";
|
||||
epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
|
||||
// open file to clean it. The file is appended to in SimulatorTwophase
|
||||
filename = output_dir + "/step_timing.param";
|
||||
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
|
||||
step_os.close();
|
||||
param.writeParam(output_dir + "/simulation.param");
|
||||
}
|
||||
|
||||
|
||||
std::cout << "\n\n================ Starting main simulation loop ===============\n"
|
||||
<< " (number of epochs: "
|
||||
<< (deck->numberOfEpochs()) << ")\n\n" << std::flush;
|
||||
|
||||
SimulatorReport rep;
|
||||
// With a deck, we may have more epochs etc.
|
||||
WellState well_state;
|
||||
int step = 0;
|
||||
SimulatorTimer simtimer;
|
||||
// Use timer for last epoch to obtain total time.
|
||||
deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
|
||||
simtimer.init(*deck);
|
||||
const double total_time = simtimer.totalTime();
|
||||
// Check for WPOLYMER presence in last epoch to decide
|
||||
// polymer injection control type.
|
||||
const bool use_wpolymer = deck->hasField("WPOLYMER");
|
||||
if (use_wpolymer) {
|
||||
if (param.has("poly_start_days")) {
|
||||
OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. "
|
||||
"You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
|
||||
}
|
||||
}
|
||||
for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
|
||||
// Set epoch index.
|
||||
deck->setCurrentEpoch(epoch);
|
||||
|
||||
// Update the timer.
|
||||
if (deck->hasField("TSTEP")) {
|
||||
simtimer.init(*deck);
|
||||
} else {
|
||||
if (epoch != 0) {
|
||||
OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch);
|
||||
}
|
||||
simtimer.init(param);
|
||||
}
|
||||
simtimer.setCurrentStepNum(step);
|
||||
simtimer.setTotalTime(total_time);
|
||||
|
||||
// Report on start of epoch.
|
||||
std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------"
|
||||
<< "\n (number of steps: "
|
||||
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
|
||||
|
||||
// Create new wells, polymer inflow controls.
|
||||
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
|
||||
boost::scoped_ptr<PolymerInflowInterface> polymer_inflow;
|
||||
if (use_wpolymer) {
|
||||
if (wells.c_wells() == 0) {
|
||||
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
|
||||
}
|
||||
polymer_inflow.reset(new PolymerInflowFromDeck(*deck, *wells.c_wells(), props->numCells()));
|
||||
} else {
|
||||
polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
|
||||
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
|
||||
param.getDefault("poly_amount", polymer_props.cMax())));
|
||||
}
|
||||
// @@@ HACK: we should really make a new well state and
|
||||
// properly transfer old well state to it every epoch,
|
||||
// since number of wells may change etc.
|
||||
if (epoch == 0) {
|
||||
well_state.init(wells.c_wells(), state);
|
||||
}
|
||||
|
||||
// Create and run simulator.
|
||||
#if 0
|
||||
std::vector<double> src(grid->c_grid()->number_of_cells, 0.0);
|
||||
src[0] = 10. / Opm::unit::day;
|
||||
src[grid->c_grid()->number_of_cells-1] = -10. / Opm::unit::day;
|
||||
PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
|
||||
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
|
||||
param.getDefault("poly_amount", polymer_props.cMax()));
|
||||
#endif
|
||||
SimulatorFullyImplicitTwophasePolymer simulator(param,
|
||||
*grid->c_grid(),
|
||||
*new_props,
|
||||
polymer_props_ad,
|
||||
linsolver,
|
||||
wells,
|
||||
*polymer_inflow,
|
||||
grav);
|
||||
if (epoch == 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;
|
||||
step = simtimer.currentStepNum();
|
||||
}
|
||||
|
||||
std::cout << "\n\n================ End of simulation ===============\n\n";
|
||||
rep.report(std::cout);
|
||||
|
||||
if (output) {
|
||||
std::string filename = output_dir + "/walltime.param";
|
||||
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
|
||||
rep.reportParam(tot_os);
|
||||
}
|
||||
|
||||
}
|
||||
catch (const std::exception &e) {
|
||||
std::cerr << "Program threw an exception: " << e.what() << "\n";
|
||||
throw;
|
||||
}
|
||||
|
||||
Binary file not shown.
Loading…
Reference in New Issue
Block a user