Implemented initialization procedure.

opm/autodiff/AQUCT_params.hpp

@@ -1,31 +0,0 @@
-#ifndef OPM_AQUCT_HEADER_INCLUDED
-#define OPM_AQUCT_HEADER_INCLUDED
-
-
-    struct AQUCT_params{
-        
-            // Aquifer ID
-            int aquiferID;
-            // Table IDs
-            int inftableID, pvttableID;
-            // Perforation cell id
-            int cell_id;
-            // Variables constants
-            double  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).
-                    std::vector<double> td, pi;
-    
-
-    };
-
-
-
-
-#endif            

opm/autodiff/AquiferCarterTracy.hpp

@@ -23,6 +23,7 @@
 
 #include <Eigen/QR>
 #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/core/props/BlackoilPhases.hpp>
@@ -30,10 +31,12 @@
 
 #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>
 
@@ -65,44 +68,30 @@ namespace Opm
             typedef double Scalar;
 
             typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
+            typedef Opm::BlackOilFluidState<Eval, FluidSystem> FluidState;
+            typedef typename FluidSystem::WaterPvt WaterPvt;
 
             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)
+            explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection,
-            : phi_aq_ (1.0), //
+                                         const int numComponents, const Scalar gravity, const Simulator& ebosSimulator    )
-              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), //
+              p0_defaulted_ (params.p0_defaulted),
+              pa0_ (params.p0),
               theta_ (params.theta), //
               c2_ (params.c2), //
               d0_ (params.d0),
@@ -111,12 +100,11 @@ namespace Opm
               aquiferID_ (params.aquiferID),
               inftableID_ (params.inftableID),
               pvttableID_ (params.pvttableID),
-              cell_idx_ (params.cell_id),
               num_components_ (numComponents),
-              gravity_ (gravity)
+              gravity_ (gravity),
+              ebos_simulator_ (ebosSimulator)
             {
-                mu_w_ = 1e-3;
+                init_quantities(connection);
-                init_quantities(aquanconParams);
             }
 
             inline const PhaseUsage&
@@ -127,79 +115,35 @@ namespace Opm
                 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));
+                    const IntensiveQuantities 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];
+                    qinflow = Qai_.at(idx);
-                    ebosResid[*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)] -= qinflow.value();
+                    ebosResid[*cellID][(FluidSystem::waterCompIdx)] -= 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);
+                        ebosJac[*cellID][*cellID][(FluidSystem::waterCompIdx)][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;
                 }
             }
 
@@ -214,50 +158,19 @@ namespace Opm
                 }
             }
 
-            inline void after_time_step()
+            inline void after_time_step(const SimulatorTimerInterface& timer)
             {
                 for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai)
                 {
-                    W_flux_ += (*Qai);
+                    W_flux_ += (*Qai)*timer.currentStepLength();
-                }
-                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 double area_fraction(const size_t i)
+            {
+                return alphai_.at(i);
+            }
+
             inline const std::vector<int> cell_id() const
             {
                 return cell_idx_;
@@ -269,9 +182,9 @@ namespace Opm
             }
 
 
-
+        private:
-        protected:
             const PhaseUsage* phase_usage_;
+            const Simulator& ebos_simulator_;
 
 
             // Aquifer ID, and other IDs
@@ -279,14 +192,17 @@ namespace Opm
             int num_components_;
 
             // Grid variables
-            std::vector<int> cell_idx_;
+
+            std::vector<size_t> cell_idx_;
+            std::vector<Scalar> faceArea_connected_;
 
             // Quantities at each grid id
             std::vector<Scalar> cell_depth_;
-            std::vector<Scalar> pressure_previous_;
+            std::vector<Eval> pressure_previous_;
-            std::vector<Scalar> pressure_current_;
+            std::vector<Eval> pressure_current_;
-            std::vector<Scalar> Qai_;
+            std::vector<Eval> Qai_;
-            std::vector<Scalar> rhow_;
+            std::vector<Eval> rhow_;
+            std::vector<Scalar> alphai_;
 
             // Variables constants
             Scalar mu_w_ , //water viscosity
