OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'atgeirr/master'

Browse Source
This commit is contained in:
Xavier Raynaud 2012-08-23 15:25:19 +02:00
commit 73fc6afab9
13 changed files with 508 additions and 159 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/Makefile.am
View File

@ -17,6 +17,7 @@ $(BOOST_SYSTEM_LIB)
noinst_PROGRAMS = \
refine_wells \
scaneclipsedeck \
sim_2p_comp_reorder \
sim_2p_incomp_reorder \
sim_wateroil \
wells_example
@ -29,6 +30,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

285
examples/sim_2p_comp_reorder.cpp Normal file
View File

@ -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()));
// 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);
}
}

4
examples/sim_wateroil.cpp
View File

@ -280,7 +280,7 @@ main(int argc, char** argv)
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
if (use_segregation_split) {
reorder_model.initGravity();
reorder_model.initGravity(grav);
}
// Column-based gravity segregation solver.
@ -409,7 +409,7 @@ main(int argc, char** argv)
// Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
if (use_segregation_split) {
reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0],
stepsize, grav, state.saturation());
stepsize, state.saturation(), state.surfacevol());
}
}
transport_timer.stop();

36
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -77,7 +77,8 @@ namespace Opm
wells_(wells),
htrans_(grid.cell_facepos[ grid.number_of_cells ]),
trans_ (grid.number_of_faces),
allcells_(grid.number_of_cells)
allcells_(grid.number_of_cells),
singular_(false)
{
if (wells_ && (wells_->number_of_phases != props.numPhases())) {
THROW("Inconsistent number of phases specified (wells vs. props): "
@ -189,6 +190,21 @@ namespace Opm
/// @brief After solve(), was the resulting pressure singular.
/// Returns true if the pressure is singular in the following
/// sense: if everything is incompressible and there are no
/// pressure conditions, the absolute values of the pressure
/// solution are arbitrary. (But the differences in pressure
/// are significant.)
bool CompressibleTpfa::singularPressure() const
{
return singular_;
}
/// Compute well potentials.
void CompressibleTpfa::computeWellPotentials(const BlackoilState& state)
{
@ -487,16 +503,20 @@ namespace Opm
cq.Af = &face_A_[0];
cq.phasemobf = &face_phasemob_[0];
cq.voldiscr = &cell_voldisc_[0];
int was_adjusted = 0;
if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) {
cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
&face_gravcap_[0], cell_press, well_bhp,
&porevol_[0], h_);
was_adjusted =
cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
&face_gravcap_[0], cell_press, well_bhp,
&porevol_[0], h_);
} else {
cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0],
&face_gravcap_[0], cell_press, well_bhp,
&porevol_[0], &initial_porevol_[0],
&rock_comp_[0], h_);
was_adjusted =
cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0],
&face_gravcap_[0], cell_press, well_bhp,
&porevol_[0], &initial_porevol_[0],
&rock_comp_[0], h_);
}
singular_ = (was_adjusted == 1);
}

18
opm/core/pressure/CompressibleTpfa.hpp
View File

@ -81,6 +81,14 @@ namespace Opm
BlackoilState& state,
WellState& well_state);
/// @brief After solve(), was the resulting pressure singular.
/// Returns true if the pressure is singular in the following
/// sense: if everything is incompressible and there are no
/// pressure conditions, the absolute values of the pressure
/// solution are arbitrary. (But the differences in pressure
/// are significant.)
bool singularPressure() const;
private:
void computePerSolveDynamicData(const double dt,
const BlackoilState& state,
@ -143,11 +151,11 @@ namespace Opm
std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
// The update to be applied to the pressures (cell and bhp).
std::vector<double> pressure_increment_;
// True if the matrix assembled would be singular but for the
// adjustment made in the cfs_*_assemble() calls. This happens
// if everything is incompressible and there are no pressure
// conditions.
bool singular_;
};
} // namespace Opm

