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

Cleaning the initialization code

Browse Source 
-remove whitespaces
-fix documentation
This commit is contained in:
Tor Harald Sandve 2017-11-21 12:27:09 +01:00
parent e613f2b8c5
commit 3486e98ae3
5 changed files with 76 additions and 128 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
examples/compute_initial_state.cpp
View File

@ -1,5 +1,6 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
@ -27,20 +28,17 @@
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <boost/filesystem.hpp>
#include <fstream>
@ -105,17 +103,8 @@ namespace
}
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
@ -134,18 +123,19 @@ try
const double grav = param.getDefault("gravity", unit::gravity);
GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *gm.c_grid();
BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param);
warnIfUnusedParams(param);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
typedef FluidSystems::BlackOil<double> FluidSystem;
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/Opm::BlackoilPhases::Aqua,
/*nonWettingPhaseIdx=*/Opm::BlackoilPhases::Liquid,
/*gasPhaseIdx=*/Opm::BlackoilPhases::Vapour> MaterialTraits;
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
MaterialLawManager materialLawManager = MaterialLawManager();
@ -154,7 +144,6 @@ try
// Initialisation.
//initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
typedef FluidSystems::BlackOil<double> FluidSystem;
FluidSystem::initFromDeck(deck, eclipseState);
PhaseUsage pu = phaseUsageFromDeck(deck);
@ -169,17 +158,8 @@ try
state.pressure() = isc.press()[ref_phase];
convertSats<FluidSystem>(state.saturation(), isc.saturation(), pu);
if (state.hasCellData(std::string("GASOILRATIO"))) {
std::vector<double>& rs = state.getCellData(std::string("GASOILRATIO"));
rs = isc.rs();
}
if (state.hasCellData(std::string("RV"))){
std::vector<double>& rv = state.getCellData(std::string("RV"));
rv = isc.rv();
}
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
// Output.
const std::string output_dir = param.getDefault<std::string>("output_dir", "output");
@ -192,5 +172,3 @@ catch (const std::exception& e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

51
opm/core/simulator/EquilibrationHelpers.hpp
View File

@ -41,36 +41,33 @@ namespace Opm
{
namespace EQUIL {
template <class Props>
class DensityCalculator;
template <>
class DensityCalculator< BlackoilPropertiesInterface >;
namespace Miscibility {
class RsFunction;
class NoMixing;
template <class FluidSystem>
class RsVD;
template <class FluidSystem>
class RsSatAtContact;
}
template <class DensCalc>
class EquilReg;
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEq;
inline double satFromPc(const BlackoilPropertiesInterface& props,
template <class FluidSystem, class MaterialLaw, class MaterialLawManager >
inline double satFromPc(const MaterialLawManager& materialLawManager,
const int phase,
const int cell,
const double target_pc,
const bool increasing = false);
struct PcEqSum
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
const bool increasing = false)
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase1,
const int phase2,
const int cell,
const double target_pc);
const double target_pc)
} // namespace Equil
} // namespace Opm
@ -183,8 +180,7 @@ namespace Opm
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] rs Dissolved gas-oil ratio at @c depth.
*/
@ -195,7 +191,6 @@ namespace Opm
, depth_(depth)
, rs_(rs)
{
}
/**
@ -249,8 +244,7 @@ namespace Opm
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] rv Dissolved gas-oil ratio at @c depth.
*/
@ -261,7 +255,6 @@ namespace Opm
, depth_(depth)
, rv_(rv)
{
}
/**
@ -324,8 +317,7 @@ namespace Opm
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
@ -394,8 +386,7 @@ namespace Opm
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
@ -470,10 +461,9 @@ namespace Opm
* Constructor.
*
* \param[in] rec Equilibration data of current region.
* \param[in] density Density calculator of current region.
* \param[in] rs Calculator of dissolved gas-oil ratio.
* \param[in] rv Calculator of vapourised oil-gas ratio.
* \param[in] pu Summary of current active phases.
* \param[in] pvtRegionIdx The pvt region index
*/
EquilReg(const EquilRecord& rec,
std::shared_ptr<Miscibility::RsFunction> rs,
@ -483,7 +473,6 @@ namespace Opm
, rs_ (rs)
, rv_ (rv)
, pvtIdx_ (pvtIdx)
{
}
@ -547,7 +536,7 @@ namespace Opm
evaporationCalculator() const { return *this->rv_; }
/**
* Retrieve active fluid phase summary.
* Retrieve pvtIdx of the region.
*/
const int
pvtIdx() const { return this->pvtIdx_; }
@ -581,10 +570,7 @@ namespace Opm
fluidState_.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
fluidState_.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
std::fill(pc_, pc_ + FluidSystem::numPhases, 0.0);
}
double operator()(double s) const
{
const auto& matParams = materialLawManager_.materialLawParams(cell_);
@ -595,7 +581,6 @@ namespace Opm
return pc - target_pc_;
}
private:
const MaterialLawManager& materialLawManager_;
const int phase_;
const int cell_;
@ -656,11 +641,9 @@ namespace Opm
default: OPM_THROW(std::runtime_error, "Unknown phaseIdx .");
}
return -1.0;
}
/// Compute saturation of some phase corresponding to a given
/// capillary pressure.
template <class FluidSystem, class MaterialLaw, class MaterialLawManager >
@ -727,7 +710,6 @@ namespace Opm
double pc1 = pc_[FluidSystem::oilPhaseIdx] + sign1 * pc_[phase1_];
double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pc2 = pc_[FluidSystem::oilPhaseIdx] + sign2 * pc_[phase2_];
return pc1 + pc2 - target_pc_;
}
private:
@ -746,7 +728,6 @@ namespace Opm
/// Compute saturation of some phase corresponding to a given
/// capillary pressure, where the capillary pressure function
/// is given as a sum of two other functions.
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase1,
@ -809,8 +790,6 @@ namespace Opm
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
}
} // namespace Equil
} // namespace Opm

