mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge branch 'tester' into master-refactor-for-cpgrid-support
Resolved conflicts: examples/sim_fibo_ad.cpp opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
This commit is contained in:
Markus Blatt
commit
327f1b003d
@ -99,51 +99,25 @@ try
|
||||
double gravity[3] = { 0.0 };
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
|
||||
#define USE_NEW_PARSER 1
|
||||
#if USE_NEW_PARSER
|
||||
Opm::ParserPtr newParser(new Opm::Parser() );
|
||||
Opm::DeckConstPtr newParserDeck = newParser->parseFile( deck_filename );
|
||||
#else
|
||||
std::shared_ptr<EclipseGridParser> deck;
|
||||
deck.reset(new EclipseGridParser(deck_filename));
|
||||
#endif
|
||||
|
||||
// Grid init
|
||||
#if USE_NEW_PARSER
|
||||
grid.reset(new GridManager(newParserDeck));
|
||||
#else
|
||||
grid.reset(new GridManager(*deck));
|
||||
#endif
|
||||
|
||||
#if USE_NEW_PARSER
|
||||
Opm::EclipseWriter outputWriter(param, newParserDeck, share_obj(*grid->c_grid()));
|
||||
#else
|
||||
Opm::EclipseWriter outputWriter(param, deck, share_obj(*grid->c_grid()));
|
||||
#endif
|
||||
|
||||
// Rock and fluid init
|
||||
#if USE_NEW_PARSER
|
||||
props.reset(new BlackoilPropertiesFromDeck(newParserDeck, *grid->c_grid(), param));
|
||||
new_props.reset(new BlackoilPropsAdFromDeck(newParserDeck, *grid->c_grid()));
|
||||
#else
|
||||
props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param));
|
||||
new_props.reset(new BlackoilPropsAdFromDeck(*deck, *grid->c_grid()));
|
||||
#endif
|
||||
|
||||
// check_well_controls = param.getDefault("check_well_controls", false);
|
||||
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
|
||||
// Rock compressibility.
|
||||
#if USE_NEW_PARSER
|
||||
rock_comp.reset(new RockCompressibility(newParserDeck));
|
||||
#else
|
||||
rock_comp.reset(new RockCompressibility(*deck));
|
||||
#endif
|
||||
|
||||
// Gravity.
|
||||
#if USE_NEW_PARSER
|
||||
gravity[2] = newParserDeck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
|
||||
#else
|
||||
gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
|
||||
#endif
|
||||
|
||||
// Init state variables (saturation and pressure).
|
||||
if (param.has("init_saturation")) {
|
||||
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
|
||||
@ -159,11 +133,7 @@ try
|
||||
}
|
||||
}
|
||||
} else {
|
||||
#if USE_NEW_PARSER
|
||||
initBlackoilStateFromDeck(*grid->c_grid(), *props, newParserDeck, gravity[2], state);
|
||||
#else
|
||||
initBlackoilStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state);
|
||||
#endif
|
||||
}
|
||||
|
||||
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
|
||||
@ -195,7 +165,6 @@ try
|
||||
param.writeParam(output_dir + "/simulation.param");
|
||||
}
|
||||
|
||||
#if USE_NEW_PARSER
|
||||
std::cout << "\n\n================ Starting main simulation loop ===============\n"
|
||||
<< std::flush;
|
||||
|
||||
@ -231,7 +200,7 @@ try
|
||||
simtimer.setCurrentStepNum(reportStepIdx);
|
||||
|
||||
if (reportStepIdx == 0) {
|
||||
outputWriter.writeInit(simtimer, state, well_state.basicWellState());
|
||||
outputWriter.writeInit(simtimer);
|
||||
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
|
||||
}
|
||||
|
||||
@ -242,7 +211,6 @@ try
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells,
|
||||
linsolver,
|
||||
grav);
|
||||
SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
|
||||
|
||||
++simtimer;
|
||||
@ -260,85 +228,6 @@ try
|
||||
fullReport.reportParam(tot_os);
|
||||
warnIfUnusedParams(param);
|
||||
}
|
||||
#else
|
||||
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.
|
||||
WellStateFullyImplicitBlackoil 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);
|
||||
}
|
||||
|
||||
if (epoch == 0)
|
||||
outputWriter.writeInit(simtimer, state, well_state.basicWellState());
|
||||
|
||||
// Create and run simulator.
|
||||
SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
|
||||
*grid->c_grid(),
|
||||
*new_props,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells,
|
||||
linsolver,
|
||||
grav);
|
||||
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
|
||||
|
||||
if (epoch == 0) {
|
||||
warnIfUnusedParams(param);
|
||||
}
|
||||
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
|
||||
if (output) {
|
||||
epoch_rep.reportParam(outStream);
|
||||
}
|
||||
// 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);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
catch (const std::exception &e) {
|
||||
std::cerr << "Program threw an exception: " << e.what() << "\n";
|
||||
|
||||
@ -288,7 +288,10 @@ spdiag(const AutoDiffBlock<double>::V& d)
|
||||
public:
|
||||
typedef AutoDiffBlock<Scalar> ADB;
|
||||
|
||||
Selector(const typename ADB::V& selection_basis)
|
||||
enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero };
|
||||
|
||||
Selector(const typename ADB::V& selection_basis,
|
||||
CriterionForLeftElement crit = GreaterEqualZero)
|
||||
{
|
||||
// Define selector structure.
|
||||
const int n = selection_basis.size();
|
||||
@ -296,10 +299,33 @@ spdiag(const AutoDiffBlock<double>::V& d)
|
||||
left_elems_.reserve(n);
|
||||
right_elems_.reserve(n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
if (selection_basis[i] < 0.0) {
|
||||
right_elems_.push_back(i);
|
||||
} else {
|
||||
bool chooseleft = false;
|
||||
switch (crit) {
|
||||
case GreaterEqualZero:
|
||||
chooseleft = selection_basis[i] >= 0.0;
|
||||
break;
|
||||
case GreaterZero:
|
||||
chooseleft = selection_basis[i] > 0.0;
|
||||
break;
|
||||
case Zero:
|
||||
chooseleft = selection_basis[i] == 0.0;
|
||||
break;
|
||||
case NotEqualZero:
|
||||
chooseleft = selection_basis[i] != 0.0;
|
||||
break;
|
||||
case LessZero:
|
||||
chooseleft = selection_basis[i] < 0.0;
|
||||
break;
|
||||
case LessEqualZero:
|
||||
chooseleft = selection_basis[i] <= 0.0;
|
||||
break;
|
||||
default:
|
||||
OPM_THROW(std::logic_error, "No such criterion: " << crit);
|
||||
}
|
||||
if (chooseleft) {
|
||||
left_elems_.push_back(i);
|
||||
} else {
|
||||
right_elems_.push_back(i);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -560,7 +560,7 @@ namespace {
|
||||
|
||||
template<class T>
|
||||
void FullyImplicitBlackoilSolver<T>::computeWellConnectionPressures(const SolutionState& state,
|
||||
const WellStateFullyImplicitBlackoil& xw)
|
||||
const WellStateFullyImplicitBlackoil& xw)
|
||||
{
|
||||
using namespace Opm::AutoDiffGrid;
|
||||
// 1. Compute properties required by computeConnectionPressureDelta().
|
||||
@ -702,7 +702,8 @@ namespace {
|
||||
|
||||
|
||||
|
||||
template<class T>
|
||||
|
||||
template <class T>
|
||||
void FullyImplicitBlackoilSolver<T>::addWellEq(const SolutionState& state,
|
||||
WellStateFullyImplicitBlackoil& xw)
|
||||
{
|
||||
@ -806,14 +807,17 @@ namespace {
|
||||
|
||||
// check for dead wells
|
||||
V isNotDeadWells = V::Constant(nw,1.0);
|
||||
for (int c = 0; c < nw; ++c){
|
||||
if (wbqt.value()[c] == 0)
|
||||
isNotDeadWells[c] = 0;
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
if (wbqt.value()[w] == 0) {
|
||||
isNotDeadWells[w] = 0;
|
||||
}
|
||||
}
|
||||
// compute wellbore mixture at std conds
|
||||
Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
|
||||
std::vector<ADB> mix_s(np, ADB::null());
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
mix_s[phase] = isNotDeadWells * wbq[phase]/wbqt;
|
||||
const int pos = pu.phase_pos[phase];
|
||||
mix_s[phase] = notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
|
||||
}
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user