42
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@ -1156,7 +1156,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
/* ---------------------------------------------------------------------- */
void
int
cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
double dt ,
struct cfs_tpfa_res_forces *forces ,
@ -1170,7 +1170,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_data *h )
/* ---------------------------------------------------------------------- */
{
int res_is_neumann, well_is_neumann, c, np2;
int res_is_neumann, well_is_neumann, c, np2, singular;
csrmatrix_zero( h->J);
vector_zero (h->J->m, h->F);
@ -1207,14 +1207,17 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
assemble_sources(dt, forces->src, h);
}
if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) {
h->J->sa[0] *= 2;
singular = res_is_neumann && well_is_neumann && h->pimpl->is_incomp;
if (singular) {
h->J->sa[0] *= 2.0;
}
return singular;
}
/* ---------------------------------------------------------------------- */
void
int
cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G ,
double dt ,
@ -1240,29 +1243,18 @@ cfs_tpfa_res_comprock_assemble(
* porevol(pressure)*rock_comp(pressure)/dt
*/
int c, w, well_is_neumann, rock_is_incomp;
int c, rock_is_incomp, singular;
size_t j;
double dpv;
const struct Wells* W;
double dpv;
/* Assemble usual system (without rock compressibility). */
cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
cpress, wpress, porevol0, h);
/* Check if we have only Neumann wells. */
well_is_neumann = 1;
W = forces->wells->W;
for (w = 0; well_is_neumann && w < W->number_of_wells; w++) {
if ((W->ctrls[w]->current >= 0) && /* OPEN? */
(W->ctrls[w]->type[ W->ctrls[w]->current ] == BHP)) {
well_is_neumann = 0;
}
}
singular = cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
cpress, wpress, porevol0, h);
/* If we made a singularity-removing adjustment in the
regular assembly, we undo it here. */
if (well_is_neumann && h->pimpl->is_incomp) {
h->J->sa[0] /= 2;
if (singular) {
h->J->sa[0] /= 2.0;
}
/* Add new terms to residual and Jacobian. */
@ -1280,9 +1272,11 @@ cfs_tpfa_res_comprock_assemble(
}
/* Re-do the singularity-removing adjustment if necessary */
if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) {
h->J->sa[0] *= 2;
if (rock_is_incomp && singular) {
h->J->sa[0] *= 2.0;
}
return rock_is_incomp && singular;
}

13
opm/core/pressure/tpfa/cfs_tpfa_residual.h
View File

@ -59,7 +59,11 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
void
cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
void
/* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible and
there are no pressure conditions on wells or boundaries.
Otherwise return 0. */
int
cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
double dt,
struct cfs_tpfa_res_forces *forces,
@ -72,7 +76,12 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
const double *porevol,
struct cfs_tpfa_res_data *h);
void
/* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible, the
rock is incompressible, and there are no pressure conditions on
wells or boundaries.
Otherwise return 0. */
int
cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G,
double dt,

50
opm/core/simulator/SimulatorCompressibleTwophase.cpp
View File

@ -17,7 +17,6 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
@ -153,7 +152,7 @@ namespace Opm
catch (...) {
THROW("Creating directories failed: " << fpath);
}
vtkfilename << "/" << std::setw(3) << std::setfill('0') << step << ".vtu";
vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
@ -228,36 +227,6 @@ namespace Opm
}
static bool allNeumannBCs(const FlowBoundaryConditions* bcs)
{
if (bcs == NULL) {
return true;
} else {
return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE)
== bcs->type + bcs->nbc;
}
}
static bool allRateWells(const Wells* wells)
{
if (wells == NULL) {
return true;
}
const int nw = wells->number_of_wells;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells->ctrls[w];
if (wc->current >= 0) {
if (wc->type[wc->current] == BHP) {
return false;
}
}
}
return true;
}
// \TODO: make CompressibleTpfa take src and bcs.
SimulatorCompressibleTwophase::Impl::Impl(const parameter::ParameterGroup& param,
@ -276,7 +245,6 @@ namespace Opm
wells_(wells_manager.c_wells()),
src_(src),
bcs_(bcs),
gravity_(gravity),
psolver_(grid, props, rock_comp, linsolver,
param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
@ -310,7 +278,7 @@ namespace Opm
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false);
if (gravity != 0 && use_segregation_split_){
tsolver_.initGravity();
tsolver_.initGravity(gravity);
extractColumn(grid_, columns_);
}
@ -402,12 +370,12 @@ namespace Opm
std::vector<double> initial_pressure = state.pressure();
psolver_.solve(timer.currentStepLength(), state, well_state);
// Renormalize pressure if rock is incompressible, and
// there are no pressure conditions (bcs or wells).
// It is deemed sufficient for now to renormalize
// using geometric volume instead of pore volume.
if ((rock_comp_ == NULL || !rock_comp_->isActive())
&& allNeumannBCs(bcs_) && allRateWells(wells_)) {
// Renormalize pressure if both fluids and rock are
// incompressible, and there are no pressure
// conditions (bcs or wells). It is deemed sufficient
// for now to renormalize using geometric volume
// instead of pore volume.
if (psolver_.singularPressure()) {
// Compute average pressures of previous and last
// step, and total volume.
double av_prev_press = 0.0;
@ -484,7 +452,7 @@ namespace Opm
transport_src, stepsize, injected, produced);
if (gravity_ != 0 && use_segregation_split_) {
tsolver_.solveGravity(columns_, &state.pressure()[0], &initial_porevol[0],
stepsize, gravity_, state.saturation());
stepsize, state.saturation(), state.surfacevol());
}
}
transport_timer.stop();

