mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-26 01:01:00 -06:00
First addition of the class BlackoilAquiferModel.
This commit is contained in:
kel85uk
parent
05a02b9d68
commit
10e896fa6b
402
opm/autodiff/AquiferCarterTracy.hpp
Normal file
402
opm/autodiff/AquiferCarterTracy.hpp
Normal file
@ -0,0 +1,402 @@
|
||||
/*
|
||||
Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface
|
||||
Copyright 2017 Statoil ASA.
|
||||
|
||||
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 OPM_AQUIFERCT_HEADER_INCLUDED
|
||||
#define OPM_AQUIFERCT_HEADER_INCLUDED
|
||||
|
||||
#include <Eigen/QR>
|
||||
#include <opm/parser/eclipse/EclipseState/AquiferCT.hpp>
|
||||
#include <opm/autodiff/BlackoilAquiferModel.hpp>
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
|
||||
#include <opm/material/densead/Math.hpp>
|
||||
#include <opm/material/densead/Evaluation.hpp>
|
||||
|
||||
#include <string>
|
||||
#include <memory>
|
||||
#include <vector>
|
||||
#include <map>
|
||||
#include <cassert>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
template<typename TypeTag>
|
||||
class AquiferCarterTracy
|
||||
{
|
||||
|
||||
public:
|
||||
typedef BlackoilModelParameters ModelParameters;
|
||||
|
||||
static const int Water = BlackoilPhases::Aqua;
|
||||
static const int Oil = BlackoilPhases::Liquid;
|
||||
static const int Gas = BlackoilPhases::Vapour;
|
||||
|
||||
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;
|
||||
|
||||
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
|
||||
|
||||
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
|
||||
|
||||
static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
|
||||
static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
|
||||
static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
|
||||
static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
|
||||
|
||||
|
||||
AquiferCarterTracy(const std::vector<int>& cell_id)
|
||||
: phi_aq_ (1.0), //
|
||||
C_t_ (1.0), //
|
||||
r_o_ (1.0), //
|
||||
k_a_ (1.0), //
|
||||
c1_ (1.0),
|
||||
h_ (1.0), //
|
||||
theta_ (1.0), //
|
||||
c2_ (1.0), //
|
||||
d0_ (1.0),
|
||||
cell_idx_ (cell_id)
|
||||
{
|
||||
mu_w_ = 1e-3;
|
||||
aqutab_td_.push_back(1.0);
|
||||
aqutab_pi_.push_back(1.0);
|
||||
aquiferID_ = 1;
|
||||
inftableID_ = 1;
|
||||
pvttableID_ = 1;
|
||||
init_quantities();
|
||||
}
|
||||
|
||||
explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const AquiferCT::AQUANCON_data& aquanconParams,
|
||||
const int numComponents, const Scalar gravity )
|
||||
: phi_aq_ (params.phi_aq), //
|
||||
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), //
|
||||
d0_ (params.d0),
|
||||
aqutab_td_ (params.td),
|
||||
aqutab_pi_ (params.pi),
|
||||
aquiferID_ (params.aquiferID),
|
||||
inftableID_ (params.inftableID),
|
||||
pvttableID_ (params.pvttableID),
|
||||
cell_idx_ (params.cell_id),
|
||||
num_components_ (numComponents),
|
||||
gravity_ (gravity)
|
||||
{
|
||||
mu_w_ = 1e-3;
|
||||
init_quantities(aquanconParams);
|
||||
}
|
||||
|
||||
inline const PhaseUsage&
|
||||
phaseUsage() const
|
||||
{
|
||||
assert(phase_usage_);
|
||||
|
||||
return *phase_usage_;
|
||||
}
|
||||
|
||||
inline int
|
||||
flowPhaseToEbosCompIdx( const int phaseIdx ) const
|
||||
{
|
||||
const auto& pu = phaseUsage();
|
||||
if (pu.phase_pos[Water] == phaseIdx)
|
||||
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
|
||||
if (pu.phase_pos[Oil] == phaseIdx)
|
||||
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
|
||||
if (pu.phase_pos[Gas] == phaseIdx)
|
||||
return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
|
||||
|
||||
// for other phases return the index
|
||||
return phaseIdx;
|
||||
}
|
||||
|
||||
inline int
|
||||
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
|
||||
{
|
||||
const auto& pu = phaseUsage();
|
||||
if (pu.phase_pos[Water] == phaseIdx) {
|
||||
return FluidSystem::waterPhaseIdx;
|
||||
}
|
||||
if (pu.phase_pos[Oil] == phaseIdx) {
|
||||
return FluidSystem::oilPhaseIdx;
|
||||
}
|
||||
if (pu.phase_pos[Gas] == phaseIdx) {
|
||||
return FluidSystem::gasPhaseIdx;
|
||||
}
|
||||
|
||||
assert(phaseIdx < 3);
|
||||
// for other phases return the index
|
||||
return phaseIdx;
|
||||
}
|
||||
|
||||
inline void calculateExplicitQuantities(const Simulator& ebosSimulator)
|
||||
{
|
||||
std::cout << "In CarterTracy<calculateExplicitQuantities>: I am aquifer #" << aquiferID_ << std::endl;
|
||||
}
|
||||
|
||||
inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer)
|
||||
{
|
||||
std::cout << "In CarterTracy<assembleAquiferEq>: I am aquifer #" << aquiferID_ << std::endl;
|
||||
// resAqui_ = 0.0;
|
||||
dt_ = timer.currentStepLength();
|
||||
auto& ebosJac = ebosSimulator.model().linearizer().matrix();
|
||||
auto& ebosResid = ebosSimulator.model().linearizer().residual();
|
||||
|
||||
// TODO: it probably can be static member for StandardWell
|
||||
const double volume = 0.002831684659200; // 0.1 cu ft;
|
||||
|
||||
auto cellID = cell_idx_.begin();
|
||||
size_t idx;
|
||||
for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx )
|
||||
{
|
||||
Eval qinflow = 0.0;
|
||||
// We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to
|
||||
// IntensiveQuantities of that particular cell_id
|
||||
const auto& intQuants = *(ebosSimulator.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);
|
||||
qinflow = Qai_[idx];
|
||||
ebosResid[*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)] -= qinflow.value();
|
||||
|
||||
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
|
||||
{
|
||||
// also need to consider the efficiency factor when manipulating the jacobians.
|
||||
ebosJac[*cellID][*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)][pvIdx] -= qinflow.derivative(pvIdx);
|
||||
}
|
||||
std::cout << "In CarterTracy<assembleAquiferEq>: I am aquifer #" << aquiferID_
|
||||
// << " -> P_wat[t+dt] = " << pressure_current_[idx] << std::endl
|
||||
<< " Qai[t+dt] = " << Qai_[idx] << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
inline void before_time_step(Simulator& ebosSimulator, 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);
|
||||
}
|
||||
}
|
||||
|
||||
inline void after_time_step()
|
||||
{
|
||||
for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai)
|
||||
{
|
||||
W_flux_ += (*Qai);
|
||||
}
|
||||
std::cout << "Aquifer # " << aquiferID_ << ": My cumulative flux = " << W_flux_ << std::endl;
|
||||
}
|
||||
|
||||
/* Made into public for testing only!!!!!!. Must be protected */
|
||||
inline const Scalar time_constant() const
|
||||
{
|
||||
Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_);
|
||||
return Tc;
|
||||
}
|
||||
|
||||
/* Made into public for testing only!!!!!!. Must be protected */
|
||||
inline const Scalar aquifer_influx_constant() const
|
||||
{
|
||||
Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_;
|
||||
return beta;
|
||||
}
|
||||
|
||||
// This is another hack to get the face area only for SPE1.
|
||||
// Ideally it should be a map which given a cell_id, it returns the area fraction
|
||||
inline const double area_fraction(const int i)
|
||||
{
|
||||
return 1000.0*20.0*0.092903/(1000.0*1000.0*0.092903*2 + 1000.0*20.0*0.092903*4);
|
||||
}
|
||||
|
||||
inline void print_private_members() const
|
||||
{
|
||||
std::cout << "Aquifer CT #" << aquiferID_ << std::endl;
|
||||
auto ita = aqutab_td_.cbegin();
|
||||
auto f_lambda = [&ita] (double i) {std::cout << *ita++ << " " << i << std::endl;};
|
||||
std::for_each( aqutab_pi_.cbegin(), aqutab_pi_.cend(), f_lambda );
|
||||
|
||||
for (auto i = coeff_.begin(); i != coeff_.end(); ++i )
|
||||
{
|
||||
std::cout << "Coeff = " << *i << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
/* Made into public for testing only!!!!!!. Must be protected */
|
||||
inline const std::vector<int> cell_id() const
|
||||
{
|
||||
return cell_idx_;
|
||||
}
|
||||
|
||||
inline const int& aquiferID() const
|
||||
{
|
||||
return aquiferID_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
protected:
|
||||
const PhaseUsage* phase_usage_;
|
||||
|
||||
|
||||
// Aquifer ID, and other IDs
|
||||
int aquiferID_, inftableID_, pvttableID_;
|
||||
int num_components_;
|
||||
|
||||
// Grid variables
|
||||
std::vector<int> cell_idx_;
|
||||
|
||||
// Quantities at each grid id
|
||||
std::vector<Scalar> cell_depth_;
|
||||
std::vector<Scalar> pressure_previous_;
|
||||
std::vector<Scalar> pressure_current_;
|
||||
std::vector<Scalar> Qai_;
|
||||
std::vector<Scalar> rhow_;
|
||||
|
||||
// Variables constants
|
||||
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).
|
||||
|
||||
// Variables for influence table
|
||||
std::vector<Scalar> aqutab_td_, aqutab_pi_;
|
||||
|
||||
// Cumulative flux
|
||||
Scalar W_flux_, dt_, pa0_, gravity_;
|
||||
|
||||
// Also return the polynomial fit
|
||||
std::vector<Scalar> coeff_;
|
||||
|
||||
|
||||
// We fit the tabular data using a polynomial fit
|
||||
// Modified from Copyright (C) 2014 Clifford Wolf <clifford@clifford.at>
|
||||
// http://svn.clifford.at/handicraft/2014/polyfit/polyfit.cc
|
||||
inline void polynomial_fit( const std::vector<Scalar> &X, const std::vector<Scalar> &y,
|
||||
std::vector<Scalar> &coeff, int order, bool bias) const
|
||||
{
|
||||
int colNum = (bias)? order + 1 : order;
|
||||
Eigen::MatrixXd A(X.size(), colNum);
|
||||
Eigen::VectorXd y_mapped = Eigen::VectorXd::Map(&y.front(), y.size());
|
||||
Eigen::VectorXd result;
|
||||
|
||||
assert(X.size() == y.size());
|
||||
assert(X.size() >= colNum);
|
||||
|
||||
// create matrix
|
||||
for (size_t i = 0; i < X.size(); i++)
|
||||
for (size_t j = 0; j < colNum; j++)
|
||||
A(i, j) = (bias)? pow(X.at(i), j) : pow(X.at(i), j+1);
|
||||
|
||||
// solve for linear least squares fit
|
||||
result = A.householderQr().solve(y_mapped);
|
||||
|
||||
coeff.resize(colNum);
|
||||
for (size_t i = 0; i < colNum; i++)
|
||||
coeff[i] = result[i];
|
||||
}
|
||||
|
||||
inline void init_quantities(const AquiferCT::AQUANCON_data& aquanconParams)
|
||||
{
|
||||
W_flux_ = 0.;
|
||||
// pa0_ is the initial aquifer water pressure. Must be calculated from equilibrium if left default,
|
||||
// or we get the information from the deck Hacked to make it at 45e6 Pa
|
||||
pa0_ = 45e6;
|
||||
|
||||
pressure_previous_.resize(cell_idx_.size(), 0.);
|
||||
pressure_current_.resize(cell_idx_.size(), 0.);
|
||||
// We hack the cell depth values for now. We can actually get it from elementcontext pos
|
||||
cell_depth_.resize(cell_idx_.size(), d0_);
|
||||
rhow_.resize(cell_idx_.size(), 998.0);
|
||||
Qai_.resize(cell_idx_.size(), 0.);
|
||||
|
||||
polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 2, true);
|
||||
}
|
||||
|
||||
inline void get_current_Pressure_cell(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
|
||||
{
|
||||
const auto& fs = intQuants.fluidState();
|
||||
pressure_water[idx] = fs.pressure(FluidSystem::waterPhaseIdx).value();
|
||||
}
|
||||
|
||||
inline void get_current_density_cell(std::vector<Scalar>& rho_water, const int idx, const IntensiveQuantities& intQuants)
|
||||
{
|
||||
const auto& fs = intQuants.fluidState();
|
||||
rho_water[idx] = fs.density(FluidSystem::waterPhaseIdx).value();
|
||||
}
|
||||
|
||||
inline Scalar dpai(int idx)
|
||||
{
|
||||
Scalar dp = pa0_ - rhow_[idx]*gravity_*(cell_depth_[idx] - d0_) - pressure_previous_[idx];
|
||||
return dp;
|
||||
}
|
||||
|
||||
inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer)
|
||||
{
|
||||
// This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
|
||||
Scalar beta = aquifer_influx_constant();
|
||||
Scalar Tc = time_constant();
|
||||
Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc;
|
||||
Scalar td = timer.simulationTimeElapsed() / Tc;
|
||||
Scalar PItdprime = coeff_[1] + 2.0*coeff_[2]*(td_plus_dt);
|
||||
Scalar PItd = coeff_[0] + coeff_[1]*td_plus_dt + coeff_[2]*td_plus_dt*td_plus_dt;
|
||||
a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_ * PItdprime) ) / ( PItd - td*PItdprime );
|
||||
b = beta / Tc / ( PItd - td*PItdprime);
|
||||
}
|
||||
|
||||
inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer)
|
||||
{
|
||||
Scalar a, b;
|
||||
calculate_a_b_constants(a,b,idx,timer);
|
||||
// This function implements Eq 5.7 of the EclipseTechnicalDescription
|
||||
Qai_[idx] = area_fraction(idx)*( a - b * ( pressure_current_[idx] - pressure_previous_[idx] ) );
|
||||
}
|
||||
|
||||
|
||||
}; // class AquiferCarterTracy
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
184
opm/autodiff/BlackoilAquiferModel.hpp
Normal file
184
opm/autodiff/BlackoilAquiferModel.hpp
Normal file
@ -0,0 +1,184 @@
|
||||
/*
|
||||
<<<<<<< HEAD
|
||||
File adapted from BlackoilWellModel.hpp
|
||||
|
||||
Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface
|
||||
Copyright 2017 Statoil ASA.
|
||||
=======
|
||||
Copyright 2016 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2016 - 2017 Statoil ASA.
|
||||
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
|
||||
Copyright 2016 - 2017 IRIS AS
|
||||
>>>>>>> 9ccee28... First addition of the class BlackoilAquiferModel.
|
||||
|
||||
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 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/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.
|
||||
template<typename TypeTag>
|
||||
class BlackoilAquiferModel {
|
||||
|
||||
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;
|
||||
|
||||
static const int numEq = BlackoilIndices::numEq;
|
||||
static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
|
||||
|
||||
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
|
||||
|
||||
typedef AquiferCarterTracy<TypeTag> Aquifer_object;
|
||||
|
||||
BlackoilAquiferModel(Simulator& ebosSimulator,
|
||||
const ModelParameters& param,
|
||||
const bool terminal_output);
|
||||
|
||||
// compute the well fluxes and assemble them in to the reservoir equations as source terms
|
||||
// and in the well equations.
|
||||
void assemble( const SimulatorTimerInterface& timer,
|
||||
const int iterationIdx );
|
||||
|
||||
// called at the beginning of a time step
|
||||
void beginTimeStep();
|
||||
// called at the end of a time step
|
||||
void timeStepSucceeded();
|
||||
|
||||
// called at the beginning of a report step
|
||||
void beginReportStep(const int time_step);
|
||||
|
||||
// called at the end of a report step
|
||||
void endReportStep();
|
||||
|
||||
const SimulatorReport& lastReport() const;
|
||||
|
||||
inline const Simulator& simulator() const
|
||||
{
|
||||
return ebosSimulator_;
|
||||
}
|
||||
|
||||
/// Hack function to get what I need from parser
|
||||
void init(const Simulator& ebosSimulator, std::vector<Aquifer_object>& aquifers);
|
||||
|
||||
protected:
|
||||
|
||||
Simulator& ebosSimulator_;
|
||||
|
||||
const ModelParameters param_;
|
||||
bool terminal_output_;
|
||||
bool has_solvent_;
|
||||
bool has_polymer_;
|
||||
std::vector<int> pvt_region_idx_;
|
||||
PhaseUsage phase_usage_;
|
||||
std::vector<bool> active_;
|
||||
size_t global_nc_;
|
||||
// the number of the cells in the local grid
|
||||
size_t number_of_cells_;
|
||||
double gravity_;
|
||||
std::vector<double> depth_;
|
||||
std::vector<Aquifer_object> aquifers_;
|
||||
|
||||
|
||||
SimulatorReport last_report_;
|
||||
|
||||
const Schedule& schedule() const
|
||||
{ return ebosSimulator_.gridManager().schedule(); }
|
||||
|
||||
void updatePrimaryVariables();
|
||||
|
||||
void initPrimaryVariablesEvaluation() const;
|
||||
|
||||
void updateConnectionIntensiveQuantities() const;
|
||||
|
||||
void calculateExplicitQuantities();
|
||||
|
||||
// The number of components in the model.
|
||||
int numComponents() const;
|
||||
|
||||
int numAquifers() const;
|
||||
|
||||
int numPhases() const;
|
||||
|
||||
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
|
||||
|
||||
void assembleAquiferEq(const SimulatorTimerInterface& timer);
|
||||
|
||||
SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer);
|
||||
|
||||
// some preparation work, mostly related to group control and RESV,
|
||||
// at the beginning of each time step (Not report step)
|
||||
void prepareTimeStep(const SimulatorTimerInterface& timer);
|
||||
|
||||
const std::vector<Aquifer_object>& aquifers();
|
||||
|
||||
};
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#include "BlackoilAquiferModel_impl.hpp"
|
||||
#endif
|
||||
287
opm/autodiff/BlackoilAquiferModel_impl.hpp
Normal file
287
opm/autodiff/BlackoilAquiferModel_impl.hpp
Normal file
@ -0,0 +1,287 @@
|
||||
namespace Opm {
|
||||
|
||||
|
||||
template<typename TypeTag>
|
||||
BlackoilAquiferModel<TypeTag>::
|
||||
BlackoilAquiferModel(Simulator& ebosSimulator,
|
||||
const ModelParameters& param,
|
||||
const bool terminal_output)
|
||||
: ebosSimulator_(ebosSimulator)
|
||||
, param_(param)
|
||||
, terminal_output_(terminal_output)
|
||||
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
|
||||
, has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
|
||||
{
|
||||
const auto& eclState = ebosSimulator_.gridManager().eclState();
|
||||
phase_usage_ = phaseUsageFromDeck(eclState);
|
||||
|
||||
active_.resize(phase_usage_.MaxNumPhases, false);
|
||||
for (int p = 0; p < phase_usage_.MaxNumPhases; ++p) {
|
||||
active_[ p ] = phase_usage_.phase_used[ p ] != 0;
|
||||
}
|
||||
|
||||
const auto& gridView = ebosSimulator_.gridView();
|
||||
|
||||
// calculate the number of elements of the compressed sequential grid. this needs
|
||||
// to be done in two steps because the dune communicator expects a reference as
|
||||
// argument for sum()
|
||||
number_of_cells_ = gridView.size(/*codim=*/0);
|
||||
global_nc_ = gridView.comm().sum(number_of_cells_);
|
||||
gravity_ = ebosSimulator_.problem().gravity()[2];
|
||||
init(ebosSimulator_, aquifers_);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// called at the beginning of a time step
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: beginTimeStep()
|
||||
{
|
||||
// Right now it doesn't do shit.
|
||||
}
|
||||
|
||||
// called at the end of a time step
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: timeStepSucceeded()
|
||||
{
|
||||
for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
|
||||
{
|
||||
aquifer->after_time_step();
|
||||
}
|
||||
}
|
||||
|
||||
// called at the beginning of a report step
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: beginReportStep(const int time_step)
|
||||
{
|
||||
// Right now it doesn't do shit.
|
||||
}
|
||||
|
||||
// called at the end of a report step
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: endReportStep()
|
||||
{
|
||||
// Right now it just spits out the constants for each aquifers
|
||||
// We are using the simple integer indexing for the aquifers
|
||||
for (int i = 0; i < numAquifers(); ++i)
|
||||
{
|
||||
std::cout << "Aquifer[" << i << "]"
|
||||
<< " : Tc = " << aquifers()[i].time_constant()
|
||||
<< ", beta = " << aquifers()[i].aquifer_influx_constant() << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
// Get the last report step
|
||||
template<typename TypeTag>
|
||||
const SimulatorReport&
|
||||
BlackoilAquiferModel<TypeTag>:: lastReport() const
|
||||
{
|
||||
for (auto i = aquifers_.begin(); i != aquifers_.end(); ++i){
|
||||
(*i).print_private_members();
|
||||
}
|
||||
return last_report_;
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>::
|
||||
assemble( const SimulatorTimerInterface& timer,
|
||||
const int iterationIdx )
|
||||
{
|
||||
last_report_ = SimulatorReport();
|
||||
|
||||
// We need to update the reservoir pressures connected to the aquifer
|
||||
updateConnectionIntensiveQuantities();
|
||||
|
||||
if (iterationIdx == 0) {
|
||||
// We can do the Table check and coefficients update in this function
|
||||
// For now, it does nothing!
|
||||
prepareTimeStep(timer);
|
||||
}
|
||||
|
||||
if (iterationIdx == 0) {
|
||||
calculateExplicitQuantities();
|
||||
}
|
||||
|
||||
if (param_.solve_aquifereq_initially_ && iterationIdx == 0) {
|
||||
// solve the aquifer equations as a pre-processing step
|
||||
last_report_ = solveAquiferEq(timer);
|
||||
}
|
||||
|
||||
assembleAquiferEq(timer);
|
||||
|
||||
last_report_.converged = true;
|
||||
}
|
||||
|
||||
// Protected function: Update the primary variables
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: updatePrimaryVariables()
|
||||
{
|
||||
// Right now it doesn't do shit.
|
||||
}
|
||||
|
||||
// Protected function: Init the primary variables
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: initPrimaryVariablesEvaluation() const
|
||||
{
|
||||
// Right now it doesn't do shit.
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: updateConnectionIntensiveQuantities() const
|
||||
{
|
||||
ElementContext elemCtx(ebosSimulator_);
|
||||
const auto& gridView = ebosSimulator_.gridView();
|
||||
const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
|
||||
for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
|
||||
elemIt != elemEndIt;
|
||||
++elemIt)
|
||||
{
|
||||
elemCtx.updatePrimaryStencil(*elemIt);
|
||||
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: calculateExplicitQuantities()
|
||||
{
|
||||
// for (auto aqui = aquifers_.begin(); aqui!= aquifers_.end(); ++aqui)
|
||||
// {
|
||||
// std::cout << "calculateExplicitQuantities: Aquifer id = " << aqui->aquiferID() << std::endl;
|
||||
// aqui->calculateExplicitQuantities(ebosSimulator_);
|
||||
// }
|
||||
}
|
||||
|
||||
|
||||
template<typename TypeTag>
|
||||
SimulatorReport
|
||||
BlackoilAquiferModel<TypeTag>:: solveAquiferEq(const SimulatorTimerInterface& timer)
|
||||
{
|
||||
// We need to solve the equilibrium equation first to
|
||||
// obtain the initial pressure of water in the aquifer
|
||||
SimulatorReport report = SimulatorReport();
|
||||
return report;
|
||||
}
|
||||
|
||||
// Protected function: Return number of components in the model.
|
||||
template<typename TypeTag>
|
||||
int
|
||||
BlackoilAquiferModel<TypeTag>:: numComponents() const
|
||||
{
|
||||
if (numPhases() == 2) {
|
||||
return 2;
|
||||
}
|
||||
int numComp = FluidSystem::numComponents;
|
||||
if (has_solvent_) {
|
||||
numComp ++;
|
||||
}
|
||||
|
||||
return numComp;
|
||||
}
|
||||
|
||||
// Protected function: Return number of aquifers in the model.
|
||||
template<typename TypeTag>
|
||||
int
|
||||
BlackoilAquiferModel<TypeTag>:: numAquifers() const
|
||||
{
|
||||
return aquifers_.size();
|
||||
}
|
||||
|
||||
// Protected function: Return number of phases in the model.
|
||||
template<typename TypeTag>
|
||||
int
|
||||
BlackoilAquiferModel<TypeTag>:: numPhases() const
|
||||
{
|
||||
// Not implemented yet!!!!!!!!!!!!
|
||||
const auto& pu = phase_usage_;
|
||||
return pu.num_phases;
|
||||
}
|
||||
|
||||
|
||||
// Protected function: returns the phase index in ebos
|
||||
template<typename TypeTag>
|
||||
int
|
||||
BlackoilAquiferModel<TypeTag>:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
|
||||
{
|
||||
const auto& pu = phase_usage_;
|
||||
if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
|
||||
return FluidSystem::waterPhaseIdx;
|
||||
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
|
||||
return FluidSystem::oilPhaseIdx;
|
||||
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
|
||||
return FluidSystem::gasPhaseIdx;
|
||||
|
||||
assert(phaseIdx < 3);
|
||||
// for other phases return the index
|
||||
return phaseIdx;
|
||||
}
|
||||
|
||||
// Protected function which calls the individual aquifer models
|
||||
template<typename TypeTag>
|
||||
void
|
||||
BlackoilAquiferModel<TypeTag>:: assembleAquiferEq(const SimulatorTimerInterface& timer)
|
||||
{
|
||||
for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
|
||||
{
|
||||
std::cout << "assembleAquiferEq: Aquifer id = " << aquifer->aquiferID() << std::endl;
|
||||
aquifer->assembleAquiferEq(ebosSimulator_, timer);
|
||||
}
|
||||
}
|
||||
|
||||
// Protected function
|
||||
// some preparation work, mostly related to group control and RESV,
|
||||
// at the beginning of each time step (Not report step)
|
||||
template<typename TypeTag>
|
||||
void BlackoilAquiferModel<TypeTag>:: prepareTimeStep(const SimulatorTimerInterface& timer)
|
||||
{
|
||||
// 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);
|
||||
}
|
||||
}
|
||||
|
||||
// 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)//, std::vector< AquiferCarterTracy<TypeTag> >& aquifers)
|
||||
{
|
||||
const auto& deck = ebosSimulator.gridManager().deck();
|
||||
const auto& eclState = ebosSimulator.gridManager().eclState();
|
||||
|
||||
// Get all the carter tracy aquifer properties data and put it in aquifers vector
|
||||
AquiferCT aquiferct = AquiferCT(eclState,deck);
|
||||
|
||||
std::vector<AquiferCT::AQUCT_data> aquifersData = aquiferct.getAquifers();
|
||||
std::vector<AquiferCT::AQUANCON_data> aquanconData = aquiferct.getAquancon();
|
||||
|
||||
// for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData)
|
||||
// {
|
||||
|
||||
// }
|
||||
|
||||
auto ita = aquifersData.cbegin();
|
||||
auto f_lambda = [&] (AquiferCT::AQUANCON_data i) {
|
||||
aquifers.push_back( AquiferCarterTracy<TypeTag> (*ita++, i, numComponents(), gravity_ ) );
|
||||
};
|
||||
std::for_each( aquanconData.cbegin(), aquanconData.cend(), f_lambda );
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
@ -30,6 +30,9 @@
|
||||
#include <opm/autodiff/BlackoilModelParameters.hpp>
|
||||
#include <opm/autodiff/BlackoilWellModel.hpp>
|
||||
#include <opm/autodiff/WellConnectionAuxiliaryModule.hpp>
|
||||
#include <opm/autodiff/BlackoilAquiferModel.hpp>
|
||||
#include <opm/autodiff/GridHelpers.hpp>
|
||||
#include <opm/autodiff/GeoProps.hpp>
|
||||
#include <opm/autodiff/BlackoilDetails.hpp>
|
||||
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
|
||||
|
||||
@ -143,6 +146,7 @@ namespace Opm {
|
||||
BlackoilModelEbos(Simulator& ebosSimulator,
|
||||
const ModelParameters& param,
|
||||
BlackoilWellModel<TypeTag>& well_model,
|
||||
BlackoilAquiferModel<TypeTag>& aquifer_model,
|
||||
const NewtonIterationBlackoilInterface& linsolver,
|
||||
const bool terminal_output
|
||||
)
|
||||
@ -157,6 +161,7 @@ namespace Opm {
|
||||
, has_energy_(GET_PROP_VALUE(TypeTag, EnableEnergy))
|
||||
, param_( param )
|
||||
, well_model_ (well_model)
|
||||
, aquifer_model_(aquifer_model)
|
||||
, terminal_output_ (terminal_output)
|
||||
, current_relaxation_(1.0)
|
||||
, dx_old_(UgGridHelpers::numCells(grid_))
|
||||
@ -204,6 +209,7 @@ namespace Opm {
|
||||
wasSwitched_.resize(numDof);
|
||||
std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
|
||||
|
||||
aquiferModel().beginTimeStep();
|
||||
wellModel().beginTimeStep();
|
||||
|
||||
if (param_.update_equations_scaling_) {
|
||||
@ -349,6 +355,7 @@ namespace Opm {
|
||||
DUNE_UNUSED_PARAMETER(well_state);
|
||||
|
||||
wellModel().timeStepSucceeded();
|
||||
aquiferModel().timeStepSucceeded();
|
||||
ebosSimulator_.problem().endTimeStep();
|
||||
|
||||
}
|
||||
@ -366,9 +373,22 @@ namespace Opm {
|
||||
ebosSimulator_.model().linearizer().linearize();
|
||||
ebosSimulator_.problem().endIteration();
|
||||
|
||||
// -------- Well equations ----------
|
||||
// -------- Well and aquifer common variables ----------
|
||||
double dt = timer.currentStepLength();
|
||||
|
||||
// -------- Aquifer models ----------
|
||||
try
|
||||
{
|
||||
// Modify the Jacobian and residuals according to the aquifer models
|
||||
aquiferModel().assemble(timer, iterationIdx);
|
||||
}
|
||||
catch( const Dune::FMatrixError& e )
|
||||
{
|
||||
OPM_THROW(Opm::NumericalProblem,"Error when assembling aquifer models");
|
||||
}
|
||||
|
||||
// -------- Well equations ----------
|
||||
|
||||
try
|
||||
{
|
||||
// assembles the well equations and applies the wells to
|
||||
@ -1081,6 +1101,9 @@ namespace Opm {
|
||||
// Well Model
|
||||
BlackoilWellModel<TypeTag>& well_model_;
|
||||
|
||||
// Aquifer Model
|
||||
BlackoilAquiferModel<TypeTag>& aquifer_model_;
|
||||
|
||||
/// \brief Whether we print something to std::cout
|
||||
bool terminal_output_;
|
||||
/// \brief The number of cells of the global grid.
|
||||
@ -1100,6 +1123,41 @@ namespace Opm {
|
||||
const BlackoilWellModel<TypeTag>&
|
||||
wellModel() const { return well_model_; }
|
||||
|
||||
BlackoilAquiferModel<TypeTag>&
|
||||
aquiferModel() { return aquifer_model_; }
|
||||
|
||||
const BlackoilAquiferModel<TypeTag>&
|
||||
aquiferModel() const { return aquifer_model_; }
|
||||
|
||||
int flowPhaseToEbosCompIdx( const int phaseIdx ) const
|
||||
{
|
||||
const auto& pu = phaseUsage_;
|
||||
if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
|
||||
return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
|
||||
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
|
||||
return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
|
||||
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
|
||||
return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
|
||||
|
||||
// for other phases return the index
|
||||
return phaseIdx;
|
||||
}
|
||||
|
||||
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
|
||||
{
|
||||
const auto& pu = phaseUsage_;
|
||||
if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
|
||||
return FluidSystem::waterPhaseIdx;
|
||||
if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
|
||||
return FluidSystem::oilPhaseIdx;
|
||||
if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
|
||||
return FluidSystem::gasPhaseIdx;
|
||||
|
||||
assert(phaseIdx < 3);
|
||||
// for other phases return the index
|
||||
return phaseIdx;
|
||||
}
|
||||
|
||||
void beginReportStep()
|
||||
{
|
||||
ebosSimulator_.problem().beginEpisode();
|
||||
|
||||
@ -63,6 +63,7 @@ namespace Opm
|
||||
param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day );
|
||||
max_strict_iter_ = param.getDefault("max_strict_iter",8);
|
||||
solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_);
|
||||
solve_aquifereq_initially_ = param.getDefault("solve_aquifereq_initially",solve_aquifereq_initially_);
|
||||
update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
|
||||
use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_);
|
||||
deck_file_name_ = param.template get<std::string>("deck_filename");
|
||||
@ -94,6 +95,7 @@ namespace Opm
|
||||
max_inner_iter_ms_wells_ = 10;
|
||||
maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day );
|
||||
solve_welleq_initially_ = true;
|
||||
solve_aquifereq_initially_ = true;
|
||||
update_equations_scaling_ = false;
|
||||
use_update_stabilization_ = true;
|
||||
use_multisegment_well_ = false;
|
||||
|
||||
@ -77,6 +77,9 @@ namespace Opm
|
||||
/// Solve well equation initially
|
||||
bool solve_welleq_initially_;
|
||||
|
||||
/// Solve aquifer equation initially
|
||||
bool solve_aquifereq_initially_;
|
||||
|
||||
/// Update scaling factors for mass balance equations
|
||||
bool update_equations_scaling_;
|
||||
|
||||
|
||||
@ -86,6 +86,7 @@ namespace Opm {
|
||||
{
|
||||
Base::prepareStep(timer, reservoir_state, well_state);
|
||||
Base::param_.solve_welleq_initially_ = false;
|
||||
Base::param_.solve_aquifereq_initially_ = false;
|
||||
SolutionState state0 = variableState(reservoir_state, well_state);
|
||||
asImpl().makeConstantState(state0);
|
||||
asImpl().computeAccum(state0, 0);
|
||||
|
||||
@ -29,6 +29,9 @@
|
||||
#include <opm/autodiff/BlackoilModelParameters.hpp>
|
||||
#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>
|
||||
@ -65,6 +68,7 @@ public:
|
||||
typedef BlackoilModelParameters ModelParameters;
|
||||
typedef NonlinearSolver<Model> Solver;
|
||||
typedef BlackoilWellModel<TypeTag> WellModel;
|
||||
typedef BlackoilAquiferModel<TypeTag> AquiferModel;
|
||||
|
||||
|
||||
/// Initialise from parameters and objects to observe.
|
||||
@ -186,6 +190,8 @@ public:
|
||||
auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(schedule(), grid());
|
||||
ebosSimulator_.model().addAuxiliaryModule(auxMod);
|
||||
}
|
||||
AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_);
|
||||
// aquifer_model.hack_init(ebosSimulator_);
|
||||
|
||||
// Main simulation loop.
|
||||
while (!timer.done()) {
|
||||
@ -202,7 +208,9 @@ public:
|
||||
|
||||
well_model.beginReportStep(timer.currentStepNum());
|
||||
|
||||
auto solver = createSolver(well_model);
|
||||
aquifer_model.beginReportStep(timer.currentStepNum());
|
||||
|
||||
auto solver = createSolver(well_model, aquifer_model);
|
||||
|
||||
// write the inital state at the report stage
|
||||
if (timer.initialStep()) {
|
||||
@ -266,6 +274,7 @@ public:
|
||||
}
|
||||
|
||||
solver->model().endReportStep();
|
||||
aquifer_model.endReportStep();
|
||||
well_model.endReportStep();
|
||||
|
||||
// take time that was used to solve system for this reportStep
|
||||
@ -308,6 +317,9 @@ public:
|
||||
total_timer.stop();
|
||||
report.total_time = total_timer.secsSinceStart();
|
||||
report.converged = true;
|
||||
|
||||
auto reportaquifer = aquifer_model.lastReport();
|
||||
|
||||
return report;
|
||||
}
|
||||
|
||||
@ -320,11 +332,12 @@ public:
|
||||
|
||||
protected:
|
||||
|
||||
std::unique_ptr<Solver> createSolver(WellModel& well_model)
|
||||
std::unique_ptr<Solver> createSolver(WellModel& well_model, AquiferModel& aquifer_model)
|
||||
{
|
||||
auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
|
||||
model_param_,
|
||||
well_model,
|
||||
aquifer_model,
|
||||
solver_,
|
||||
terminal_output_));
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user