mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-26 17:20:59 -06:00
Cleans up the Aquifer model classes.
This commit is contained in:
kel85uk
parent
2fe24e16cf
commit
b7fbe66c91
@ -24,21 +24,14 @@
|
||||
#include <opm/parser/eclipse/EclipseState/AquiferCT.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
|
||||
#include <opm/autodiff/BlackoilAquiferModel.hpp>
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
#include <opm/common/utility/numeric/linearInterpolation.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
|
||||
#include <opm/material/densead/Math.hpp>
|
||||
#include <opm/material/densead/Evaluation.hpp>
|
||||
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
|
||||
|
||||
#include <string>
|
||||
#include <memory>
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
#include <map>
|
||||
#include <cassert>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
@ -48,17 +41,11 @@ namespace Opm
|
||||
{
|
||||
|
||||
public:
|
||||
typedef BlackoilModelParameters ModelParameters;
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
|
||||
typedef typename GridView::template Codim<0>::Entity Element;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||
|
||||
static const int numEq = BlackoilIndices::numEq;
|
||||
typedef double Scalar;
|
||||
@ -71,73 +58,56 @@ namespace Opm
|
||||
|
||||
|
||||
|
||||
explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection,
|
||||
const Scalar gravity, const Simulator& ebosSimulator )
|
||||
AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data,
|
||||
const Aquancon::AquanconOutput& connection,
|
||||
Simulator& ebosSimulator )
|
||||
: ebos_simulator_ (ebosSimulator),
|
||||
aquiferID_ (params.aquiferID),
|
||||
inftableID_ (params.inftableID),
|
||||
pvttableID_ (params.pvttableID),
|
||||
phi_aq_ (params.phi_aq), //
|
||||
d0_ (params.d0),
|
||||
C_t_ (params.C_t), //
|
||||
r_o_ (params.r_o), //
|
||||
k_a_ (params.k_a), //
|
||||
c1_ (params.c1),
|
||||
h_ (params.h), //
|
||||
theta_ (params.theta), //
|
||||
c2_ (params.c2), //
|
||||
aqutab_td_ (params.td),
|
||||
aqutab_pi_ (params.pi),
|
||||
pa0_ (params.p0),
|
||||
gravity_ (gravity),
|
||||
p0_defaulted_ (params.p0_defaulted)
|
||||
aquct_data_ (aquct_data),
|
||||
gravity_ (ebos_simulator_.problem().gravity()[2])
|
||||
{
|
||||
init_quantities(connection);
|
||||
initQuantities(connection);
|
||||
}
|
||||
|
||||
inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer)
|
||||
inline void assembleAquiferEq(const SimulatorTimerInterface& timer)
|
||||
{
|
||||
dt_ = timer.currentStepLength();
|
||||
auto& ebosJac = ebosSimulator.model().linearizer().matrix();
|
||||
auto& ebosResid = ebosSimulator.model().linearizer().residual();
|
||||
auto& ebosJac = ebos_simulator_.model().linearizer().matrix();
|
||||
auto& ebosResid = ebos_simulator_.model().linearizer().residual();
|
||||
|
||||
|
||||
auto cellID = cell_idx_.begin();
|
||||
|
||||
size_t idx;
|
||||
for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx )
|
||||
size_t cellID;
|
||||
for ( size_t idx = 0; idx < cell_idx_.size(); ++idx )
|
||||
{
|
||||
Eval qinflow = 0.0;
|
||||
cellID = cell_idx_.at(idx);
|
||||
// We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to
|
||||
// IntensiveQuantities of that particular cell_id
|
||||
const IntensiveQuantities intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0));
|
||||
const IntensiveQuantities intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellID, /*timeIdx=*/ 0));
|
||||
// This is the pressure at td + dt
|
||||
get_current_Pressure_cell(pressure_current_,idx,intQuants);
|
||||
get_current_density_cell(rhow_,idx,intQuants);
|
||||
calculate_inflow_rate(idx, timer);
|
||||
updateCellPressure(pressure_current_,idx,intQuants);
|
||||
updateCellDensity(idx,intQuants);
|
||||
calculateInflowRate(idx, timer);
|
||||
qinflow = Qai_.at(idx);
|
||||
ebosResid[*cellID][waterCompIdx] -= qinflow.value();
|
||||
ebosResid[cellID][waterCompIdx] -= qinflow.value();
|
||||
|
||||
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
|
||||
{
|
||||
// also need to consider the efficiency factor when manipulating the jacobians.
|
||||
ebosJac[*cellID][*cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx);
|
||||
ebosJac[cellID][cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
inline void before_time_step(Simulator& ebosSimulator, const SimulatorTimerInterface& timer)
|
||||
inline void beforeTimeStep(const SimulatorTimerInterface& timer)
|
||||
{
|
||||
auto cellID = cell_idx_.begin();
|
||||
size_t idx;
|
||||
for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx )
|
||||
{
|
||||
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0));
|
||||
get_current_Pressure_cell(pressure_previous_ ,idx,intQuants);
|
||||
const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0));
|
||||
updateCellPressure(pressure_previous_ ,idx,intQuants);
|
||||
}
|
||||
}
|
||||
|
||||
inline void after_time_step(const SimulatorTimerInterface& timer)
|
||||
inline void afterTimeStep(const SimulatorTimerInterface& timer)
|
||||
{
|
||||
for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai)
|
||||
{
|
||||
@ -146,139 +116,117 @@ namespace Opm
|
||||
}
|
||||
|
||||
|
||||
inline const std::vector<int> cell_id() const
|
||||
{
|
||||
return cell_idx_;
|
||||
}
|
||||
|
||||
inline const int& aquiferID() const
|
||||
{
|
||||
return aquiferID_;
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
const Simulator& ebos_simulator_;
|
||||
|
||||
|
||||
// Aquifer ID, and other IDs
|
||||
int aquiferID_, inftableID_, pvttableID_;
|
||||
Simulator& ebos_simulator_;
|
||||
|
||||
// Grid variables
|
||||
|
||||
std::vector<size_t> cell_idx_;
|
||||
std::vector<Scalar> faceArea_connected_;
|
||||
|
||||
// Quantities at each grid id
|
||||
std::vector<Scalar> cell_depth_;
|
||||
std::vector<Eval> pressure_previous_;
|
||||
std::vector<Scalar> pressure_previous_;
|
||||
std::vector<Eval> pressure_current_;
|
||||
std::vector<Eval> Qai_;
|
||||
std::vector<Eval> rhow_;
|
||||
std::vector<Scalar> alphai_;
|
||||
|
||||
// Variables constants
|
||||
const AquiferCT::AQUCT_data aquct_data_;
|
||||
|
||||
Scalar mu_w_ , //water viscosity
|
||||
phi_aq_ , //aquifer porosity
|
||||
d0_, // aquifer datum depth
|
||||
C_t_ , //total compressibility
|
||||
r_o_ , //aquifer inner radius
|
||||
k_a_ , //aquifer permeability
|
||||
c1_, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB)
|
||||
h_ , //aquifer thickness
|
||||
theta_ , //angle subtended by the aquifer boundary
|
||||
c2_ ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB).
|
||||
beta_ , // Influx constant
|
||||
Tc_ , // Time constant
|
||||
pa0_ , // initial aquifer pressure
|
||||
gravity_ ; // gravitational acceleration
|
||||
|
||||
// Variables for influence table
|
||||
std::vector<Scalar> aqutab_td_, aqutab_pi_;
|
||||
|
||||
// Cumulative flux
|
||||
Scalar dt_, pa0_, gravity_;
|
||||
bool p0_defaulted_;
|
||||
Eval W_flux_;
|
||||
|
||||
|
||||
inline const double area_fraction(const size_t i)
|
||||
{
|
||||
return alphai_.at(i);
|
||||
}
|
||||
|
||||
inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
|
||||
inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
|
||||
{
|
||||
// We use the opm-common numeric linear interpolator
|
||||
pitd = Opm::linearInterpolation(aqutab_td_, aqutab_pi_, td);
|
||||
pitd_prime = Opm::linearInterpolationDerivative(aqutab_td_, aqutab_pi_, td);
|
||||
pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td);
|
||||
pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td);
|
||||
}
|
||||
|
||||
inline void init_quantities(const Aquancon::AquanconOutput& connection)
|
||||
inline void initQuantities(const Aquancon::AquanconOutput& connection)
|
||||
{
|
||||
// We reset the cumulative flux at the start of any simulation, so, W_flux = 0
|
||||
W_flux_ = 0.;
|
||||
|
||||
// We next get our connections to the aquifer and initialize these quantities using the initialize_connections function
|
||||
initialize_connections(connection);
|
||||
initializeConnections(connection);
|
||||
|
||||
calculate_aquifer_condition();
|
||||
calculateAquiferCondition();
|
||||
|
||||
calculateAquiferConstants();
|
||||
|
||||
pressure_previous_.resize(cell_idx_.size(), 0.);
|
||||
pressure_current_.resize(cell_idx_.size(), 0.);
|
||||
Qai_.resize(cell_idx_.size(), 0.0);
|
||||
}
|
||||
|
||||
inline void get_current_Pressure_cell(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
|
||||
inline void updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
|
||||
{
|
||||
const auto& fs = intQuants.fluidState();
|
||||
pressure_water.at(idx) = fs.pressure(waterPhaseIdx);
|
||||
}
|
||||
|
||||
inline void get_current_density_cell(std::vector<Eval>& rho_water, const int idx, const IntensiveQuantities& intQuants)
|
||||
inline void updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
|
||||
{
|
||||
const auto& fs = intQuants.fluidState();
|
||||
rho_water.at(idx) = fs.density(waterPhaseIdx);
|
||||
pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value();
|
||||
}
|
||||
|
||||
inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants)
|
||||
{
|
||||
const auto& fs = intQuants.fluidState();
|
||||
rhow_.at(idx) = fs.density(waterPhaseIdx);
|
||||
}
|
||||
|
||||
inline Scalar dpai(int idx)
|
||||
{
|
||||
Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_) - pressure_previous_.at(idx).value();
|
||||
Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0) - pressure_previous_.at(idx);
|
||||
return dp;
|
||||
}
|
||||
|
||||
// This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
|
||||
inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer)
|
||||
inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer)
|
||||
{
|
||||
Scalar beta = aquifer_influx_constant();
|
||||
Scalar Tc = time_constant();
|
||||
Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc;
|
||||
Scalar td = timer.simulationTimeElapsed() / Tc;
|
||||
const Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc_;
|
||||
const Scalar td = timer.simulationTimeElapsed() / Tc_;
|
||||
Scalar PItdprime = 0.;
|
||||
Scalar PItd = 0.;
|
||||
get_influence_table_values(PItd, PItdprime, td_plus_dt);
|
||||
a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime );
|
||||
b = beta / (Tc * ( PItd - td*PItdprime));
|
||||
getInfluenceTableValues(PItd, PItdprime, td_plus_dt);
|
||||
a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime );
|
||||
b = beta_ / (Tc_ * ( PItd - td*PItdprime));
|
||||
}
|
||||
|
||||
// This function implements Eq 5.7 of the EclipseTechnicalDescription
|
||||
inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer)
|
||||
inline void calculateInflowRate(int idx, const SimulatorTimerInterface& timer)
|
||||
{
|
||||
Scalar a, b;
|
||||
calculate_a_b_constants(a,b,idx,timer);
|
||||
Qai_.at(idx) = area_fraction(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx).value() ) );
|
||||
calculateEqnConstants(a,b,idx,timer);
|
||||
Qai_.at(idx) = alphai_.at(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx) ) );
|
||||
}
|
||||
|
||||
inline const Scalar time_constant() const
|
||||
inline void calculateAquiferConstants()
|
||||
{
|
||||
Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_);
|
||||
return Tc;
|
||||
}
|
||||
|
||||
inline const Scalar aquifer_influx_constant() const
|
||||
{
|
||||
Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_;
|
||||
return beta;
|
||||
// We calculate the influx constant
|
||||
beta_ = aquct_data_.c2 * aquct_data_.h
|
||||
* aquct_data_.theta * aquct_data_.phi_aq
|
||||
* aquct_data_.C_t
|
||||
* aquct_data_.r_o * aquct_data_.r_o;
|
||||
// We calculate the time constant
|
||||
Tc_ = mu_w_ * aquct_data_.phi_aq
|
||||
* aquct_data_.C_t
|
||||
* aquct_data_.r_o * aquct_data_.r_o
|
||||
/ ( aquct_data_.k_a * aquct_data_.c1 );
|
||||
}
|
||||
|
||||
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
|
||||
inline void initialize_connections(const Aquancon::AquanconOutput& connection)
|
||||
inline void initializeConnections(const Aquancon::AquanconOutput& connection)
|
||||
{
|
||||
const auto& eclState = ebos_simulator_.vanguard().eclState();
|
||||
const auto& ugrid = ebos_simulator_.vanguard().grid();
|
||||
@ -293,7 +241,7 @@ namespace Opm
|
||||
assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) );
|
||||
|
||||
// We hack the cell depth values for now. We can actually get it from elementcontext pos
|
||||
cell_depth_.resize(cell_idx_.size(), d0_);
|
||||
cell_depth_.resize(cell_idx_.size(), aquct_data_.d0);
|
||||
alphai_.resize(cell_idx_.size(), 1.0);
|
||||
faceArea_connected_.resize(cell_idx_.size(),0.0);
|
||||
Scalar faceArea;
|
||||
@ -301,14 +249,6 @@ namespace Opm
|
||||
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
|
||||
auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid);
|
||||
|
||||
for (auto influxCoeff: connection.influx_coeff){
|
||||
std::cout << "influx_coeff = " << influxCoeff << std::endl;
|
||||
}
|
||||
|
||||
for (auto influxMult: connection.influx_multiplier){
|
||||
std::cout << "influx_multiplier = " << influxMult << std::endl;
|
||||
}
|
||||
|
||||
// Translate the C face tag into the enum used by opm-parser's TransMult class
|
||||
Opm::FaceDir::DirEnum faceDirection;
|
||||
|
||||
@ -326,19 +266,22 @@ namespace Opm
|
||||
// the logically-Cartesian direction of the face
|
||||
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
|
||||
|
||||
|
||||
if (faceTag == 0) // left
|
||||
faceDirection = Opm::FaceDir::XMinus;
|
||||
else if (faceTag == 1) // right
|
||||
faceDirection = Opm::FaceDir::XPlus;
|
||||
else if (faceTag == 2) // back
|
||||
faceDirection = Opm::FaceDir::YMinus;
|
||||
else if (faceTag == 3) // front
|
||||
faceDirection = Opm::FaceDir::YPlus;
|
||||
else if (faceTag == 4) // bottom
|
||||
faceDirection = Opm::FaceDir::ZMinus;
|
||||
else if (faceTag == 5) // top
|
||||
faceDirection = Opm::FaceDir::ZPlus;
|
||||
switch(faceTag)
|
||||
{
|
||||
case 0: faceDirection = Opm::FaceDir::XMinus;
|
||||
break;
|
||||
case 1: faceDirection = Opm::FaceDir::XPlus;
|
||||
break;
|
||||
case 2: faceDirection = Opm::FaceDir::YMinus;
|
||||
break;
|
||||
case 3: faceDirection = Opm::FaceDir::YPlus;
|
||||
break;
|
||||
case 4: faceDirection = Opm::FaceDir::ZMinus;
|
||||
break;
|
||||
case 5: faceDirection = Opm::FaceDir::ZPlus;
|
||||
break;
|
||||
default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined");
|
||||
}
|
||||
|
||||
if (faceDirection == connection.reservoir_face_dir.at(idx))
|
||||
{
|
||||
@ -360,16 +303,20 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
|
||||
inline void calculate_aquifer_condition()
|
||||
inline void calculateAquiferCondition()
|
||||
{
|
||||
|
||||
int pvttableIdx = pvttableID_ - 1;
|
||||
int pvttableIdx = aquct_data_.pvttableID - 1;
|
||||
|
||||
rhow_.resize(cell_idx_.size(),0.);
|
||||
|
||||
if (p0_defaulted_)
|
||||
if (aquct_data_.p0 < 1.0)
|
||||
{
|
||||
pa0_ = calculate_reservoir_equilibrium(rhow_);
|
||||
pa0_ = calculateReservoirEquilibrium();
|
||||
}
|
||||
else
|
||||
{
|
||||
pa0_ = aquct_data_.p0;
|
||||
}
|
||||
|
||||
// Initialize a FluidState object first
|
||||
@ -388,10 +335,11 @@ namespace Opm
|
||||
}
|
||||
|
||||
// This function is for calculating the aquifer properties from equilibrium state with the reservoir
|
||||
inline Scalar calculate_reservoir_equilibrium(std::vector<Eval>& rho_water_reservoir)
|
||||
inline Scalar calculateReservoirEquilibrium()
|
||||
{
|
||||
// Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
|
||||
std::vector<Scalar> water_pressure_reservoir, pw_aquifer;
|
||||
std::vector<Scalar> pw_aquifer;
|
||||
Scalar water_pressure_reservoir;
|
||||
|
||||
for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
|
||||
{
|
||||
@ -399,9 +347,9 @@ namespace Opm
|
||||
const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0));
|
||||
const auto& fs = intQuants.fluidState();
|
||||
|
||||
water_pressure_reservoir.push_back( fs.pressure(waterPhaseIdx).value() );
|
||||
rho_water_reservoir.at(idx) = fs.density(waterPhaseIdx);
|
||||
pw_aquifer.push_back( (water_pressure_reservoir.at(idx) - rho_water_reservoir.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_))*area_fraction(idx) );
|
||||
water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
|
||||
rhow_.at(idx) = fs.density(waterPhaseIdx);
|
||||
pw_aquifer.push_back( (water_pressure_reservoir - rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0))*alphai_.at(idx) );
|
||||
}
|
||||
|
||||
// We take the average of the calculated equilibrium pressures.
|
||||
|
||||
@ -24,40 +24,12 @@
|
||||
#ifndef OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED
|
||||
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
|
||||
#include <opm/common/utility/platform_dependent/disable_warnings.h>
|
||||
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
|
||||
|
||||
#include <cassert>
|
||||
#include <tuple>
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/AquiferCT.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Aquancon.hpp>
|
||||
|
||||
#include <opm/core/simulator/SimulatorReport.hpp>
|
||||
|
||||
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
|
||||
|
||||
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
|
||||
#include <opm/autodiff/BlackoilDetails.hpp>
|
||||
#include <opm/autodiff/BlackoilModelParameters.hpp>
|
||||
#include <opm/autodiff/RateConverter.hpp>
|
||||
#include <opm/autodiff/AquiferCarterTracy.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
||||
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
|
||||
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
|
||||
|
||||
#include <dune/common/fmatrix.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/matrixmatrix.hh>
|
||||
|
||||
#include <opm/material/densead/Math.hpp>
|
||||
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/// Class for handling the blackoil well model.
|
||||
@ -67,22 +39,13 @@ namespace Opm {
|
||||
public:
|
||||
|
||||
// --------- Types ---------
|
||||
typedef BlackoilModelParameters ModelParameters;
|
||||
|
||||
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
|
||||
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
|
||||
|
||||
typedef AquiferCarterTracy<TypeTag> Aquifer_object;
|
||||
|
||||
BlackoilAquiferModel(Simulator& ebosSimulator,
|
||||
const ModelParameters& param,
|
||||
const bool terminal_output);
|
||||
explicit BlackoilAquiferModel(Simulator& ebosSimulator);
|
||||
|
||||
// compute the well fluxes and assemble them in to the reservoir equations as source terms
|
||||
// and in the well equations.
|
||||
@ -92,36 +55,22 @@ namespace Opm {
|
||||
// called at the end of a time step
|
||||
void timeStepSucceeded(const SimulatorTimerInterface& timer);
|
||||
|
||||
inline const Simulator& simulator() const
|
||||
{
|
||||
return ebosSimulator_;
|
||||
}
|
||||
|
||||
// This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
|
||||
void init(const Simulator& ebosSimulator, std::vector<Aquifer_object>& aquifers);
|
||||
|
||||
protected:
|
||||
|
||||
Simulator& ebosSimulator_;
|
||||
|
||||
const ModelParameters param_;
|
||||
bool terminal_output_;
|
||||
|
||||
double gravity_;
|
||||
std::vector<Aquifer_object> aquifers_;
|
||||
|
||||
// This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
|
||||
void init();
|
||||
|
||||
void updateConnectionIntensiveQuantities() const;
|
||||
|
||||
int numAquifers() const;
|
||||
|
||||
void assembleAquiferEq(const SimulatorTimerInterface& timer);
|
||||
|
||||
// at the beginning of each time step (Not report step)
|
||||
void prepareTimeStep(const SimulatorTimerInterface& timer);
|
||||
|
||||
const std::vector<Aquifer_object>& aquifers();
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -3,19 +3,10 @@ namespace Opm {
|
||||
|
||||
template<typename TypeTag>
|
||||
BlackoilAquiferModel<TypeTag>::
|
||||
BlackoilAquiferModel(Simulator& ebosSimulator,
|
||||
const ModelParameters& param,
|
||||
const bool terminal_output)
|
||||
BlackoilAquiferModel(Simulator& ebosSimulator)
|
||||
: ebosSimulator_(ebosSimulator)
|
||||
, param_(param)
|
||||
, terminal_output_(terminal_output)
|
||||
{
|
||||
// const auto& gridView = ebosSimulator_.gridView();
|
||||
|
||||
// number_of_cells_ = gridView.size(/*codim=*/0);
|
||||
// global_nc_ = gridView.comm().sum(number_of_cells_);
|
||||
gravity_ = ebosSimulator_.problem().gravity()[2];
|
||||
init(ebosSimulator_, aquifers_);
|
||||
init();
|
||||
}
|
||||
|
||||
|
||||
@ -26,7 +17,7 @@ namespace Opm {
|
||||
{
|
||||
for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
|
||||
{
|
||||
aquifer->after_time_step(timer);
|
||||
aquifer->afterTimeStep(timer);
|
||||
}
|
||||
}
|
||||
|
||||
@ -65,16 +56,6 @@ namespace Opm {
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Protected function: Return number of aquifers in the model.
|
||||
template<typename TypeTag>
|
||||
int
|
||||
BlackoilAquiferModel<TypeTag>:: numAquifers() const
|
||||
{
|
||||
return aquifers_.size();
|
||||
}
|
||||
|
||||
|
||||
// Protected function which calls the individual aquifer models
|
||||
template<typename TypeTag>
|
||||
void
|
||||
@ -82,7 +63,7 @@ namespace Opm {
|
||||
{
|
||||
for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
|
||||
{
|
||||
aquifer->assembleAquiferEq(ebosSimulator_, timer);
|
||||
aquifer->assembleAquiferEq(timer);
|
||||
}
|
||||
}
|
||||
|
||||
@ -95,42 +76,33 @@ namespace Opm {
|
||||
// Here we can ask each carter tracy aquifers to get the current previous time step's pressure
|
||||
for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
|
||||
{
|
||||
aquifer->before_time_step(ebosSimulator_, timer);
|
||||
aquifer->beforeTimeStep(timer);
|
||||
}
|
||||
}
|
||||
|
||||
// Protected function: Returns a reference to the aquifers members in the model
|
||||
template<typename TypeTag>
|
||||
const std::vector< AquiferCarterTracy<TypeTag> >&
|
||||
BlackoilAquiferModel<TypeTag>:: aquifers()
|
||||
{
|
||||
return aquifers_;
|
||||
}
|
||||
|
||||
|
||||
// Initialize the aquifers in the deck
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy<TypeTag> >& aquifers)
|
||||
BlackoilAquiferModel<TypeTag>:: init()
|
||||
{
|
||||
updateConnectionIntensiveQuantities();
|
||||
const auto& deck = ebosSimulator.vanguard().deck();
|
||||
const auto& eclState = ebosSimulator.vanguard().eclState();
|
||||
const auto& deck = ebosSimulator_.vanguard().deck();
|
||||
const auto& eclState = ebosSimulator_.vanguard().eclState();
|
||||
|
||||
// Get all the carter tracy aquifer properties data and put it in aquifers vector
|
||||
AquiferCT aquiferct = AquiferCT(eclState,deck);
|
||||
Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
|
||||
const AquiferCT aquiferct = AquiferCT(eclState,deck);
|
||||
const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
|
||||
|
||||
std::vector<AquiferCT::AQUCT_data> aquifersData = aquiferct.getAquifers();
|
||||
std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput();
|
||||
|
||||
assert( aquifersData.size() == aquifer_connect.size() );
|
||||
assert( aquifersData.size() == aquifer_connection.size() );
|
||||
|
||||
|
||||
for (size_t i = 0; i < aquifersData.size(); ++i)
|
||||
{
|
||||
aquifers.push_back(
|
||||
AquiferCarterTracy<TypeTag> (aquifersData.at(i), aquifer_connection.at(i), gravity_, ebosSimulator_)
|
||||
aquifers_.push_back(
|
||||
AquiferCarterTracy<TypeTag> (aquifersData.at(i), aquifer_connection.at(i), ebosSimulator_)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
@ -30,8 +30,6 @@
|
||||
#include <opm/autodiff/BlackoilModelParameters.hpp>
|
||||
#include <opm/autodiff/BlackoilWellModel.hpp>
|
||||
#include <opm/autodiff/BlackoilAquiferModel.hpp>
|
||||
#include <opm/autodiff/GridHelpers.hpp>
|
||||
#include <opm/autodiff/GeoProps.hpp>
|
||||
#include <opm/autodiff/WellConnectionAuxiliaryModule.hpp>
|
||||
#include <opm/autodiff/BlackoilDetails.hpp>
|
||||
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
|
||||
@ -349,7 +347,7 @@ namespace Opm {
|
||||
const ReservoirState& reservoir_state,
|
||||
WellState& well_state)
|
||||
{
|
||||
// DUNE_UNUSED_PARAMETER(timer);
|
||||
DUNE_UNUSED_PARAMETER(timer);
|
||||
DUNE_UNUSED_PARAMETER(reservoir_state);
|
||||
DUNE_UNUSED_PARAMETER(well_state);
|
||||
|
||||
@ -372,8 +370,8 @@ namespace Opm {
|
||||
ebosSimulator_.model().linearizer().linearize();
|
||||
ebosSimulator_.problem().endIteration();
|
||||
|
||||
// -------- Well and aquifer common variables ----------
|
||||
double dt = timer.currentStepLength();
|
||||
// -------- Current time step length ----------
|
||||
const double dt = timer.currentStepLength();
|
||||
|
||||
// -------- Aquifer models ----------
|
||||
try
|
||||
@ -381,9 +379,9 @@ namespace Opm {
|
||||
// Modify the Jacobian and residuals according to the aquifer models
|
||||
aquiferModel().assemble(timer, iterationIdx);
|
||||
}
|
||||
catch( const Dune::FMatrixError& e )
|
||||
catch( ... )
|
||||
{
|
||||
OPM_THROW(Opm::NumericalProblem,"Error when assembling aquifer models");
|
||||
OPM_THROW(Opm::NumericalIssue,"Error when assembling aquifer models");
|
||||
}
|
||||
|
||||
// -------- Well equations ----------
|
||||
|
||||
@ -30,8 +30,6 @@
|
||||
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
|
||||
#include <opm/autodiff/BlackoilWellModel.hpp>
|
||||
#include <opm/autodiff/BlackoilAquiferModel.hpp>
|
||||
#include <opm/autodiff/RateConverter.hpp>
|
||||
#include <opm/autodiff/SimFIBODetails.hpp>
|
||||
#include <opm/autodiff/moduleVersion.hpp>
|
||||
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
|
||||
#include <opm/grid/utility/StopWatch.hpp>
|
||||
@ -191,7 +189,7 @@ public:
|
||||
ebosSimulator_.model().addAuxiliaryModule(auxMod);
|
||||
}
|
||||
|
||||
AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_);
|
||||
AquiferModel aquifer_model(ebosSimulator_);
|
||||
|
||||
// Main simulation loop.
|
||||
while (!timer.done()) {
|
||||
|
||||
Loading…
Reference in New Issue
Block a user