7
opm/core/simulator/SimulatorIncompTwophase.cpp
View File

@ -335,7 +335,7 @@ namespace Opm
computePorevolume(grid_, props_.porosity(), porevol);
}
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
std::vector<double> initial_porevol = porevol;
// Main simulation loop.
Opm::time::StopWatch pressure_timer;
@ -452,6 +452,7 @@ namespace Opm
// Update pore volumes if rock is compressible.
if (rock_comp_ && rock_comp_->isActive()) {
initial_porevol = porevol;
computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol);
}
@ -467,11 +468,11 @@ namespace Opm
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
}
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0],
tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0],
stepsize, state.saturation());
Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced);
if (use_segregation_split_) {
tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation());
tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation());
}
}
transport_timer.stop();

38
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
View File

@ -24,6 +24,7 @@
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/utility/RootFinders.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <fstream>
@ -38,7 +39,8 @@ namespace Opm
typedef RegulaFalsi<WarnAndContinueOnError> RootFinder;
TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid,
TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(
const UnstructuredGrid& grid,
const Opm::BlackoilPropertiesInterface& props,
const double tol,
const int maxit)
@ -51,6 +53,7 @@ namespace Opm
dt_(0.0),
saturation_(grid.number_of_cells, -1.0),
fractionalflow_(grid.number_of_cells, -1.0),
gravity_(0),
mob_(2*grid.number_of_cells, -1.0),
ia_upw_(grid.number_of_cells + 1, -1),
ja_upw_(grid.number_of_faces, -1),
@ -93,6 +96,11 @@ namespace Opm
props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL);
// Check non-miscibility requirement (only done for first cell).
if (A_[1] != 0.0 || A_[2] != 0.0) {
THROW("TransportModelCompressibleTwophase requires a property object without miscibility.");
}
std::vector<int> seq(grid_.number_of_cells);
std::vector<int> comp(grid_.number_of_cells + 1);
int ncomp;
@ -107,15 +115,9 @@ namespace Opm
&ia_downw_[0], &ja_downw_[0]);
reorderAndTransport(grid_, darcyflux);
toBothSat(saturation_, saturation);
// Compute surface volume as a postprocessing step from saturation and A_
surfacevol = saturation;
const int np = props_.numPhases();
for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
for (int phase = 0; phase < np; ++phase) {
surfacevol[np*cell + phase] *= A_[np*np*cell + np*phase + phase];
}
}
computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]);
}
// Residual function r(s) for a single-cell implicit Euler transport
@ -384,7 +386,7 @@ namespace Opm
void TransportModelCompressibleTwophase::initGravity()
void TransportModelCompressibleTwophase::initGravity(const double* grav)
{
// Set up transmissibilities.
std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
@ -393,12 +395,15 @@ namespace Opm
gravflux_.resize(nf);
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &trans_[0]);
// Remember gravity vector.
gravity_ = grav;
}
void TransportModelCompressibleTwophase::initGravityDynamic(const double* grav)
void TransportModelCompressibleTwophase::initGravityDynamic()
{
// Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f)
// + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ]
@ -420,7 +425,7 @@ namespace Opm
for (int ci = 0; ci < 2; ++ci) {
double gdz = 0.0;
for (int d = 0; d < dim; ++d) {
gdz += grav[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]);
gdz += gravity_[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]);
}
gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]);
}
@ -504,11 +509,11 @@ namespace Opm
const double* pressure,
const double* porevolume0,
const double dt,
const double* grav,
std::vector<double>& saturation)
std::vector<double>& saturation,
std::vector<double>& surfacevol)
{
// Assume that solve() has already been called, so that A_ is current.
initGravityDynamic(grav);
initGravityDynamic();
// Initialize mobilities.
const int nc = grid_.number_of_cells;
@ -541,6 +546,9 @@ namespace Opm
std::cout << "Gauss-Seidel column solver average iterations: "
<< double(num_iters)/double(columns.size()) << std::endl;
toBothSat(saturation_, saturation);
// Compute surface volume as a postprocessing step from saturation and A_
computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]);
}
} // namespace Opm

