diff --git a/.gitignore b/.gitignore
index 585d691bc..3f98f891f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,3 +32,5 @@ test_vec
# emacs directory setting:
.dir-locals.el
+
+build
diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake
old mode 100644
new mode 100755
index b01e3c3fd..c4175c03c
--- a/compareECLFiles.cmake
+++ b/compareECLFiles.cmake
@@ -195,6 +195,13 @@ add_test_compareECLFiles(CASENAME spe1_thermal
REL_TOL ${rel_tol}
DIR spe1)
+add_test_compareECLFiles(CASENAME ctaquifer_2d_oilwater
+ FILENAME 2D_OW_CTAQUIFER
+ SIMULATOR flow
+ ABS_TOL ${abs_tol}
+ REL_TOL ${rel_tol}
+ DIR aquifer-oilwater)
+
foreach(SIM flow flow_legacy)
add_test_compareECLFiles(CASENAME spe3
FILENAME SPE3CASE1
diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp
new file mode 100644
index 000000000..ff9841016
--- /dev/null
+++ b/opm/autodiff/AquiferCarterTracy.hpp
@@ -0,0 +1,366 @@
+/*
+ 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 .
+*/
+
+#ifndef OPM_AQUIFERCT_HEADER_INCLUDED
+#define OPM_AQUIFERCT_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+
+namespace Opm
+{
+
+ template
+ class AquiferCarterTracy
+ {
+
+ public:
+
+ 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;
+
+ static const int numEq = BlackoilIndices::numEq;
+ typedef double Scalar;
+
+ typedef DenseAd::Evaluation Eval;
+ typedef Opm::BlackOilFluidState FluidState;
+
+ static const auto waterCompIdx = FluidSystem::waterCompIdx;
+ static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
+
+
+
+ AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data,
+ const Aquancon::AquanconOutput& connection,
+ Simulator& ebosSimulator )
+ : ebos_simulator_ (ebosSimulator),
+ aquct_data_ (aquct_data),
+ gravity_ (ebos_simulator_.problem().gravity()[2])
+ {
+ initQuantities(connection);
+ }
+
+ inline void assembleAquiferEq(const SimulatorTimerInterface& timer)
+ {
+ auto& ebosJac = ebos_simulator_.model().linearizer().matrix();
+ auto& ebosResid = ebos_simulator_.model().linearizer().residual();
+
+ 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 = *(ebos_simulator_.model().cachedIntensiveQuantities(cellID, /*timeIdx=*/ 0));
+ // This is the pressure at td + dt
+ updateCellPressure(pressure_current_,idx,intQuants);
+ updateCellDensity(idx,intQuants);
+ calculateInflowRate(idx, timer);
+ qinflow = Qai_.at(idx);
+ 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);
+ }
+ }
+ }
+
+ 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 = *(ebos_simulator_.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0));
+ updateCellPressure(pressure_previous_ ,idx,intQuants);
+ }
+ }
+
+ inline void afterTimeStep(const SimulatorTimerInterface& timer)
+ {
+ for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai)
+ {
+ W_flux_ += (*Qai)*timer.currentStepLength();
+ }
+ }
+
+
+ private:
+ Simulator& ebos_simulator_;
+
+ // Grid variables
+ std::vector cell_idx_;
+ std::vector faceArea_connected_;
+
+ // Quantities at each grid id
+ std::vector cell_depth_;
+ std::vector pressure_previous_;
+ std::vector pressure_current_;
+ std::vector Qai_;
+ std::vector rhow_;
+ std::vector alphai_;
+
+ // Variables constants
+ const AquiferCT::AQUCT_data aquct_data_;
+
+ Scalar mu_w_ , //water viscosity
+ beta_ , // Influx constant
+ Tc_ , // Time constant
+ pa0_ , // initial aquifer pressure
+ gravity_ ; // gravitational acceleration
+
+ Eval W_flux_;
+
+
+ inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
+ {
+ // We use the opm-common numeric linear interpolator
+ pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td);
+ pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td);
+ }
+
+ 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
+ initializeConnections(connection);
+
+ 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 updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants)
+ {
+ const auto& fs = intQuants.fluidState();
+ pressure_water.at(idx) = fs.pressure(waterPhaseIdx);
+ }
+
+ inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants)
+ {
+ const auto& fs = intQuants.fluidState();
+ 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) - aquct_data_.d0) - pressure_previous_.at(idx);
+ return dp;
+ }
+
+ // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription
+ inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer)
+ {
+ const Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc_;
+ const Scalar td = timer.simulationTimeElapsed() / Tc_;
+ Scalar PItdprime = 0.;
+ Scalar PItd = 0.;
+ 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 calculateInflowRate(int idx, const SimulatorTimerInterface& timer)
+ {
+ Scalar a, b;
+ calculateEqnConstants(a,b,idx,timer);
+ Qai_.at(idx) = alphai_.at(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx) ) );
+ }
+
+ inline void calculateAquiferConstants()
+ {
+ // 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 initializeConnections(const Aquancon::AquanconOutput& connection)
+ {
+ const auto& eclState = ebos_simulator_.vanguard().eclState();
+ const auto& ugrid = ebos_simulator_.vanguard().grid();
+ const auto& grid = eclState.getInputGrid();
+
+ cell_idx_ = connection.global_index;
+ auto globalCellIdx = ugrid.globalCell();
+
+ assert( cell_idx_ == connection.global_index);
+ assert( (cell_idx_.size() == connection.influx_coeff.size()) );
+ assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) );
+ 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(), aquct_data_.d0);
+ alphai_.resize(cell_idx_.size(), 1.0);
+ faceArea_connected_.resize(cell_idx_.size(),0.0);
+ Scalar faceArea;
+
+ auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
+ auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid);
+
+ // Translate the C face tag into the enum used by opm-parser's TransMult class
+ Opm::FaceDir::DirEnum faceDirection;
+
+ // denom_face_areas is the sum of the areas connected to an aquifer
+ Scalar denom_face_areas = 0.;
+ for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
+ {
+ auto cellFacesRange = cell2Faces[cell_idx_.at(idx)];
+
+ for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter)
+ {
+ // The index of the face in the compressed grid
+ const int faceIdx = *cellFaceIter;
+
+ // the logically-Cartesian direction of the face
+ const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
+
+ 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))
+ {
+ // Check now if the face is outside of the reservoir, or if it adjoins an inactive cell
+ // Do not make the connection if the product of the two cellIdx > 0. This is because the
+ // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining)
+ faceArea = (faceCells(faceIdx,0)*faceCells(faceIdx,1) > 0)? 0. : Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
+ faceArea_connected_.at(idx) = faceArea;
+ denom_face_areas += ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) );
+ }
+ }
+ auto cellCenter = grid.getCellCenter(cell_idx_.at(idx));
+ cell_depth_.at(idx) = cellCenter[2];
+ }
+
+ for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
+ {
+ alphai_.at(idx) = ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas;
+ }
+ }
+
+ inline void calculateAquiferCondition()
+ {
+
+ int pvttableIdx = aquct_data_.pvttableID - 1;
+
+ rhow_.resize(cell_idx_.size(),0.);
+
+ if (aquct_data_.p0 < 1.0)
+ {
+ pa0_ = calculateReservoirEquilibrium();
+ }
+ else
+ {
+ pa0_ = aquct_data_.p0;
+ }
+
+ // Initialize a FluidState object first
+ FluidState fs_aquifer;
+ // We use the temperature of the first cell connected to the aquifer
+ // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs
+ fs_aquifer.assign( ebos_simulator_.model().cachedIntensiveQuantities(cell_idx_.at(0), /*timeIdx=*/ 0)->fluidState() );
+ Eval temperature_aq, pa0_mean;
+ temperature_aq = fs_aquifer.temperature(0);
+ pa0_mean = pa0_;
+
+ Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
+
+ mu_w_ = mu_w_aquifer.value();
+
+ }
+
+ // This function is for calculating the aquifer properties from equilibrium state with the reservoir
+ inline Scalar calculateReservoirEquilibrium()
+ {
+ // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
+ std::vector pw_aquifer;
+ Scalar water_pressure_reservoir;
+
+ for (size_t idx = 0; idx < cell_idx_.size(); ++idx)
+ {
+ size_t cellIDx = cell_idx_.at(idx);
+ const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0));
+ const auto& fs = intQuants.fluidState();
+
+ 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.
+ Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.)/pw_aquifer.size();
+ return aquifer_pres_avg;
+ }
+
+
+ }; // class AquiferCarterTracy
+
+
+} // namespace Opm
+
+#endif
\ No newline at end of file
diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp
new file mode 100644
index 000000000..f904ec7bf
--- /dev/null
+++ b/opm/autodiff/BlackoilAquiferModel.hpp
@@ -0,0 +1,82 @@
+/*
+ File adapted from BlackoilWellModel.hpp
+
+ 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 .
+*/
+
+
+#ifndef OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED
+#define OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+
+namespace Opm {
+
+ /// Class for handling the blackoil well model.
+ template
+ class BlackoilAquiferModel {
+
+ public:
+
+ // --------- Types ---------
+ typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+ typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ typedef AquiferCarterTracy Aquifer_object;
+
+ explicit BlackoilAquiferModel(Simulator& ebosSimulator);
+
+ // 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 end of a time step
+ void timeStepSucceeded(const SimulatorTimerInterface& timer);
+
+ protected:
+
+ Simulator& ebosSimulator_;
+
+ std::vector aquifers_;
+
+ // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
+ void init();
+
+ void updateConnectionIntensiveQuantities() const;
+
+ void assembleAquiferEq(const SimulatorTimerInterface& timer);
+
+ // at the beginning of each time step (Not report step)
+ void prepareTimeStep(const SimulatorTimerInterface& timer);
+
+ bool aquiferActive() const;
+
+ };
+
+
+} // namespace Opm
+
+#include "BlackoilAquiferModel_impl.hpp"
+#endif
\ No newline at end of file
diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp
new file mode 100644
index 000000000..e92c952cc
--- /dev/null
+++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp
@@ -0,0 +1,130 @@
+namespace Opm {
+
+
+ template
+ BlackoilAquiferModel::
+ BlackoilAquiferModel(Simulator& ebosSimulator)
+ : ebosSimulator_(ebosSimulator)
+ {
+ init();
+ }
+
+
+ // called at the end of a time step
+ template
+ void
+ BlackoilAquiferModel:: timeStepSucceeded(const SimulatorTimerInterface& timer)
+ {
+ if ( !aquiferActive() ) {
+ return;
+ }
+
+ for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
+ {
+ aquifer->afterTimeStep(timer);
+ }
+ }
+
+ template
+ void
+ BlackoilAquiferModel::
+ assemble( const SimulatorTimerInterface& timer,
+ const int iterationIdx )
+ {
+ if ( !aquiferActive() ) {
+ return;
+ }
+
+ // 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);
+ }
+
+ assembleAquiferEq(timer);
+ }
+
+
+ template
+ void
+ BlackoilAquiferModel:: updateConnectionIntensiveQuantities() const
+ {
+ ElementContext elemCtx(ebosSimulator_);
+ const auto& gridView = ebosSimulator_.gridView();
+ const auto& elemEndIt = gridView.template end();
+ for (auto elemIt = gridView.template begin();
+ elemIt != elemEndIt;
+ ++elemIt)
+ {
+ elemCtx.updatePrimaryStencil(*elemIt);
+ elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+ }
+ }
+
+ // Protected function which calls the individual aquifer models
+ template
+ void
+ BlackoilAquiferModel:: assembleAquiferEq(const SimulatorTimerInterface& timer)
+ {
+ for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
+ {
+ aquifer->assembleAquiferEq(timer);
+ }
+ }
+
+ // Protected function
+ // some preparation work, mostly related to group control and RESV,
+ // at the beginning of each time step (Not report step)
+ template
+ void BlackoilAquiferModel:: 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->beforeTimeStep(timer);
+ }
+ }
+
+ // Initialize the aquifers in the deck
+ template
+ void
+ BlackoilAquiferModel:: init()
+ {
+ const auto& deck = ebosSimulator_.vanguard().deck();
+
+ if ( !deck.hasKeyword("AQUCT") ) {
+ return ;
+ }
+
+ updateConnectionIntensiveQuantities();
+ const auto& eclState = ebosSimulator_.vanguard().eclState();
+
+ // Get all the carter tracy aquifer properties data and put it in aquifers vector
+ const AquiferCT aquiferct = AquiferCT(eclState,deck);
+ const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck);
+
+ std::vector aquifersData = aquiferct.getAquifers();
+ std::vector aquifer_connection = aquifer_connect.getAquOutput();
+
+ assert( aquifersData.size() == aquifer_connection.size() );
+
+
+ for (size_t i = 0; i < aquifersData.size(); ++i)
+ {
+ aquifers_.push_back(
+ AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), ebosSimulator_)
+ );
+ }
+ }
+
+ template
+ bool
+ BlackoilAquiferModel:: aquiferActive() const
+ {
+ return !aquifers_.empty();
+ }
+
+} // namespace Opm
diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp
index c16461496..9f37945d0 100644
--- a/opm/autodiff/BlackoilModelEbos.hpp
+++ b/opm/autodiff/BlackoilModelEbos.hpp
@@ -29,6 +29,7 @@
#include
#include
+#include
#include
#include
#include
@@ -143,6 +144,7 @@ namespace Opm {
BlackoilModelEbos(Simulator& ebosSimulator,
const ModelParameters& param,
BlackoilWellModel& well_model,
+ BlackoilAquiferModel& aquifer_model,
const NewtonIterationBlackoilInterface& linsolver,
const bool terminal_output
)
@@ -157,6 +159,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_))
@@ -349,6 +352,7 @@ namespace Opm {
DUNE_UNUSED_PARAMETER(well_state);
wellModel().timeStepSucceeded();
+ aquiferModel().timeStepSucceeded(timer);
ebosSimulator_.problem().endTimeStep();
}
@@ -365,9 +369,22 @@ namespace Opm {
ebosSimulator_.problem().beginIteration();
ebosSimulator_.model().linearizer().linearize();
ebosSimulator_.problem().endIteration();
+
+ // -------- Aquifer models ----------
+ try
+ {
+ // Modify the Jacobian and residuals according to the aquifer models
+ aquiferModel().assemble(timer, iterationIdx);
+ }
+ catch( ... )
+ {
+ OPM_THROW(Opm::NumericalIssue,"Error when assembling aquifer models");
+ }
+
+ // -------- Current time step length ----------
+ const double dt = timer.currentStepLength();
// -------- Well equations ----------
- double dt = timer.currentStepLength();
try
{
@@ -409,13 +426,13 @@ namespace Opm {
if (elem.partitionType() != Dune::InteriorEntity)
continue;
- unsigned globalElemIdx = elemMapper.index(elem);
+ unsigned globalElemIdx = elemMapper.index(elem);
const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
Scalar pressureNew;
- pressureNew = priVarsNew[Indices::pressureSwitchIdx];
+ pressureNew = priVarsNew[Indices::pressureSwitchIdx];
- Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
+ Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 };
Scalar oilSaturationNew = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
@@ -458,7 +475,7 @@ namespace Opm {
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
- resultDelta += tmp*tmp;
+ resultDelta += tmp*tmp;
resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
}
}
@@ -466,9 +483,9 @@ namespace Opm {
resultDelta = gridView.comm().sum(resultDelta);
resultDenom = gridView.comm().sum(resultDenom);
- if (resultDenom > 0.0)
- return resultDelta/resultDenom;
- return 0.0;
+ if (resultDenom > 0.0)
+ return resultDelta/resultDenom;
+ return 0.0;
}
@@ -1081,6 +1098,9 @@ namespace Opm {
// Well Model
BlackoilWellModel& well_model_;
+ // Aquifer Model
+ BlackoilAquiferModel& aquifer_model_;
+
/// \brief Whether we print something to std::cout
bool terminal_output_;
/// \brief The number of cells of the global grid.
@@ -1100,6 +1120,9 @@ namespace Opm {
const BlackoilWellModel&
wellModel() const { return well_model_; }
+ BlackoilAquiferModel&
+ aquiferModel() { return aquifer_model_; }
+
void beginReportStep()
{
ebosSimulator_.problem().beginEpisode();
@@ -1112,7 +1135,6 @@ namespace Opm {
private:
-
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }
double drMaxRel() const { return param_.dr_max_rel_; }
@@ -1123,4 +1145,4 @@ namespace Opm {
};
} // namespace Opm
-#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
+#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
\ No newline at end of file
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
index 4445eb390..674f14918 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
@@ -29,6 +29,7 @@
#include
#include
#include
+#include
#include
#include
#include
@@ -65,6 +66,7 @@ public:
typedef BlackoilModelParameters ModelParameters;
typedef NonlinearSolver Solver;
typedef BlackoilWellModel WellModel;
+ typedef BlackoilAquiferModel AquiferModel;
/// Initialise from parameters and objects to observe.
@@ -187,6 +189,8 @@ public:
ebosSimulator_.model().addAuxiliaryModule(auxMod);
}
+ AquiferModel aquifer_model(ebosSimulator_);
+
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
@@ -202,7 +206,7 @@ public:
well_model.beginReportStep(timer.currentStepNum());
- auto solver = createSolver(well_model);
+ auto solver = createSolver(well_model, aquifer_model);
// write the inital state at the report stage
if (timer.initialStep()) {
@@ -308,6 +312,7 @@ public:
total_timer.stop();
report.total_time = total_timer.secsSinceStart();
report.converged = true;
+
return report;
}
@@ -320,11 +325,12 @@ public:
protected:
- std::unique_ptr createSolver(WellModel& well_model)
+ std::unique_ptr createSolver(WellModel& well_model, AquiferModel& aquifer_model)
{
auto model = std::unique_ptr(new Model(ebosSimulator_,
model_param_,
well_model,
+ aquifer_model,
solver_,
terminal_output_));
diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh
index ab88a21c9..8d384e50d 100755
--- a/tests/update_reference_data.sh
+++ b/tests/update_reference_data.sh
@@ -20,7 +20,7 @@ copyToReferenceDir () {
}
tests=${@:2}
-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"
+test -z "$tests" && tests="spe11 spe12 spe12p spe1oilgas spe1nowells spe1thermal ctaquifer_2d_oilwater 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-tests/norne/flow
@@ -98,6 +98,15 @@ for test_name in ${tests}; do
EGRID INIT SMSPEC UNRST UNSMRY
fi
+ if grep -q "ctaquifer_2d_oilwater" <<< $test_name
+ then
+ copyToReferenceDir \
+ $configuration/build-opm-simulators/tests/results/flow+ctaquifer_2d_oilwater/ \
+ $OPM_TESTS_ROOT/aquifer-oilwater/opm-simulation-reference/flow \
+ 2D_OW_CTAQUIFER \
+ EGRID INIT SMSPEC UNRST UNSMRY
+ fi
+
if grep -q "msw_2d_h" <<< $test_name
then
copyToReferenceDir \