Merge pull request #1571 from kel85uk/fix_#14_aquifer_making_connections

Fixes for carter tracy aquifer model

opm/autodiff/AquiferCarterTracy.hpp

@@ -235,6 +235,31 @@ namespace Opm
                       / ( aquct_data_.k_a * aquct_data_.c1 );
             }
 
+            template<class faceCellType, class ugridType>
+            inline const double getFaceArea(const faceCellType& faceCells, const ugridType& ugrid,
+                                            const int faceIdx, const int idx,
+                                            const Aquancon::AquanconOutput& connection) const
+            {
+                // 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)
+                double faceArea = 0.;
+                const auto cellNeighbour0 = faceCells(faceIdx,0);
+                const auto cellNeighbour1 = faceCells(faceIdx,1);
+                const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
+                const auto calculatedFaceArea = (!connection.influx_coeff.at(idx))? 
+                                                defaultFaceArea : 
+                                                *(connection.influx_coeff.at(idx));
+                faceArea = (cellNeighbour0 * cellNeighbour1 > 0)? 0. : calculatedFaceArea;
+                if (cellNeighbour1 == 0){
+                    faceArea = (cellNeighbour0 < 0)? faceArea : 0.;
+                }
+                else if (cellNeighbour0 == 0){
+                    faceArea = (cellNeighbour1 < 0)? faceArea : 0.;
+                }
+                return faceArea;
+            }
+
             // 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)
             {
@@ -246,7 +271,7 @@ namespace Opm
                 auto globalCellIdx = ugrid.globalCell();
 
                 assert( cell_idx_ == connection.global_index);
-                assert( (cell_idx_.size() == connection.influx_coeff.size()) );
+                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()) );
 
@@ -297,11 +322,7 @@ namespace Opm
 
                         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
+                            faceArea_connected_.at(idx) = getFaceArea(faceCells, ugrid, faceIdx, idx, connection);
-                            // 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) );
                         }
                     }
@@ -309,9 +330,13 @@ namespace Opm
                     cell_depth_.at(idx) = cellCenter[2];
                 }
 
+                const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
+
                 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;
+                    alphai_.at(idx) = (denom_face_areas < eps_sqrt)? // Prevent no connection NaNs due to division by zero
+                                      0.
+                                    : ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas;
                 }
             }
 
@@ -322,13 +347,13 @@ namespace Opm
                 
                 rhow_.resize(cell_idx_.size(),0.);
                 
-                if (aquct_data_.p0 < 1.0)
+                if (!aquct_data_.p0)
                 {
                    pa0_ = calculateReservoirEquilibrium();
                 }
                 else
                 {
-                   pa0_ = aquct_data_.p0;
+                   pa0_ = *(aquct_data_.p0);
                 }
 
                 // use the thermodynamic state of the first active cell as a