15
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
View File

@ -30,6 +30,8 @@ namespace Opm
class BlackoilPropertiesInterface;
/// Implements a reordering transport solver for compressible,
/// non-miscible two-phase flow.
class TransportModelCompressibleTwophase : public TransportModelInterface
{
public:
@ -52,6 +54,7 @@ namespace Opm
/// \param[in] source Transport source term.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
/// \param[in, out] surfacevol Surface volume densities for each phase.
void solve(const double* darcyflux,
const double* pressure,
const double* porevolume0,
@ -62,7 +65,8 @@ namespace Opm
std::vector<double>& surfacevol);
/// Initialise quantities needed by gravity solver.
void initGravity();
/// \param[in] grav Gravity vector
void initGravity(const double* grav);
/// Solve for gravity segregation.
/// This uses a column-wise nonlinear Gauss-Seidel approach.
@ -72,14 +76,14 @@ namespace Opm
/// \param[in] columns Vector of cell-columns.
/// \param[in] porevolume0 Array of pore volumes at start of timestep.
/// \param[in] dt Time step.
/// \param[in] grav Gravity vector.
/// \param[in, out] saturation Phase saturations.
/// \param[in, out] surfacevol Surface volume densities for each phase.
void solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure,
const double* porevolume0,
const double dt,
const double* grav,
std::vector<double>& saturation);
std::vector<double>& saturation,
std::vector<double>& surfacevol);
private:
virtual void solveSingleCell(const int cell);
@ -88,7 +92,7 @@ namespace Opm
const int pos,
const double* gravflux);
int solveGravityColumn(const std::vector<int>& cells);
void initGravityDynamic(const double* grav);
void initGravityDynamic();
private:
const UnstructuredGrid& grid_;
@ -110,6 +114,7 @@ namespace Opm
std::vector<double> saturation_; // P (= num. phases) per cell
std::vector<double> fractionalflow_; // = m[0]/(m[0] + m[1]) per cell
// For gravity segregation.
const double* gravity_;
std::vector<double> trans_;
std::vector<double> density_;
std::vector<double> gravflux_;

112
opm/core/utility/miscUtilitiesBlackoil.cpp
View File

