OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1471 from totto82/implement_blackoil_energy_PR

Browse Source 
Add blackoil + energy specialization to Flow
This commit is contained in:
Andreas Lauser 2018-05-08 11:30:12 +02:00 committed by GitHub
commit 483e04c5da
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
10 changed files with 221 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -31,6 +31,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/flow_ebos_oilwater.cpp
opm/simulators/flow_ebos_polymer.cpp
opm/simulators/flow_ebos_solvent.cpp
opm/simulators/flow_ebos_energy.cpp
opm/simulators/flow_ebos_oilwater_polymer.cpp
opm/autodiff/Compat.cpp
@ -446,6 +447,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/flow_ebos_oilwater.hpp
opm/simulators/flow_ebos_polymer.hpp
opm/simulators/flow_ebos_solvent.hpp
opm/simulators/flow_ebos_energy.hpp
opm/simulators/flow_ebos_oilwater_polymer.hpp
opm/simulators/ensureDirectoryExists.hpp
opm/simulators/ParallelFileMerger.hpp

7
compareECLFiles.cmake
View File

@ -188,6 +188,13 @@ add_test_compareECLFiles(CASENAME spe1_nowells
REL_TOL ${rel_tol}
DIR spe1)
add_test_compareECLFiles(CASENAME spe1_thermal
FILENAME SPE1CASE2_THERMAL
SIMULATOR flow
ABS_TOL ${abs_tol}
REL_TOL ${rel_tol}
DIR spe1)
foreach(SIM flow flow_legacy)
add_test_compareECLFiles(CASENAME spe3
FILENAME SPE3CASE1

8
examples/flow.cpp
View File

@ -25,6 +25,7 @@
#include <opm/simulators/flow_ebos_oilwater.hpp>
#include <opm/simulators/flow_ebos_solvent.hpp>
#include <opm/simulators/flow_ebos_polymer.hpp>
#include <opm/simulators/flow_ebos_energy.hpp>
#include <opm/simulators/flow_ebos_oilwater_polymer.hpp>
#include <opm/autodiff/MissingFeatures.hpp>
@ -188,6 +189,11 @@ int main(int argc, char** argv)
Opm::flowEbosSolventSetDeck(*deck, *eclipseState, *schedule, *summary_config);
return Opm::flowEbosSolventMain(argc, argv);
}
// Energy case
else if ( phases.active( Opm::Phase::ENERGY ) ) {
Opm::flowEbosEnergySetDeck(*deck, *eclipseState, *schedule, *summary_config);
return Opm::flowEbosEnergyMain(argc, argv);
}
// Blackoil case
else if( phases.size() == 3 ) {
Opm::flowEbosBlackoilSetDeck(*deck, *eclipseState, *schedule, *summary_config);
@ -196,7 +202,7 @@ int main(int argc, char** argv)
else
{
if (outputCout)
std::cerr << "No suitable configuration found, valid are Twophase, polymer, solvent, or blackoil" << std::endl;
std::cerr << "No suitable configuration found, valid are Twophase, polymer, solvent, energy, or blackoil" << std::endl;
return EXIT_FAILURE;
}
}

