Change constructor to take ParameterGroup argument.

Also make tof to limit against >= 0.0, for case when upstream cell values go
below zero.

Disabled some debug output.

examples/compute_tof.cpp

@@ -169,16 +169,19 @@ main(int argc, char** argv)
     // Choice of tof solver.
     bool use_dg = param.getDefault("use_dg", false);
     int dg_degree = -1;
-    bool use_cvi = false;
-    bool use_limiter = false;
     bool use_multidim_upwind = false;
+    // Need to initialize dg solver here, since it uses parameters now.
+    boost::scoped_ptr<Opm::TransportModelTracerTofDiscGal> dg_solver;
     if (use_dg) {
         dg_degree = param.getDefault("dg_degree", 0);
-        use_cvi = param.getDefault("use_cvi", false);
+        dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(*grid->c_grid(), param));
-        use_limiter = param.getDefault("use_limiter", false);
     } else {
         use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
     }
+    bool compute_tracer = param.getDefault("compute_tracer", false);
+    if (use_dg && compute_tracer) {
+        THROW("DG for tracer not yet implemented.");
+    }
 
     // Write parameters used for later reference.
     bool output = param.getDefault("output", true);
@@ -232,12 +235,16 @@ main(int argc, char** argv)
     // Solve time-of-flight.
     transport_timer.start();
     std::vector<double> tof;
+    std::vector<double> tracer;
     if (use_dg) {
-        Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi, use_limiter);
+        dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
-        tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
     } else {
         Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind);
-        tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
+        if (compute_tracer) {
+            tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tof, tracer);
+        } else {
+            tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
+        }
     }
     transport_timer.stop();
     double tt = transport_timer.secsSinceStart();
@@ -249,7 +256,17 @@ main(int argc, char** argv)
     if (output) {
         std::string tof_filename = output_dir + "/tof.txt";
         std::ofstream tof_stream(tof_filename.c_str());
+        tof_stream.precision(16);
         std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
+        if (compute_tracer) {
+            std::string tracer_filename = output_dir + "/tracer.txt";
+            std::ofstream tracer_stream(tracer_filename.c_str());
+            tracer_stream.precision(16);
+            const int nt = tracer.size()/num_cells;
+            for (int i = 0; i < nt*num_cells; ++i) {
+                tracer_stream << tracer[i] << (((i + 1) % nt == 0) ? '\n' : ' ');
+            }
+        }
     }
 
     std::cout << "\n\n================    End of simulation     ===============\n"

examples/compute_tof_from_files.cpp

@@ -123,13 +123,12 @@ main(int argc, char** argv)
     // Choice of tof solver.
     bool use_dg = param.getDefault("use_dg", false);
     int dg_degree = -1;
-    bool use_cvi = false;
-    bool use_limiter = false;
     bool use_multidim_upwind = false;
+    // Need to initialize dg solver here, since it uses parameters now.
+    boost::scoped_ptr<Opm::TransportModelTracerTofDiscGal> dg_solver;
     if (use_dg) {
         dg_degree = param.getDefault("dg_degree", 0);
-        use_cvi = param.getDefault("use_cvi", false);
+        dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(grid, param));
-        use_limiter = param.getDefault("use_limiter", false);
     } else {
         use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
     }
@@ -164,8 +163,7 @@ main(int argc, char** argv)
     std::vector<double> tof;
     std::vector<double> tracer;
     if (use_dg) {
-        Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi, use_limiter);
+        dg_solver->solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
-        tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
     } else {
         Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind);
         if (compute_tracer) {

opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp

@@ -23,6 +23,7 @@
 #include <opm/core/grid.h>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/VelocityInterpolation.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/linalg/blas_lapack.h>
 #include <algorithm>
 #include <cmath>
@@ -128,27 +129,50 @@ namespace Opm
     /// \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,
+                                                                   const parameter::ParameterGroup& param)
-                                                                   const bool use_limiter)
         : grid_(grid),