@@ -304,8 +220,9 @@ namespace Opm
             std::vector<Scalar> aqutab_td_, aqutab_pi_;
 
             // Cumulative flux
-            Scalar W_flux_, dt_, pa0_, gravity_;
+            Scalar dt_, pa0_, gravity_;
-
+            bool p0_defaulted_;
+            Eval W_flux_;
             // Also return the polynomial fit
             std::vector<Scalar> coeff_;
             
@@ -316,7 +233,7 @@ namespace Opm
             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;
+                size_t 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;
@@ -337,60 +254,200 @@ namespace Opm
                     coeff[i] = result[i];
             }
 
-            inline void init_quantities(const AquiferCT::AQUANCON_data& aquanconParams)
+            inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td)
             {
+                // http://kluge.in-chemnitz.de/opensource/spline/
+            }
+
+            inline void init_quantities(const Aquancon::AquanconOutput& connection)
+            {
+                // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
                 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
+                // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function
-                pa0_ = 45e6;
+                initialize_connections(connection);
+
+                calculate_aquifer_condition();
 
                 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
+                Qai_.resize(cell_idx_.size(), 0.0);
-                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);
+                polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 1, true);
             }
 
-            inline void get_current_Pressure_cell(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
+            inline void get_current_Pressure_cell(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
             {
                 const auto& fs = intQuants.fluidState();
-                pressure_water[idx] = fs.pressure(FluidSystem::waterPhaseIdx).value();
+                pressure_water.at(idx) = fs.pressure(FluidSystem::waterPhaseIdx);
             }
 
-            inline void get_current_density_cell(std::vector<Scalar>& rho_water, const int idx, const IntensiveQuantities& intQuants)
+            inline void get_current_density_cell(std::vector<Eval>& rho_water, const int idx, const IntensiveQuantities& intQuants)
             {
                 const auto& fs = intQuants.fluidState();
-                rho_water[idx] = fs.density(FluidSystem::waterPhaseIdx).value();
+                rho_water.at(idx) = fs.density(FluidSystem::waterPhaseIdx);
             }
 
             inline Scalar dpai(int idx)
             {
-                Scalar dp = pa0_ - rhow_[idx]*gravity_*(cell_depth_[idx] - d0_) - pressure_previous_[idx];
+                Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_) - pressure_previous_.at(idx).value();
                 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)
             {
-                // 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;
+                Scalar PItdprime = coeff_.at(1);
-                a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_ * PItdprime) ) / ( PItd - td*PItdprime );
+                Scalar PItd = coeff_.at(0) + coeff_.at(1)*td_plus_dt;
-                b = beta / Tc / ( PItd - td*PItdprime);
+                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)
             {
                 Scalar a, b;
                 calculate_a_b_constants(a,b,idx,timer);
-                // This function implements Eq 5.7 of the EclipseTechnicalDescription
+                Qai_.at(idx) = area_fraction(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx).value() ) );
-                Qai_[idx] = area_fraction(idx)*( a - b * ( pressure_current_[idx] - pressure_previous_[idx] ) );
+            }
+
+            inline const Scalar time_constant() const
+            {
+                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;
+            }
+
+            // 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)
+            {
+                const auto& eclState = ebos_simulator_.vanguard().eclState();
+                const auto& ugrid = ebos_simulator_.vanguard().grid();
+                const auto& grid = eclState.getInputGrid();
+
+                cell_idx_ = connection.global_index;
+
+                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(), d0_);
+                alphai_.resize(cell_idx_.size(), 1.0);
+                faceArea_connected_.resize(cell_idx_.size(),0.0);
+
+                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);
+
+                        
+                        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;
+
+                        if (faceDirection == connection.reservoir_face_dir.at(idx))
+                        {
+                            faceArea_connected_.at(idx) = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
+                            denom_face_areas += 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) = faceArea_connected_.at(idx)/denom_face_areas;
+                }
+            }
+
+            inline void calculate_aquifer_condition()
+            {
+
+                int pvttableIdx = pvttableID_ - 1;
+                
+                rhow_.resize(cell_idx_.size(),0.);
+                
+                if (p0_defaulted_)
+                {
+                   pa0_ = calculate_reservoir_equilibrium(rhow_);
+                }
+
+                // 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_;
+
+                // rho_mean = FluidSystem::referenceDensity(waterPhaseIdx, pvttableIdx)
+                //             *FluidSystem::waterPvt().inverseFormationVolumeFactor(pvttableIdx, temperature_aq, pa0_mean);
+
+                Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean);
+                std::cout << "Pa0 = " << pa0_mean << ", viscosity = " << mu_w_aquifer.value() << std::endl;
+
+                mu_w_ = mu_w_aquifer.value();
+
+            }
+
+            // 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)
+            {
+                // 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;
+
+                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.push_back( fs.pressure(FluidSystem::waterPhaseIdx).value() );
+                    rho_water_reservoir.at(idx) = fs.density(FluidSystem::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) );
+                }
+
+                // 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;
             }
 
 