@ -48,39 +48,39 @@ namespace Opm
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
const std::vector<double>& press,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& src,
const double dt,
double* injected,
double* produced)
const std::vector<double>& s,
const std::vector<double>& src,
const double dt,
double* injected,
double* produced)
{
const int num_cells = src.size();
const int np = s.size()/src.size();
if (int(s.size()) != num_cells*np) {
THROW("Sizes of s and src vectors do not match.");
}
std::fill(injected, injected + np, 0.0);
std::fill(produced, produced + np, 0.0);
const int num_cells = src.size();
const int np = s.size()/src.size();
if (int(s.size()) != num_cells*np) {
THROW("Sizes of s and src vectors do not match.");
}
std::fill(injected, injected + np, 0.0);
std::fill(produced, produced + np, 0.0);
std::vector<double> visc(np);
std::vector<double> mob(np);
for (int c = 0; c < num_cells; ++c) {
if (src[c] > 0.0) {
injected[0] += src[c]*dt;
} else if (src[c] < 0.0) {
const double flux = -src[c]*dt;
const double* sat = &s[np*c];
props.relperm(1, sat, &c, &mob[0], 0);
std::vector<double> mob(np);
for (int c = 0; c < num_cells; ++c) {
if (src[c] > 0.0) {
injected[0] += src[c]*dt;
} else if (src[c] < 0.0) {
const double flux = -src[c]*dt;
const double* sat = &s[np*c];
props.relperm(1, sat, &c, &mob[0], 0);
props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
double totmob = 0.0;
for (int p = 0; p < np; ++p) {
mob[p] /= visc[p];
totmob += mob[p];
}
for (int p = 0; p < np; ++p) {
produced[p] += (mob[p]/totmob)*flux;
}
}
}
double totmob = 0.0;
for (int p = 0; p < np; ++p) {
mob[p] /= visc[p];
totmob += mob[p];
}
for (int p = 0; p < np; ++p) {
produced[p] += (mob[p]/totmob)*flux;
}
}
}
}
@ -93,11 +93,11 @@ namespace Opm
/// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobilities.
void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<int>& cells,
const std::vector<double>& press,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& totmob)
const std::vector<double>& s,
std::vector<double>& totmob)
{
std::vector<double> pmobc;
@ -126,12 +126,12 @@ namespace Opm
/// @param[out] totmob total mobility
/// @param[out] omega fractional-flow weighted fluid densities.
void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega)
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega)
{
std::vector<double> pmobc;
@ -185,10 +185,10 @@ namespace Opm
props.relperm(nc, &s[0], &cells[0],
&pmobc[0], dpmobc);
std::transform(pmobc.begin(), pmobc.end(),
mu.begin(),
pmobc.begin(),
std::divides<double>());
std::transform(pmobc.begin(), pmobc.end(),
mu.begin(),
pmobc.begin(),
std::divides<double>());
}
/// Computes the fractional flow for each cell in the cells argument
@ -220,5 +220,35 @@ namespace Opm
}
}
/// Computes the surface volume densities from saturations by the formula
/// z = A s
/// for a number of data points, where z is the surface volume density,
/// s is the saturation (both as column vectors) and A is the
/// phase-to-component relation matrix.
/// @param[in] n number of data points
/// @param[in] np number of phases, must be 2 or 3
/// @param[in] A array containing n square matrices of size num_phases^2,
/// in Fortran ordering, typically the output of a call
/// to the matrix() method of a BlackoilProperties* class.
/// @param[in] saturation concatenated saturation values (for all P phases)
/// @param[out] surfacevol concatenated surface-volume values (for all P phases)
void computeSurfacevol(const int n,
const int np,
const double* A,
const double* saturation,
double* surfacevol)
{
// Note: since this is a simple matrix-vector product, it can
// be done by a BLAS call, but then we have to reorder the A
// matrix data.
std::fill(surfacevol, surfacevol + n*np, 0.0);
for (int i = 0; i < n; ++i) {
for (int row = 0; row < np; ++row) {
for (int col = 0; col < np; ++col) {
surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col];
}
}
}
}
} // namespace Opm

45
opm/core/utility/miscUtilitiesBlackoil.hpp
View File

@ -46,11 +46,11 @@ namespace Opm
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& src,
const double dt,
double* injected,
double* produced);
const std::vector<double>& s,
const std::vector<double>& src,
const double dt,
double* injected,
double* produced);
/// @brief Computes total mobility for a set of saturation values.
/// @param[in] props rock and fluid properties
@ -60,11 +60,11 @@ namespace Opm
/// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobilities.
void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& totmob);
const std::vector<double>& s,
std::vector<double>& totmob);
/// @brief Computes total mobility and omega for a set of saturation values.
/// @param[in] props rock and fluid properties
@ -75,12 +75,12 @@ namespace Opm
/// @param[out] totmob total mobility
/// @param[out] omega fractional-flow weighted fluid densities.
void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props,
const std::vector<int>& cells,
const std::vector<int>& cells,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega);
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega);
/// @brief Computes phase mobilities for a set of saturation values.
@ -96,7 +96,7 @@ namespace Opm
const std::vector<double>& z,
const std::vector<double>& s,
std::vector<double>& pmobc);
/// Computes the fractional flow for each cell in the cells argument
/// @param[in] props rock and fluid properties
@ -112,6 +112,25 @@ namespace Opm
const std::vector<double>& s,
std::vector<double>& fractional_flows);
/// Computes the surface volume densities from saturations by the formula
/// z = A s
/// for a number of data points, where z is the surface volume density,
/// s is the saturation (both as column vectors) and A is the
/// phase-to-component relation matrix.
/// @param[in] n number of data points
/// @param[in] np number of phases, must be 2 or 3
/// @param[in] A array containing n square matrices of size num_phases,
/// in Fortran ordering, typically the output of a call
/// to the matrix() method of a BlackoilProperties* class.
/// @param[in] saturation concatenated saturation values (for all P phases)
/// @param[out] surfacevol concatenated surface-volume values (for all P phases)
void computeSurfacevol(const int n,
const int np,
const double* A,
const double* saturation,
double* surfacevol);
} // namespace Opm
#endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED