Merge pull request #110 from atgeirr/dg-improvements

Improvements for time-of-flight and tracer computations

examples/compute_tof_from_files.cpp

@@ -133,6 +133,10 @@ main(int argc, char** argv)
     } 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);
@@ -158,12 +162,17 @@ main(int argc, char** argv)
     Opm::time::StopWatch transport_timer;
     transport_timer.start();
     std::vector<double> tof;
+    std::vector<double> tracer;
     if (use_dg) {
         Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi, use_limiter);
         tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
     } else {
         Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind);
-        tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof);
+        if (compute_tracer) {
+            tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tof, tracer);
+        } else {
+            tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof);
+        }
     }
     transport_timer.stop();
     double tt = transport_timer.secsSinceStart();
@@ -175,5 +184,14 @@ main(int argc, char** argv)
         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()/grid.number_of_cells;
+            for (int i = 0; i < nt*grid.number_of_cells; ++i) {
+                tracer_stream << tracer[i] << (((i + 1) % nt == 0) ? '\n' : ' ');
+            }
+        }
     }
 }

opm/core/transport/reorder/TransportModelTracerTof.cpp

@@ -32,7 +32,14 @@ namespace Opm
     /// \param[in] grid      A 2d or 3d grid.
     TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid,
                                                      const bool use_multidim_upwind)
-        : grid_(grid), use_multidim_upwind_(use_multidim_upwind)
+        : grid_(grid),
+          darcyflux_(0),
+          porevolume_(0),
+          source_(0),
+          tof_(0),
+          tracer_(0),
+          num_tracers_(0),
+          use_multidim_upwind_(use_multidim_upwind)
     {
     }
 
@@ -68,6 +75,62 @@ namespace Opm
             face_tof_.resize(grid_.number_of_faces);
             std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
         }
+        num_tracers_ = 0;
+        reorderAndTransport(grid_, darcyflux);
+    }
+
+
+
+
+    /// Solve for time-of-flight and a number of tracers.
+    /// One tracer will be used for each inflow flux specified in
+    /// the source parameter.
+    /// \param[in]  darcyflux         Array of signed face fluxes.
+    /// \param[in]  porevolume        Array of pore volumes.
+    /// \param[in]  source            Source term. Sign convention is:
+    ///                                 (+) inflow flux,
+    ///                                 (-) outflow flux.
+    /// \param[out] tof               Array of time-of-flight values (1 per cell).
+    /// \param[out] tracer            Array of tracer values (N per cell, where N is
+    ///                               the number of cells c for which source[c] > 0.0).
+    void TransportModelTracerTof::solveTofTracer(const double* darcyflux,
+                                                 const double* porevolume,
+                                                 const double* source,
+                                                 std::vector<double>& tof,
+                                                 std::vector<double>& tracer)
+    {
+        darcyflux_ = darcyflux;
+        porevolume_ = porevolume;
+        source_ = source;
+#ifndef NDEBUG
+        // Sanity check for sources.
+        const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
+        if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
+            THROW("Sources do not sum to zero: " << cum_src);
+        }
+#endif
+        tof.resize(grid_.number_of_cells);
+        std::fill(tof.begin(), tof.end(), 0.0);
+        tof_ = &tof[0];
+        // Find the tracer heads (injectors).
+        std::vector<int> tracerheads;
+        for (int c = 0; c < grid_.number_of_cells; ++c) {
+            if (source[c] > 0.0) {
+                tracerheads.push_back(c);
+            }
+        }
+        num_tracers_ = tracerheads.size();
+        tracer.resize(grid_.number_of_cells*num_tracers_);
+        std::fill(tracer.begin(), tracer.end(), 0.0);
+        for (int tr = 0; tr < num_tracers_; ++tr) {
+            tracer[tracerheads[tr]*num_tracers_ + tr] = 1.0;
+        }
+        tracer_ = &tracer[0];
+        if (use_multidim_upwind_) {
+            face_tof_.resize(grid_.number_of_faces);
+            std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
+            THROW("Multidimensional upwind not yet implemented for tracer.");
+        }
         reorderAndTransport(grid_, darcyflux);
     }
 
@@ -107,6 +170,9 @@ namespace Opm
                 // face.
                 if (other != -1) {
                     upwind_term += flux*tof_[other];
+                    for (int tr = 0; tr < num_tracers_; ++tr) {
+                        tracer_[num_tracers_*cell + tr] += flux*tracer_[num_tracers_*other + tr];
+                    }
                 }
             } else {
                 downwind_flux += flux;
@@ -115,6 +181,14 @@ namespace Opm
 
         // Compute tof.
         tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
+
+        // Compute tracers (if any).
+        // Do not change tracer solution in source cells.
+        if (source_[cell] <= 0.0) {
+            for (int tr = 0; tr < num_tracers_; ++tr) {
+                tracer_[num_tracers_*cell + tr] *= -1.0/downwind_flux;
+            }
+        }
     }
 
 