-          use_cvi_(use_cvi),
+          use_cvi_(false),
-          use_limiter_(use_limiter),
+          use_limiter_(false),
           limiter_relative_flux_threshold_(1e-3),
           limiter_method_(MinUpwindFace),
           limiter_usage_(DuringComputations),
           coord_(grid.dimensions),
           velocity_(grid.dimensions)
     {
+        use_cvi_ = param.getDefault("use_cvi", use_cvi_);
+        use_limiter_ = param.getDefault("use_limiter", use_limiter_);
+        if (use_limiter_) {
+            limiter_relative_flux_threshold_ = param.getDefault("limiter_relative_flux_threshold",
+                                                                limiter_relative_flux_threshold_);
+            const std::string limiter_method_str = param.getDefault<std::string>("limiter_method", "MinUpwindFace");
+            if (limiter_method_str == "MinUpwindFace") {
+                limiter_method_ = MinUpwindFace;
+            } else if (limiter_method_str == "MinUpwindAverage") {
+                limiter_method_ = MinUpwindAverage;
+            } else {
+                THROW("Unknown limiter method: " << limiter_method_str);
+            }
+            const std::string limiter_usage_str = param.getDefault<std::string>("limiter_usage", "DuringComputations");
+            if (limiter_usage_str == "DuringComputations") {
+                limiter_usage_ = DuringComputations;
+            } else if (limiter_usage_str == "AsPostProcess") {
+                limiter_usage_ = AsPostProcess;
+            } else if (limiter_usage_str == "AsSimultaneousPostProcess") {
+                limiter_usage_ = AsSimultaneousPostProcess;
+            } else {
+                THROW("Unknown limiter usage spec: " << limiter_usage_str);
+            }
+        }
         // A note about the use_cvi_ member variable:
         // In principle, we should not need it, since the choice of velocity
         // interpolation is made below, but we may need to use higher order
         // quadrature to exploit CVI, so we store the choice.
         // An alternative would be to add a virtual method isConstant() to
         // the VelocityInterpolationInterface.
-        if (use_cvi) {
+        if (use_cvi_) {
-            velocity_interpolation_.reset(new VelocityInterpolationECVI(grid));
+            velocity_interpolation_.reset(new VelocityInterpolationECVI(grid_));
         } else {
-            velocity_interpolation_.reset(new VelocityInterpolationConstant(grid));
+            velocity_interpolation_.reset(new VelocityInterpolationConstant(grid_));
         }
     }
 
@@ -456,7 +480,7 @@ namespace Opm
         // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
         const double total_flux = std::max(-upstream_flux, downstream_flux);
 
-        // Find minimum tof on upstream faces.
+        // Find minimum tof on upstream faces and for this cell.
         const int dim = grid_.dimensions;
         const int num_basis = DGBasis::numBasisFunc(dim, degree_);
         double min_upstream_tof = 1e100;
@@ -506,6 +530,9 @@ namespace Opm
             min_upstream_tof = 0.0;
             min_here_tof = 0.0;
         }
+        if (min_upstream_tof < 0.0) {
+            min_upstream_tof = 0.0;
+        }
         const double tof_c = tof_coeff_[num_basis*cell];
         double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
         if (tof_c < min_upstream_tof) {
@@ -518,12 +545,12 @@ namespace Opm
 
         // Actually do the limiting (if applicable).
         if (limiter < 1.0) {
-            std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
+            // std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
             for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) {
                 tof[i] *= limiter;
             }
         } else {
-            std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
+            // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
         }
     }
 

opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp

@@ -33,6 +33,7 @@ namespace Opm
 
     class IncompPropertiesInterface;
     class VelocityInterpolationInterface;
+    namespace parameter { class ParameterGroup; }
 
     /// Implements a discontinuous Galerkin solver for
     /// (single-phase) time-of-flight using reordering.
@@ -51,8 +52,7 @@ namespace Opm
         /// \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,
+                                       const parameter::ParameterGroup& param);
-                                       const bool use_limiter = false);
 
 
         /// Solve for time-of-flight.