opm/autodiff/BlackoilAquiferModel.hpp

@@ -3,11 +3,6 @@
 
   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
-
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -40,6 +35,7 @@
 #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>
 
@@ -49,7 +45,6 @@
 #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>
@@ -86,8 +81,6 @@ namespace Opm {
             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,
@@ -102,7 +95,7 @@ namespace Opm {
             // called at the beginning of a time step
             void beginTimeStep();
             // called at the end of a time step
-            void timeStepSucceeded();
+            void timeStepSucceeded(const SimulatorTimerInterface& timer);
 
             // called at the beginning of a report step
             void beginReportStep(const int time_step);
@@ -117,7 +110,7 @@ namespace Opm {
                 return ebosSimulator_;
             }
 
-            /// Hack function to get what I need from parser
+            // 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:
@@ -141,17 +134,9 @@ namespace Opm {
 
             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;
 
@@ -159,8 +144,6 @@ namespace Opm {
 
             int numPhases() const;
 
-            int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
-
             void assembleAquiferEq(const SimulatorTimerInterface& timer);
 
             SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer);

opm/autodiff/BlackoilAquiferModel_impl.hpp

@@ -12,7 +12,7 @@ namespace Opm {
         , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
         , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
     {
-        const auto& eclState = ebosSimulator_.gridManager().eclState();
+        const auto& eclState = ebosSimulator_.vanguard().eclState();
         phase_usage_ = phaseUsageFromDeck(eclState);
 
         active_.resize(phase_usage_.MaxNumPhases, false);
@@ -38,17 +38,17 @@ namespace Opm {
     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()
+    BlackoilAquiferModel<TypeTag>:: timeStepSucceeded(const SimulatorTimerInterface& timer)
     {
         for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
         {
-            aquifer->after_time_step();
+            aquifer->after_time_step(timer);
         }
     }
 
@@ -57,7 +57,7 @@ namespace Opm {
     void
     BlackoilAquiferModel<TypeTag>:: beginReportStep(const int time_step) 
     {
-        // Right now it doesn't do shit.
+
     }
 
     // called at the end of a report step
@@ -65,14 +65,7 @@ namespace Opm {
     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
@@ -80,9 +73,6 @@ namespace Opm {
     const SimulatorReport& 
     BlackoilAquiferModel<TypeTag>:: lastReport() const 
     {
-        for (auto i = aquifers_.begin(); i != aquifers_.end(); ++i){
-            (*i).print_private_members();
-        }
         return last_report_;
     }
 
@@ -93,7 +83,6 @@ namespace Opm {
               const int iterationIdx                )
     {
         last_report_ = SimulatorReport();
-        
         // We need to update the reservoir pressures connected to the aquifer
         updateConnectionIntensiveQuantities();
 
@@ -103,35 +92,14 @@ namespace Opm {
             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
@@ -149,17 +117,6 @@ namespace Opm {
         }
     }
 
-    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 
@@ -200,30 +157,10 @@ namespace Opm {
     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
@@ -231,7 +168,6 @@ namespace Opm {
     {
         for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer)
         {
-            std::cout << "assembleAquiferEq: Aquifer id = " << aquifer->aquiferID() << std::endl;
             aquifer->assembleAquiferEq(ebosSimulator_, timer);
         }
     }
@@ -261,47 +197,27 @@ namespace Opm {
     // 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)
+    BlackoilAquiferModel<TypeTag>:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy<TypeTag> >& aquifers)
     {
-        const auto& deck = ebosSimulator.gridManager().deck();
+        updateConnectionIntensiveQuantities();
-        const auto& eclState = ebosSimulator.gridManager().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);
 
         std::vector<AquiferCT::AQUCT_data> aquifersData = aquiferct.getAquifers();
-        std::vector<AquiferCT::AQUANCON_data> aquanconData = aquiferct.getAquancon();
+        std::vector<Aquancon::AquanconOutput> aquifer_connection = aquifer_connect.getAquOutput();
 
-        // for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData)
+        assert( aquifersData.size() == aquifer_connect.size() );
-        // {
 
-        // }
 
-        auto ita = aquifersData.cbegin();
+        for (int i = 0; i < aquifersData.size(); ++i)
-        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 );
-    }
-
-    // Begin the hack to initialize the aquifers in the deck
-    template<typename TypeTag>
-    std::vector< AquiferCarterTracy<TypeTag> >
-    BlackoilAquiferModel<TypeTag>:: hack_init(const Simulator& ebosSimulator)//, std::vector< AquiferCarterTracy<TypeTag> >& aquifers)
-    {
-        std::vector< AquiferCarterTracy<TypeTag> > aquifers;
-        /** Begin hack!!!!! */
-        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();
-
-        for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData)
         {
-            aquifers.push_back( AquiferCarterTracy<TypeTag> (*aquiferData, numComponents(), gravity_ ) );
+            aquifers.push_back( 
+                                 AquiferCarterTracy<TypeTag> (aquifersData.at(i), aquifer_connection.at(i), numComponents(), gravity_, ebosSimulator_) 
+                              );
         }
     }
 

opm/autodiff/BlackoilModelEbos.hpp

@@ -29,10 +29,10 @@
 
 #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/WellConnectionAuxiliaryModule.hpp>
 #include <opm/autodiff/BlackoilDetails.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 
@@ -350,12 +350,12 @@ 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);
 
             wellModel().timeStepSucceeded();
-            aquiferModel().timeStepSucceeded();
+            aquiferModel().timeStepSucceeded(timer);
             ebosSimulator_.problem().endTimeStep();
 
         }
@@ -429,13 +429,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];
@@ -478,7 +478,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];
                 }
             }
@@ -486,9 +486,9 @@ namespace Opm {
             resultDelta = gridView.comm().sum(resultDelta);
             resultDenom = gridView.comm().sum(resultDenom);
 
-	    if (resultDenom > 0.0)
+        if (resultDenom > 0.0)
-	      return resultDelta/resultDenom;
+          return resultDelta/resultDenom;
-	    return 0.0;
+        return 0.0;
         }
 
 
@@ -599,7 +599,6 @@ namespace Opm {
           virtual void apply( const X& x, Y& y ) const
           {
             A_.mv( x, y );
-
             // add well model modification to y
             wellMod_.apply(x, y );
 
@@ -613,7 +612,6 @@ namespace Opm {
           virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
           {
             A_.usmv(alpha,x,y);
-
             // add scaled well model modification to y
             wellMod_.applyScaleAdd( alpha, x, y );
 
@@ -1129,29 +1127,14 @@ namespace Opm {
         const BlackoilAquiferModel<TypeTag>&
         aquiferModel() const { return aquifer_model_; }
 
-        int flowPhaseToEbosCompIdx( const int phaseIdx ) const
+        int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const
         {
-            const auto& pu = phaseUsage_;
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx)
-            if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
+                return Water;
-                return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx)
-            if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
+                return Oil;
-                return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx)
-            if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
+                return Gas;
-                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
@@ -1170,7 +1153,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_; }

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -190,8 +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()) {