Merge branch 'master' into ert
Conflicts: Makefile.am configure.ac examples/Makefile.am opm/core/GridManager.cpp opm/core/eclipse/EclipseGridParser.cpp opm/core/grid/cpgpreprocess/preprocess.h tests/Makefile.am
This commit is contained in:
commit
00bbfa06a8
119 files changed
+7557
-2902
No files matched your search
+15
-13
@@ -1,15 +1,15 @@
|
||||
# Build-time flags needed to form example programs
|
||||
ERT_INCLUDE_PATH = $(ERT_ROOT)/include
|
||||
|
||||
AM_CPPFLAGS = \
|
||||
-I$(top_srcdir) \
|
||||
$(BOOST_CPPFLAGS) \
|
||||
AM_CPPFLAGS = \
|
||||
-I$(top_srcdir) \
|
||||
$(OPM_BOOST_CPPFLAGS) \
|
||||
-I$(ERT_INCLUDE_PATH)
|
||||
|
||||
# All targets link to the library
|
||||
LDADD = \
|
||||
$(top_builddir)/libopmcore.la \
|
||||
$(BOOST_FILESYSTEM_LIB) \
|
||||
LDADD = \
|
||||
$(top_builddir)/lib/libopmcore.la \
|
||||
$(BOOST_FILESYSTEM_LIB) \
|
||||
$(BOOST_SYSTEM_LIB)
|
||||
|
||||
# ----------------------------------------------------------------------
|
||||
@@ -17,11 +17,12 @@ $(BOOST_SYSTEM_LIB)
|
||||
#
|
||||
# Please keep the list sorted.
|
||||
|
||||
noinst_PROGRAMS = \
|
||||
refine_wells \
|
||||
scaneclipsedeck \
|
||||
sim_2p_incomp_reorder \
|
||||
sim_wateroil \
|
||||
noinst_PROGRAMS = \
|
||||
refine_wells \
|
||||
scaneclipsedeck \
|
||||
sim_2p_comp_reorder \
|
||||
sim_2p_incomp_reorder \
|
||||
sim_wateroil \
|
||||
wells_example
|
||||
|
||||
# ----------------------------------------------------------------------
|
||||
@@ -32,6 +33,7 @@ wells_example
|
||||
# Please maintain sort order from "noinst_PROGRAMS".
|
||||
|
||||
refine_wells_SOURCES = refine_wells.cpp
|
||||
sim_2p_comp_reorder_SOURCES = sim_2p_comp_reorder.cpp
|
||||
sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp
|
||||
sim_wateroil_SOURCES = sim_wateroil.cpp
|
||||
wells_example_SOURCES = wells_example.cpp
|
||||
@@ -43,7 +45,7 @@ if UMFPACK
|
||||
noinst_PROGRAMS += spu_2p
|
||||
|
||||
spu_2p_SOURCES = spu_2p.cpp
|
||||
spu_2p_LDADD = \
|
||||
$(LDADD) \
|
||||
spu_2p_LDADD = \
|
||||
$(LDADD) \
|
||||
$(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS)
|
||||
endif
|
||||
@@ -0,0 +1,285 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#if HAVE_CONFIG_H
|
||||
#include "config.h"
|
||||
#endif // HAVE_CONFIG_H
|
||||
|
||||
#include <opm/core/pressure/FlowBCManager.hpp>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/GridManager.hpp>
|
||||
#include <opm/core/newwells.h>
|
||||
#include <opm/core/wells/WellsManager.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/initState.hpp>
|
||||
#include <opm/core/simulator/SimulatorReport.hpp>
|
||||
#include <opm/core/simulator/SimulatorTimer.hpp>
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/miscUtilities.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
#include <opm/core/fluid/BlackoilPropertiesBasic.hpp>
|
||||
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
|
||||
#include <opm/core/fluid/RockCompressibility.hpp>
|
||||
|
||||
#include <opm/core/linalg/LinearSolverFactory.hpp>
|
||||
|
||||
#include <opm/core/simulator/BlackoilState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
#include <opm/core/simulator/SimulatorCompressibleTwophase.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)
|
||||
{
|
||||
using namespace Opm;
|
||||
|
||||
std::cout << "\n================ Test program for weakly compressible 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<BlackoilPropertiesInterface> props;
|
||||
boost::scoped_ptr<RockCompressibility> rock_comp;
|
||||
BlackoilState state;
|
||||
// bool check_well_controls = false;
|
||||
// int max_well_control_iterations = 0;
|
||||
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 BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param));
|
||||
// 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));
|
||||
// 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);
|
||||
}
|
||||
initBlackoilSurfvol(*grid->c_grid(), *props, state);
|
||||
} 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);
|
||||
initBlackoilSurfvol(*grid->c_grid(), *props, state);
|
||||
}
|
||||
|
||||
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
|
||||
int num_cells = grid->c_grid()->number_of_cells;
|
||||
std::vector<double> src(num_cells, 0.0);
|
||||
if (use_deck) {
|
||||
// Do nothing, wells will be the driving force, not source terms.
|
||||
} else {
|
||||
// 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;
|
||||
src[0] = flow_per_sec;
|
||||
src[num_cells - 1] = -flow_per_sec;
|
||||
}
|
||||
|
||||
// Boundary conditions.
|
||||
FlowBCManager bcs;
|
||||
if (param.getDefault("use_pside", false)) {
|
||||
int pside = param.get<int>("pside");
|
||||
double pside_pressure = param.get<double>("pside_pressure");
|
||||
bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
|
||||
}
|
||||
|
||||
// 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 (...) {
|
||||
THROW("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.
|
||||
WellsManager wells; // no wells.
|
||||
SimulatorCompressibleTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells,
|
||||
src,
|
||||
bcs.c_bcs(),
|
||||
linsolver,
|
||||
grav);
|
||||
SimulatorTimer simtimer;
|
||||
simtimer.init(param);
|
||||
warnIfUnusedParams(param);
|
||||
WellState 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;
|
||||
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) {
|
||||
THROW("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.
|
||||
SimulatorCompressibleTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells,
|
||||
src,
|
||||
bcs.c_bcs(),
|
||||
linsolver,
|
||||
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);
|
||||
}
|
||||
|
||||
}
|
||||
@@ -43,9 +43,10 @@
|
||||
|
||||
#include <opm/core/simulator/TwophaseState.hpp>
|
||||
#include <opm/core/simulator/WellState.hpp>
|
||||
#include <opm/core/simulator/SimulatorTwophase.hpp>
|
||||
#include <opm/core/simulator/SimulatorIncompTwophase.hpp>
|
||||
|
||||
#include <boost/scoped_ptr.hpp>
|
||||
#include <boost/filesystem.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
@@ -53,6 +54,18 @@
|
||||
#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 -----------------
|
||||
@@ -81,9 +94,7 @@ main(int argc, char** argv)
|
||||
// Grid init
|
||||
grid.reset(new GridManager(*deck));
|
||||
// Rock and fluid init
|
||||
const int* gc = grid->c_grid()->global_cell;
|
||||
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
|
||||
props.reset(new IncompPropertiesFromDeck(*deck, global_cell));
|
||||
props.reset(new IncompPropertiesFromDeck(*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.
|
||||
@@ -157,17 +168,28 @@ main(int argc, char** argv)
|
||||
// Linear solver.
|
||||
LinearSolverFactory linsolver(param);
|
||||
|
||||
// Warn if any parameters are unused.
|
||||
// if (param.anyUnused()) {
|
||||
// std::cout << "-------------------- Unused parameters: --------------------\n";
|
||||
// param.displayUsage();
|
||||
// std::cout << "----------------------------------------------------------------" << std::endl;
|
||||
// }
|
||||
|
||||
// Write parameters used for later reference.
|
||||
// if (output) {
|
||||
// param.writeParam(output_dir + "/spu_2p.param");
|
||||
// }
|
||||
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 (...) {
|
||||
THROW("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"
|
||||
@@ -177,17 +199,19 @@ main(int argc, char** argv)
|
||||
SimulatorReport rep;
|
||||
if (!use_deck) {
|
||||
// Simple simulation without a deck.
|
||||
SimulatorTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
0, // wells
|
||||
src,
|
||||
bcs.c_bcs(),
|
||||
linsolver,
|
||||
grav);
|
||||
WellsManager wells; // no wells.
|
||||
SimulatorIncompTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells,
|
||||
src,
|
||||
bcs.c_bcs(),
|
||||
linsolver,
|
||||
grav);
|
||||
SimulatorTimer simtimer;
|
||||
simtimer.init(param);
|
||||
warnIfUnusedParams(param);
|
||||
WellState well_state;
|
||||
well_state.init(0, state);
|
||||
rep = simulator.run(simtimer, state, well_state);
|
||||
@@ -221,7 +245,7 @@ main(int argc, char** argv)
|
||||
<< "\n (number of steps: "
|
||||
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
|
||||
|
||||
// Create new wells, well_satate
|
||||
// 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,
|
||||
@@ -231,17 +255,22 @@ main(int argc, char** argv)
|
||||
}
|
||||
|
||||
// Create and run simulator.
|
||||
SimulatorTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells.c_wells(),
|
||||
src,
|
||||
bcs.c_bcs(),
|
||||
linsolver,
|
||||
grav);
|
||||
SimulatorIncompTwophase simulator(param,
|
||||
*grid->c_grid(),
|
||||
*props,
|
||||
rock_comp->isActive() ? rock_comp.get() : 0,
|
||||
wells,
|
||||
src,
|
||||
bcs.c_bcs(),
|
||||
linsolver,
|
||||
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();
|
||||
@@ -250,4 +279,11 @@ main(int argc, char** argv)
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
}
|
||||
+19
-39
@@ -68,8 +68,6 @@
|
||||
#include <vector>
|
||||
#include <numeric>
|
||||
|
||||
#define TRANSPORT_SOLVER_FIXED 1
|
||||
|
||||
|
||||
template <class State>
|
||||
static void outputState(const UnstructuredGrid& grid,
|
||||
@@ -143,7 +141,6 @@ main(int argc, char** argv)
|
||||
std::cout << "--------------- Reading parameters ---------------" << std::endl;
|
||||
|
||||
// Reading various control parameters.
|
||||
const bool use_reorder = param.getDefault("use_reorder", true);
|
||||
const bool output = param.getDefault("output", true);
|
||||
std::string output_dir;
|
||||
int output_interval = 1;
|
||||
@@ -178,9 +175,7 @@ main(int argc, char** argv)
|
||||
// Grid init
|
||||
grid.reset(new Opm::GridManager(deck));
|
||||
// Rock and fluid init
|
||||
const int* gc = grid->c_grid()->global_cell;
|
||||
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
|
||||
props.reset(new Opm::BlackoilPropertiesFromDeck(deck, global_cell));
|
||||
props.reset(new BlackoilPropertiesFromDeck(deck, *grid->c_grid(), param));
|
||||
// Wells init.
|
||||
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
|
||||
check_well_controls = param.getDefault("check_well_controls", false);
|
||||
@@ -233,28 +228,8 @@ main(int argc, char** argv)
|
||||
}
|
||||
}
|
||||
bool use_segregation_split = false;
|
||||
bool use_column_solver = false;
|
||||
bool use_gauss_seidel_gravity = false;
|
||||
if (use_gravity && use_reorder) {
|
||||
if (use_gravity) {
|
||||
use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
|
||||
if (use_segregation_split) {
|
||||
use_column_solver = param.getDefault("use_column_solver", use_column_solver);
|
||||
if (use_column_solver) {
|
||||
use_gauss_seidel_gravity = param.getDefault("use_gauss_seidel_gravity", use_gauss_seidel_gravity);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check that rock compressibility is not used with solvers that do not handle it.
|
||||
// int nl_pressure_maxiter = 0;
|
||||
// double nl_pressure_tolerance = 0.0;
|
||||
if (rock_comp->isActive()) {
|
||||
THROW("No rock compressibility in comp. pressure solver yet.");
|
||||
if (!use_reorder) {
|
||||
THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
|
||||
}
|
||||
// nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10);
|
||||
// nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
|
||||
}
|
||||
|
||||
// Source-related variables init.
|
||||
@@ -266,11 +241,11 @@ main(int argc, char** argv)
|
||||
// Extra rock init.
|
||||
std::vector<double> porevol;
|
||||
if (rock_comp->isActive()) {
|
||||
THROW("CompressibleTpfa solver does not handle this.");
|
||||
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
|
||||
} else {
|
||||
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
|
||||
}
|
||||
std::vector<double> initial_porevol = porevol;
|
||||
double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
|
||||
|
||||
|
||||
@@ -297,21 +272,20 @@ main(int argc, char** argv)
|
||||
const double nl_press_change_tol = param.getDefault("nl_press_change_tol", 10.0);
|
||||
const int nl_press_maxiter = param.getDefault("nl_press_maxiter", 20);
|
||||
const double *grav = use_gravity ? &gravity[0] : 0;
|
||||
Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, linsolver,
|
||||
Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
|
||||
nl_press_res_tol, nl_press_change_tol, nl_press_maxiter,
|
||||
grav, wells->c_wells());
|
||||
// Reordering solver.
|
||||
#if TRANSPORT_SOLVER_FIXED
|
||||
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
|
||||
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
|
||||
Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
|
||||
if (use_gauss_seidel_gravity) {
|
||||
if (use_segregation_split) {
|
||||
reorder_model.initGravity(grav);
|
||||
}
|
||||
#endif // TRANSPORT_SOLVER_FIXED
|
||||
|
||||
// Column-based gravity segregation solver.
|
||||
std::vector<std::vector<int> > columns;
|
||||
if (use_column_solver) {
|
||||
if (use_segregation_split) {
|
||||
Opm::extractColumn(*grid->c_grid(), columns);
|
||||
}
|
||||
|
||||
@@ -414,24 +388,30 @@ main(int argc, char** argv)
|
||||
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
|
||||
wells->c_wells(), well_state.perfRates(), reorder_src);
|
||||
|
||||
// Compute new porevolumes after pressure solve, if necessary.
|
||||
if (rock_comp->isActive()) {
|
||||
initial_porevol = porevol;
|
||||
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
|
||||
}
|
||||
// Solve transport.
|
||||
transport_timer.start();
|
||||
#if TRANSPORT_SOLVER_FIXED
|
||||
double stepsize = simtimer.currentStepLength();
|
||||
if (num_transport_substeps != 1) {
|
||||
stepsize /= double(num_transport_substeps);
|
||||
std::cout << "Making " << num_transport_substeps << " transport substeps." << std::endl;
|
||||
}
|
||||
for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
|
||||
reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0],
|
||||
&porevol[0], &reorder_src[0], stepsize, state.saturation());
|
||||
// Note that for now we do not handle rock compressibility,
|
||||
// although the transport solver should be able to.
|
||||
reorder_model.solve(&state.faceflux()[0], &state.pressure()[0],
|
||||
&porevol[0], &initial_porevol[0], &reorder_src[0], stepsize,
|
||||
state.saturation(), state.surfacevol());
|
||||
// Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
|
||||
if (use_segregation_split) {
|
||||
THROW("Segregation not implemented yet.");
|
||||
// reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation());
|
||||
reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0],
|
||||
stepsize, state.saturation(), state.surfacevol());
|
||||
}
|
||||
}
|
||||
#endif // TRANSPORT_SOLVER_FIXED
|
||||
transport_timer.stop();
|
||||
double tt = transport_timer.secsSinceStart();
|
||||
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
|
||||
|
||||
+2
-4
@@ -317,9 +317,7 @@ main(int argc, char** argv)
|
||||
// Grid init
|
||||
grid.reset(new Opm::GridManager(deck));
|
||||
// Rock and fluid init
|
||||
const int* gc = grid->c_grid()->global_cell;
|
||||
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
|
||||
props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
|
||||
props.reset(new Opm::IncompPropertiesFromDeck(deck, *grid->c_grid()));
|
||||
// Wells init.
|
||||
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
|
||||
check_well_controls = param.getDefault("check_well_controls", false);
|
||||
@@ -503,7 +501,7 @@ main(int argc, char** argv)
|
||||
|
||||
// Write parameters used for later reference.
|
||||
if (output) {
|
||||
param.writeParam(output_dir + "/spu_2p.param");
|
||||
param.writeParam(output_dir + "/simulation.param");
|
||||
}
|
||||
|
||||
// Main simulation loop.
|
||||
|
||||
@@ -38,10 +38,8 @@ int main(int argc, char** argv)
|
||||
// Finally handle the wells
|
||||
WellsManager wells(parser, *grid.c_grid(), NULL);
|
||||
|
||||
std::vector<int> global_cells(grid.c_grid()->global_cell, grid.c_grid()->global_cell + grid.c_grid()->number_of_cells);
|
||||
|
||||
double gravity[3] = {0.0, 0.0, parameters.getDefault<double>("gravity", 0.0)};
|
||||
IncompPropertiesFromDeck incomp_properties(parser, global_cells);
|
||||
IncompPropertiesFromDeck incomp_properties(parser, *grid.c_grid());
|
||||
|
||||
RockCompressibility rock_comp(parser);
|
||||
|
||||
|
||||
Reference in new issue
Block a user