Merge remote-tracking branch 'upstream/master'

opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp

@@ -20,10 +20,11 @@
 #include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/VelocityInterpolation.hpp>
 #include <opm/core/linalg/blas_lapack.h>
 #include <algorithm>
-#include <numeric>
 #include <cmath>
+#include <numeric>
 
 namespace Opm
 {
@@ -403,32 +404,6 @@ namespace Opm
 
 
 
-    // Initial version: only a constant interpolation.
-    static void interpolateVelocity(const UnstructuredGrid& grid,
-                                    const int cell,
-                                    const double* darcyflux,
-                                    const double* /*x*/,
-                                    double* v)
-    {
-        const int dim = grid.dimensions;
-        std::fill(v, v + dim, 0.0);
-        const double* cc = grid.cell_centroids + cell*dim;
-        for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) {
-            const int face = grid.cell_faces[hface];
-            const double* fc = grid.face_centroids + face*dim;
-            double flux = 0.0;
-            if (cell == grid.face_cells[2*face]) {
-                flux = darcyflux[face];
-            } else {
-                flux = -darcyflux[face];
-            }
-            for (int dd = 0; dd < dim; ++dd) {
-                v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell];
-            }
-        }
-    }
-
-
 
     // ---------------   Methods of TransportModelTracerTofDiscGal ---------------
 
@@ -436,11 +411,19 @@ namespace Opm
 
     /// Construct solver.
     /// \param[in] grid      A 2d or 3d grid.
-    TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid)
+    /// \param[in] use_cvi   If true, use corner point velocity interpolation.
+    ///                      Otherwise, use the basic constant interpolation.
+    TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
+                                                                   const bool use_cvi)
         : grid_(grid),
           coord_(grid.dimensions),
           velocity_(grid.dimensions)
     {
+        if (use_cvi) {
+            velocity_interpolation_.reset(new VelocityInterpolationECVI(grid));
+        } else {
+            velocity_interpolation_.reset(new VelocityInterpolationConstant(grid));
+        }
     }
 
 
@@ -482,9 +465,11 @@ namespace Opm
         tof_coeff_ = &tof_coeff[0];
         rhs_.resize(num_basis);
         jac_.resize(num_basis*num_basis);
+        orig_jac_.resize(num_basis*num_basis);
         basis_.resize(num_basis);
         basis_nb_.resize(num_basis);
         grad_basis_.resize(num_basis*grid_.dimensions);
+        velocity_interpolation_->setupFluxes(darcyflux);
         reorderAndTransport(grid_, darcyflux);
     }
 
@@ -560,7 +545,7 @@ namespace Opm
                 quad.quadPtCoord(quad_pt, &coord_[0]);
                 DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
                 DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]);
