finishing the function addWellContributions for StandardWellV

opm/autodiff/MatrixBlock.hpp

@@ -329,6 +329,61 @@ namespace Detail
             }
         }
     }
+
+    //! calculates ret = A * B
+    template< class K>
+    static inline void multMatrix( const Dune::DynamicMatrix<K> &A,
+                                   const Dune::DynamicMatrix<K> &B,
+                                   Dune::DynamicMatrix<K> &ret )
+    {
+        typedef typename Dune::DynamicMatrix<K> :: size_type size_type;
+
+        const size_type m = A.rows();
+        const size_type n = A.cols();
+
+        assert(n == B.rows() );
+
+        const size_type p = B.cols();
+
+        ret.resize(m, p);
+
+        for( size_type i = 0; i < m; ++i )
+        {
+            for( size_type j = 0; j < p; ++j )
+            {
+                ret[ i ][ j ] = K( 0 );
+                for( size_type k = 0; k < n; ++k )
+                    ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
+            }
+        }
+    }
+
+
+    //! calculates ret = A^T * B
+    template< class K, int m, int p >
+    static inline void multMatrixTransposed ( const Dune::DynamicMatrix<K> &A,
+                                              const Dune::DynamicMatrix<K> &B,
+                                              Dune::FieldMatrix< K, m, p> &ret )
+    {
+        typedef typename Dune::DynamicMatrix<K> :: size_type size_type;
+
+        // A is a tranpose matrix
+        const size_type n = A.rows();
+        assert(m == A.cols() );
+
+        assert(n == B.rows() );
+        assert(p == B.cols() );
+
+        for( size_type i = 0; i < m; ++i )
+        {
+            for( size_type j = 0; j < p; ++j )
+            {
+                ret[ i ][ j ] = K( 0 );
+                for( size_type k = 0; k < n; ++k )
+                    ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ];
+            }
+        }
+    }
 } // namespace Detail
 } // namespace Opm
 

opm/autodiff/StandardWellV_impl.hpp

@@ -2741,7 +2741,7 @@ namespace Opm
         // at (0,j) only if this well has a perforation at cell j.
 
         // TODO: we have some matrix type problem here.
-        /* for ( auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC )
+        for ( auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC )
         {
             const auto row_index = colC.index();
             auto& row = mat[row_index];
@@ -2754,13 +2754,13 @@ namespace Opm
                 while ( col != row.end() && col.index() < col_index ) ++col;
                 assert(col != row.end() && col.index() == col_index);
 
-                Dune::FieldMatrix<Scalar, numWellEq, numEq> tmp;
+                Dune::DynamicMatrix<Scalar> tmp;
+                Detail::multMatrix(invDuneD_[0][0],  (*colB), tmp);
                 typename Mat::block_type tmp1;
-                Dune::FMatrixHelp::multMatrix(invDuneD_[0][0],  (*colB), tmp);
                 Detail::multMatrixTransposed((*colC), tmp, tmp1);
                 (*col) -= tmp1;
             }
-        } */
+        }
     }
 
 

opm/autodiff/VFPHelpers.hpp

@@ -62,7 +62,9 @@ inline EvalWell zeroIfNanInf(const EvalWell& value) {
         OpmLog::warning("NAN_OR_INF_VFP_EVAL", "NAN or INF Evalution encountered during VFP calculation, the Evalution is set to zero");
     }
 
-    return nan_or_inf ? 0.0 * value : value;
+    using Toolbox = MathToolbox<EvalWell>;
+
+    return nan_or_inf ? Toolbox::createBlank(value) : value;
 }