opm/core/transport/reorder/TransportModelTracerTof.hpp

@@ -59,6 +59,23 @@ namespace Opm
                       const double* source,
                       std::vector<double>& tof);
 
+        /// Solve for time-of-flight and a number of tracers.
+        /// One tracer will be used for each inflow flux specified in
+        /// the source parameter.
+        /// \param[in]  darcyflux         Array of signed face fluxes.
+        /// \param[in]  porevolume        Array of pore volumes.
+        /// \param[in]  source            Source term. Sign convention is:
+        ///                                 (+) inflow flux,
+        ///                                 (-) outflow flux.
+        /// \param[out] tof               Array of time-of-flight values (1 per cell).
+        /// \param[out] tracer            Array of tracer values (N per cell, where N is
+        ///                               the number of cells c for which source[c] > 0.0).
+        void solveTofTracer(const double* darcyflux,
+                            const double* porevolume,
+                            const double* source,
+                            std::vector<double>& tof,
+                            std::vector<double>& tracer);
+
     private:
         virtual void solveSingleCell(const int cell);
         void solveSingleCellMultidimUpwind(const int cell);
@@ -73,6 +90,8 @@ namespace Opm
         const double* porevolume_;  // one volume per cell
         const double* source_;      // one volumetric source term per cell
         double* tof_;
+        double* tracer_;
+        int num_tracers_;
         bool use_multidim_upwind_;
         std::vector<double> face_tof_;       // For multidim upwind face tofs.
         mutable std::vector<int> adj_faces_; // For multidim upwind logic.

opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp

@@ -404,9 +404,10 @@ namespace Opm
         const int dim = grid_.dimensions;
         const int num_basis = DGBasis::numBasisFunc(dim, degree_);
 
-        double limiter = 1e100;
+        double max_slope_mult = 0.0;
+        int num_upstream_faces = 0;
         // For inflow faces, ensure that cell tof does not dip below
-        // the minimum value from upstream (for that face).
+        // the minimum value from upstream (for all faces).
         for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
             const int face = grid_.cell_faces[hface];
             double flux = 0.0;
@@ -418,6 +419,11 @@ namespace Opm
                 flux = -darcyflux_[face];
                 upstream_cell = grid_.face_cells[2*face];
             }
+            if (flux >= 0.0) {
+                // This is a downstream face.
+                continue;
+            }
+            ++num_upstream_faces;
 
             // Evaluate the solution in all corners, and find the appropriate limiter.
             bool upstream = (upstream_cell >= 0 && flux < 0.0);
@@ -437,26 +443,23 @@ namespace Opm
                     min_upstream = std::min(min_upstream, tof_upstream);
                 }
             }
-            if (min_here < min_upstream) {
-                // Must limit slope.
-                const double tof_c = tof_coeff_[num_basis*cell];
-                if (tof_c < min_upstream) {
-                    // Handle by setting a flat solution.
-                    std::cout << "Trouble in cell " << cell << std::endl;
-                    limiter = 0.0;
-                    tof_coeff_[num_basis*cell] = min_upstream;
-                    break;
-                }
-                const double face_limit = (tof_c - min_upstream)/(tof_c - min_here);
-                limiter = std::min(limiter, face_limit);
+            // Compute maximum slope multiplier.
+            const double tof_c = tof_coeff_[num_basis*cell];
+            if (tof_c < min_upstream) {
+                // Handle by setting a flat solution.
+                std::cout << "Trouble in cell " << cell << std::endl;
+                max_slope_mult = 0.0;
+                tof_coeff_[num_basis*cell] = min_upstream;
+                break;
             }
+            const double face_mult = (tof_c - min_upstream)/(tof_c - min_here);
+            max_slope_mult = std::max(max_slope_mult, face_mult);
         }
+        ASSERT(max_slope_mult >= 0.0);
 
-        if (limiter < 0.0) {
-            THROW("Error in limiter.");
-        }
-
-        if (limiter < 1.0) {
+        // Actually do the limiting (if applicable).
+        const double limiter = max_slope_mult;
+        if (num_upstream_faces > 0 && limiter < 1.0) {
             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_coeff_[i] *= limiter;