-                interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]);
+                velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]);
                 const double w = quad.quadPtWeight(quad_pt);
                 for (int j = 0; j < num_basis; ++j) {
                     for (int i = 0; i < num_basis; ++i) {
@@ -633,6 +618,7 @@ namespace Opm
         std::vector<MAT_SIZE_T> piv(num_basis);
         MAT_SIZE_T ldb = num_basis;
         MAT_SIZE_T info = 0;
+        orig_jac_ = jac_;
         dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
         if (info != 0) {
             // Print the local matrix and rhs.
@@ -640,7 +626,7 @@ namespace Opm
                       << " with A = \n";
             for (int row = 0; row < n; ++row) {
                 for (int col = 0; col < n; ++col) {
-                    std::cerr << "    " << jac_[row + n*col];
+                    std::cerr << "    " << orig_jac_[row + n*col];
                 }
                 std::cerr << '\n';
             }

opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp

@@ -21,15 +21,18 @@
 #define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
 
 #include <opm/core/transport/reorder/TransportModelInterface.hpp>
+#include <boost/shared_ptr.hpp>
 #include <vector>
 #include <map>
 #include <ostream>
+
 struct UnstructuredGrid;
 
 namespace Opm
 {
 
     class IncompPropertiesInterface;
+    class VelocityInterpolationInterface;
 
     /// Implements a discontinuous Galerkin solver for
     /// (single-phase) time-of-flight using reordering.
@@ -45,7 +48,10 @@ namespace Opm
     public:
         /// Construct solver.
         /// \param[in] grid      A 2d or 3d grid.
-        TransportModelTracerTofDiscGal(const UnstructuredGrid& grid);
+        /// \param[in] use_cvi   If true, use corner point velocity interpolation.
+        ///                      Otherwise, use the basic constant interpolation.
+        TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
+                                       const bool use_cvi);
 
 
         /// Solve for time-of-flight.
@@ -72,7 +78,12 @@ namespace Opm
         virtual void solveMultiCell(const int num_cells, const int* cells);
 
     private:
+        // Disable copying and assignment.
+        TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&);
+        TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&);
+
         const UnstructuredGrid& grid_;
+        boost::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
         const double* darcyflux_;   // one flux per grid face
         const double* porevolume_;  // one volume per cell
         const double* source_;      // one volumetric source term per cell
@@ -80,6 +91,7 @@ namespace Opm
         double* tof_coeff_;
         std::vector<double> rhs_;   // single-cell right-hand-side
         std::vector<double> jac_;   // single-cell jacobian
+        std::vector<double> orig_jac_;   // single-cell jacobian (copy)
         // Below: storage for quantities needed by solveSingleCell().
         std::vector<double> coord_;
         std::vector<double> basis_;

opm/core/utility/writeVtkData.hpp

@@ -26,25 +26,23 @@
 #include <vector>
 #include <tr1/array>
 #include <iosfwd>
+#include <opm/core/utility/DataMap.hpp>
 
 struct UnstructuredGrid;
 
 namespace Opm
 {
 
-    /// Intended to map strings (giving the output field names) to data.
-    typedef std::map<std::string, const std::vector<double>*> DataMap;
-
     /// Vtk output for cartesian grids.
     void writeVtkData(const std::tr1::array<int, 3>& dims,
-		      const std::tr1::array<double, 3>& cell_size,
+                      const std::tr1::array<double, 3>& cell_size,
-		      const DataMap& data,
+                      const DataMap& data,
-		      std::ostream& os);
+                      std::ostream& os);
 
     /// Vtk output for general grids.
     void writeVtkData(const UnstructuredGrid& grid,
-		      const DataMap& data,
+                      const DataMap& data,
-		      std::ostream& os);
+                      std::ostream& os);
 } // namespace Opm
 
 #endif // OPM_WRITEVTKDATA_HEADER_INCLUDED

opm/core/wells/WellsGroup.cpp

@@ -761,16 +761,17 @@ namespace Opm
             wells_->ctrls[self_index_]->current = ~ wells_->ctrls[self_index_]->current;
         }
     }
-    
+
     std::pair<WellNode*, double> WellNode::getWorstOffending(const std::vector<double>& well_reservoirrates_phase,
                                                              const std::vector<double>& well_surfacerates_phase,
                                                              ProductionSpecification::ControlMode mode)
     {
         const int np = phaseUsage().num_phases;
         const int index = self_index_*np;
-        return std::make_pair<WellNode*, double>(this, rateByMode(&well_reservoirrates_phase[index],
+        return std::pair<WellNode*, double>(this,
-                                                                  &well_surfacerates_phase[index],
+                                            rateByMode(&well_reservoirrates_phase[index],
-                                                                  mode));
+                                                       &well_surfacerates_phase[index],
+                                                       mode));
     }
     
     void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
@@ -782,7 +783,7 @@ namespace Opm
              && (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) {
             return;
         }
-        if (!wells_->type[self_index_] == INJECTOR) {
+        if (wells_->type[self_index_] != INJECTOR) {
             ASSERT(target == 0.0);
             return;
         }
@@ -859,7 +860,7 @@ namespace Opm
             std::cout << "Returning" << std::endl;
             return;
         }
-        if (!wells_->type[self_index_] == PRODUCER) {
+        if (wells_->type[self_index_] != PRODUCER) {
             ASSERT(target == 0.0);
             return;
         }