27
opm/core/simulator/initStateEquil.hpp
View File

@ -2,6 +2,7 @@
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
@ -39,12 +40,10 @@
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <array>
#include <cassert>
#include <utility>
@ -70,7 +69,6 @@ namespace Opm
*/
namespace EQUIL {
/**
* Compute initial phase pressures by means of equilibration.
*
@ -121,6 +119,11 @@ namespace Opm
/**
* Compute initial phase saturations by means of equilibration.
*
* \tparam FluidSystem The FluidSystem from opm-material
* Must be initialized before used.
*
* \tparam Grid Type of the grid
*
* \tparam Region Type of an equilibration region information
* base. Typically an instance of the EquilReg
* class template.
@ -132,10 +135,15 @@ namespace Opm
* as well as provide an inner type,
* const_iterator, to traverse the range.
*
* \tparam MaterialLawManager The MaterialLawManager from opm-material
*
* \param[in] G Grid.
* \param[in] reg Current equilibration region.
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] props Property object, needed for capillary functions.
* \param[in] materialLawManager The MaterialLawManager from opm-material
* \param[in] swat_init A vector of initial water saturations.
* The capillary pressure is scaled to fit these values
* \param[in] phase_pressures Phase pressures, one vector for each active phase,
* of pressure values in each cell in the current
* equilibration region.
@ -201,7 +209,6 @@ namespace Opm
const Grid& G )
{
std::vector<int> eqlnum;
if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) {
const int nc = UgGridHelpers::numCells(G);
eqlnum.resize(nc);
@ -239,7 +246,6 @@ namespace Opm
rs_(UgGridHelpers::numCells(G)),
rv_(UgGridHelpers::numCells(G))
{
//Check for presence of kw SWATINIT
if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatInit) {
const std::vector<double>& swat_init_ecl = eclipseState.
@ -252,13 +258,11 @@ namespace Opm
swat_init_[c] = swat_init_ecl[deck_pos];
}
}
// Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(eclipseState);
const auto& tables = eclipseState.getTableManager();
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(eclipseState, G));
setRegionPvtIdx(G, eclipseState, eqlmap);
// Create Rs functions.
@ -272,7 +276,6 @@ namespace Opm
continue;
}
const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].liveOilInitConstantRs()) {
if (rsvdTables.size() <= 0 ) {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available.");
@ -355,14 +358,10 @@ namespace Opm
const Vec& rv() const { return rv_; }
private:
typedef EquilReg EqReg;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
std::vector<int> regionPvtIdx_;
PVec pp_;
PVec sat_;
Vec rs_;
@ -410,10 +409,8 @@ namespace Opm
copyFromRegion(pressures[p], cells, pp_[p]);
copyFromRegion(sat[p], cells, sat_[p]);
}
const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
if (oil && gas) {
const int oilpos = FluidSystem::oilPhaseIdx;
const int gaspos = FluidSystem::gasPhaseIdx;

80
opm/core/simulator/initStateEquil_impl.hpp
View File

@ -2,6 +2,7 @@
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
@ -35,9 +36,6 @@
namespace Opm
{
namespace Details {
template <class RHS>
class RK4IVP {
public:
@ -108,8 +106,8 @@ namespace Opm
stepsize() const { return (span_[1] - span_[0]) / N_; }
};
namespace PhasePressODE {
template <class FluidSystem>
namespace PhasePressODE {
template <class FluidSystem>
class Water {
public:
Water(const double temp,
@ -317,7 +315,6 @@ namespace Opm
class Grid,
class Region,
class CellRange>
void
oil(const Grid& G ,
const Region& reg ,
@ -332,7 +329,6 @@ namespace Opm
typedef Oil<FluidSystem, typename Region::CalcDissolution> ODE;
const double T = 273.15 + 20; // standard temperature for now
ODE drho(T, reg.dissolutionCalculator(),
reg.pvtIdx(), grav);
@ -690,7 +686,6 @@ namespace Opm
sg = 1.0 - sw;
phase_saturations[waterpos][local_index] = sw;
phase_saturations[gaspos][local_index] = sg;
if ( water ) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
}
@ -712,42 +707,41 @@ namespace Opm
double so = 1.0;
double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 };
if (water) {
double swu = scaledDrainageInfo.Swu;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu);
so -= swu;
}
if (gas) {
double sgu = scaledDrainageInfo.Sgu;
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu);
so-= sgu;
}
fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
if (water) {
double swu = scaledDrainageInfo.Swu;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu);
so -= swu;
}
if (gas) {
double sgu = scaledDrainageInfo.Sgu;
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu);
so-= sgu;
}
fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
if (water && sw > scaledDrainageInfo.Swu-threshold_sat ) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pcWat;
} else if (gas && sg > scaledDrainageInfo.Sgu-threshold_sat) {
fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgu);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx];
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas;
}
if (gas && sg < scaledDrainageInfo.Sgl+threshold_sat) {
fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgl);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx];
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pcGas;
}
if (water && sw < scaledDrainageInfo.Swl+threshold_sat) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swl);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pcWat;
}
if (water && sw > scaledDrainageInfo.Swu-threshold_sat ) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pcWat;
} else if (gas && sg > scaledDrainageInfo.Sgu-threshold_sat) {
fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgu);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx];
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas;
}
if (gas && sg < scaledDrainageInfo.Sgl+threshold_sat) {
fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgl);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx];
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pcGas;
}
if (water && sw < scaledDrainageInfo.Swl+threshold_sat) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swl);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pcWat;
}
}
return phase_saturations;
}

8
tests/test_equil.cpp
View File

@ -210,10 +210,10 @@ BOOST_AUTO_TEST_CASE (PhasePressure)
initDefaultFluidSystem();
Opm::EQUIL::EquilReg
region(record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0);
region(record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0);
std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0);