21
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -81,6 +81,7 @@ SET_BOOL_PROP(EclFlowProblem, UseVolumetricResidual, false);
// necessary but it makes things a bit more explicit
SET_BOOL_PROP(EclFlowProblem, EnablePolymer, false);
SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false);
SET_BOOL_PROP(EclFlowProblem, EnableTemperature, true);
SET_BOOL_PROP(EclFlowProblem, EnableEnergy, false);
}}
@ -114,8 +115,10 @@ namespace Opm {
static const int numEq = Indices::numEq;
static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
static const int solventSaturationIdx = Indices::solventSaturationIdx;
static const int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
static const int temperatureIdx = Indices::temperatureIdx;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
@ -151,6 +154,7 @@ namespace Opm {
, has_vapoil_(FluidSystem::enableVaporizedOil())
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
, has_energy_(GET_PROP_VALUE(TypeTag, EnableEnergy))
, param_( param )
, well_model_ (well_model)
, terminal_output_ (terminal_output)
@ -713,6 +717,12 @@ namespace Opm {
c = std::max(c, 0.0);
}
if (has_energy_) {
double& T = priVars[Indices::temperatureIdx];
const double dT = dx[cell_idx][Indices::temperatureIdx];
T -= dT;
}
// Update rs and rv
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ) {
unsigned pvtRegionIdx = ebosSimulator_.problem().pvtRegionIndex(cell_idx);
@ -888,6 +898,12 @@ namespace Opm {
R_sum[ contiPolymerEqIdx ] += R2;
maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue );
}
if (has_energy_ ) {
B_avg[ contiEnergyEqIdx ] += 1.0;
const auto R2 = ebosResid[cell_idx][contiEnergyEqIdx];
R_sum[ contiEnergyEqIdx ] += R2;
maxCoeff[ contiEnergyEqIdx ] = std::max( maxCoeff[ contiEnergyEqIdx ], std::abs( R2 ) / pvValue );
}
}
@ -956,6 +972,10 @@ namespace Opm {
key[ polymerConcentrationIdx ] = "P";
}
if (has_energy_) {
key[ temperatureIdx ] = "E";
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " MB(" + key[ compIdx ] + ") ";
}
@ -1053,6 +1073,7 @@ namespace Opm {
const bool has_vapoil_;
const bool has_solvent_;
const bool has_polymer_;
const bool has_energy_;
ModelParameters param_;
SimulatorReport failureReport_;

18
opm/autodiff/StandardWell.hpp
View File

@ -53,6 +53,10 @@ namespace Opm
using Base::numEq;
using Base::has_solvent;
using Base::has_polymer;
using Base::has_energy;
// the positions of the primary variables for StandardWell
// there are three primary variables, the second and the third ones are F_w and F_g
// the first one can be total rate (G_t) or bhp, based on the control
@ -61,22 +65,22 @@ namespace Opm
static const int XvarWell = 0;
static const int WFrac = gasoil? -1000: 1;
static const int GFrac = gasoil? 1: 2;
static const int SFrac = 3;
static const int SFrac = !has_solvent ? -1000 : 3;
using typename Base::Scalar;
using typename Base::ConvergenceReport;
using Base::has_solvent;
using Base::has_polymer;
using Base::name;
using Base::Water;
using Base::Oil;
using Base::Gas;
// TODO: with flow_ebos，for a 2P deck, // TODO: for the 2p deck, numEq will be 3, a dummy phase is already added from the reservoir side.
// it will cause problem here without processing the dummy phase.
static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // number of wellEq is only numEq - 1 for polymer
// polymer concentration and temperature are already known by the well, so
// polymer and energy conservation do not need to be considered explicitly
static const int numPolymerEq = has_polymer ? 1 : 0;
static const int numEnergyEq = has_energy ? 1 : 0;
static const int numWellEq =numEq - numPolymerEq - numEnergyEq;
using typename Base::Mat;
using typename Base::BVector;
using typename Base::Eval;
@ -102,7 +106,7 @@ namespace Opm
using Base::contiSolventEqIdx;
using Base::contiPolymerEqIdx;
static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
StandardWell(const Well* well, const int time_step, const Wells* wells,
const ModelParameters& param,

69
opm/autodiff/StandardWell_impl.hpp
View File

@ -634,6 +634,71 @@ namespace Opm
well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
}
}
if (has_energy) {
auto fs = intQuants.fluidState();
const int reportStepIdx = ebosSimulator.episodeIndex();
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
// convert to reservoar conditions
EvalWell cq_r_thermal = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
if(FluidSystem::waterPhaseIdx == phaseIdx)
cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx));
// remove dissolved gas and vapporized oil
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * g_g + rs * q_or * b_o
// d = 1.0 - rs * rv
const EvalWell d = extendEval(1.0 - fs.Rv() * fs.Rs());
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
if(FluidSystem::gasPhaseIdx == phaseIdx)
cq_r_thermal = (cq_s[gasCompIdx] - extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) );
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
if(FluidSystem::oilPhaseIdx == phaseIdx)
cq_r_thermal = (cq_s[oilCompIdx] - extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * extendEval(fs.invB(phaseIdx)) );
} else {
cq_r_thermal = cq_s[activeCompIdx] / extendEval(fs.invB(phaseIdx));
}
// change temperature for injecting fluids
if (well_type_ == INJECTOR && cq_s[activeCompIdx] > 0.0){
const auto& injProps = this->well_ecl_->getInjectionProperties(reportStepIdx);
fs.setTemperature(injProps.temperature);
typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
typename FluidSystem::template ParameterCache<FsScalar> paramCache;
const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
paramCache.setRegionIndex(pvtRegionIdx);
paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx));
paramCache.updatePhase(fs, phaseIdx);
const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
fs.setDensity(phaseIdx, rho);
const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
fs.setEnthalpy(phaseIdx, h);
}
// compute the thermal flux
cq_r_thermal *= extendEval(fs.enthalpy(phaseIdx)) * extendEval(fs.density(phaseIdx));
// scale the flux by the scaling factor for the energy equation
cq_r_thermal *= GET_PROP_VALUE(TypeTag, BlackOilEnergyScalingFactor);
if (!only_wells) {
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
ebosJac[cell_idx][cell_idx][contiEnergyEqIdx][pvIdx] -= cq_r_thermal.derivative(pvIdx);
}
ebosResid[cell_idx][contiEnergyEqIdx] -= cq_r_thermal.value();
}
}
}
if (has_polymer) {
// TODO: the application of well efficiency factor has not been tested with an example yet
@ -1430,8 +1495,8 @@ namespace Opm
getWellConvergence(const std::vector<double>& B_avg) const
{
// the following implementation assume that the polymer is always after the w-o-g phases
// For the polymer case, there is one more mass balance equations of reservoir than wells
assert((int(B_avg.size()) == num_components_) || has_polymer);
// For the polymer case and the energy case, there is one more mass balance equations of reservoir than wells
assert((int(B_avg.size()) == num_components_) || has_polymer || has_energy);
const double tol_wells = param_.tolerance_wells_;
const double maxResidualAllowed = param_.max_residual_allowed_;

1
opm/autodiff/WellInterface.hpp
View File

@ -92,6 +92,7 @@ namespace Opm
static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
static const bool has_energy = GET_PROP_VALUE(TypeTag, EnableEnergy);
static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;

65
opm/simulators/flow_ebos_energy.cpp Normal file
View File

@ -0,0 +1,65 @@
/*
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 "config.h"
#include <opm/simulators/flow_ebos_energy.hpp>
#include <opm/material/common/ResetLocale.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
#include <opm/autodiff/FlowMainEbos.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
#include <dune/common/parallel/mpihelper.hh>
#endif
namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(EclFlowEnergyProblem, INHERITS_FROM(EclFlowProblem));
SET_BOOL_PROP(EclFlowEnergyProblem, EnableEnergy, true);
}}
namespace Opm {
void flowEbosEnergySetDeck(Deck &deck, EclipseState& eclState, Schedule& schedule, SummaryConfig& summaryConfig)
{
typedef TTAG(EclFlowEnergyProblem) TypeTag;
typedef GET_PROP_TYPE(TypeTag, Vanguard) Vanguard;
Vanguard::setExternalDeck(&deck, &eclState, &schedule, &summaryConfig);
}
// ----------------- Main program -----------------
int flowEbosEnergyMain(int argc, char** argv)
{
// we always want to use the default locale, and thus spare us the trouble
// with incorrect locale settings.
Opm::resetLocale();
// initialize MPI, finalize is done automatically on exit
#if HAVE_DUNE_FEM
Dune::Fem::MPIManager::initialize(argc, argv);
#else
Dune::MPIHelper::instance(argc, argv).rank();
#endif
Opm::FlowMainEbos<TTAG(EclFlowEnergyProblem)> mainfunc;
return mainfunc.execute(argc, argv);
}
}

30
opm/simulators/flow_ebos_energy.hpp Normal file
View File

@ -0,0 +1,30 @@
/*
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/>.
*/
#ifndef FLOW_EBOS_ENERGY_HPP
#define FLOW_EBOS_ENERGY_HPP
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
namespace Opm {
void flowEbosEnergySetDeck(Deck &deck, EclipseState& eclState, Schedule& schedule, SummaryConfig& summaryConfig);
int flowEbosEnergyMain(int argc, char** argv);
}
#endif // FLOW_EBOS_ENERGY_HPP

11
tests/update_reference_data.sh
View File

@ -20,7 +20,7 @@ copyToReferenceDir () {
}
tests=${@:2}
test -z "$tests" && tests="spe11 spe12 spe12p spe1oilgas spe1nowells spe3 spe5 spe9 norne_init msw_2d_h msw_3d_hfa polymer2d spe9group polymer_oilwater"
test -z "$tests" && tests="spe11 spe12 spe12p spe1oilgas spe1nowells spe1thermal spe3 spe5 spe9 norne_init msw_2d_h msw_3d_hfa polymer2d spe9group polymer_oilwater"
if grep -q -i "norne " <<< $ghprbCommentBody
then
if test -d $WORKSPACE/deps/opm-data/norne/flow
@ -89,6 +89,15 @@ for test_name in ${tests}; do
EGRID INIT SMSPEC UNRST UNSMRY
fi
if grep -q "spe1thermal" <<< $test_name
then
copyToReferenceDir \
$configuration/build-opm-simulators/tests/results/flow+spe1_thermal/ \
$OPM_TESTS_ROOT/spe1/opm-simulation-reference/flow \
SPE1CASE2_THERMAL \
EGRID INIT SMSPEC UNRST UNSMRY
fi
if grep -q "msw_2d_h" <<< $test_name
then
copyToReferenceDir \