Working and building skeleton which works for multiple aquifers. Also gets the properties of the aquifers using AQUCT keyword.

opm/autodiff/AQUCT_params.hpp

@@ -0,0 +1,31 @@
+#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/BlackoilAquiferModel.hpp

@@ -1,15 +1,13 @@
 /*
-<<<<<<< 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).
 
@@ -62,7 +60,6 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmatrix.hh>
 
-
 #include <opm/material/densead/Math.hpp>
 
 
@@ -74,7 +71,6 @@ namespace Opm {
         
         public:
 
-
             // ---------      Types      ---------
             typedef BlackoilModelParameters ModelParameters;
 
@@ -156,6 +152,7 @@ namespace Opm {
 
             void calculateExplicitQuantities();
 
+
             // The number of components in the model.
             int numComponents() const;
 

opm/autodiff/BlackoilAquiferModel_impl.hpp

@@ -284,4 +284,70 @@ namespace Opm {
         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();
+        
+        if (!deck.hasKeyword("AQUCT")){
+            std::cout << "Nothing is shown! Where is AQUCT!????" << std::endl;
+        }
+
+        const auto& aquctKeyword = deck.getKeyword("AQUCT");
+        std::vector<AQUCT_params> aquctParams;
+        // Resize the parameter vector container based on row entries in aquct
+        // We do the same for aquifers too because number of aquifers is assumed to be for each entry in aquct
+        aquctParams.resize(aquctKeyword.size());
+        // aquifers.resize(aquctKeyword.size());
+
+        const int tableID = 0;
+
+        std::cout << "Parsing AQUCT stuff" << std::endl;
+        for (size_t aquctRecordIdx = 0; aquctRecordIdx < aquctKeyword.size(); ++ aquctRecordIdx) 
+        {
+            const auto& aquctRecord = aquctKeyword.getRecord(aquctRecordIdx);
+
+            aquctParams.at(aquctRecordIdx).aquiferID = aquctRecord.getItem("AQUIFER_ID").template get<int>(0);
+            aquctParams.at(aquctRecordIdx).h = aquctRecord.getItem("THICKNESS_AQ").template get<double>(0);
+            aquctParams.at(aquctRecordIdx).phi_aq = aquctRecord.getItem("PORO_AQ").template get<double>(0);
+            aquctParams.at(aquctRecordIdx).d0 = aquctRecord.getItem("DAT_DEPTH").getSIDouble(0);
+            aquctParams.at(aquctRecordIdx).C_t = aquctRecord.getItem("C_T").template get<double>(0);
+            aquctParams.at(aquctRecordIdx).r_o = aquctRecord.getItem("RAD").getSIDouble(0);
+            aquctParams.at(aquctRecordIdx).k_a = aquctRecord.getItem("PERM_AQ").getSIDouble(0);
+            aquctParams.at(aquctRecordIdx).theta = aquctRecord.getItem("INFLUENCE_ANGLE").template get<double>(0);
+            aquctParams.at(aquctRecordIdx).c1 = 0.008527; // We are using SI
+            aquctParams.at(aquctRecordIdx).c2 = 6.283;
+            aquctParams.at(aquctRecordIdx).inftableID = aquctRecord.getItem("TABLE_NUM_INFLUENCE_FN").template get<int>(0) - 1;
+            aquctParams.at(aquctRecordIdx).pvttableID = aquctRecord.getItem("TABLE_NUM_WATER_PRESS").template get<int>(0) - 1;
+
+            std::cout << aquctParams.at(aquctRecordIdx).inftableID << std::endl;
+            // Get the correct influence table values
+            const auto& aqutabTable = eclState.getTableManager().getAqutabTables().getTable(aquctParams.at(aquctRecordIdx).inftableID);
+            const auto& aqutab_tdColumn = aqutabTable.getColumn(0);
+            const auto& aqutab_piColumn = aqutabTable.getColumn(1);
+            aquctParams.at(aquctRecordIdx).td = aqutab_tdColumn.vectorCopy();
+            aquctParams.at(aquctRecordIdx).pi = aqutab_piColumn.vectorCopy();
+
+            // We determine the cell perforation here.
+            int cellID = 10 + aquctRecordIdx;
+
+            aquctParams.at(aquctRecordIdx).cell_id = cellID;
+
+            // We do not have mu_w as it has to be calculated from pvttable
+            aquifers.push_back(Aquifer_object( aquctParams.at(aquctRecordIdx) ));
+        }
+
+        // I want to deliberately add another aquifer
+        aquifers.push_back( Aquifer_object(99) );
+
+        // aquifers_ = aquifers;
+
+        return aquifers;
+    }
+
 } // namespace Opm