From b9a2c141136b1c87e5acf84a7fcc59258f070de7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Sep 2012 14:00:17 +0200 Subject: [PATCH] Add proper support for source terms. This fixes the problem with infinite tofs at sinks. --- .../reorder/TransportModelTracerTof.cpp | 25 ++++++++++++++++--- .../reorder/TransportModelTracerTof.hpp | 9 ++++--- 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 0825706e3..70c2fc400 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -21,6 +21,8 @@ #include #include #include +#include +#include namespace Opm { @@ -37,15 +39,25 @@ namespace Opm /// Solve for time-of-flight at next timestep. - /// \param[in] darcyflux Array of signed face fluxes. - /// \param[in] porevolume Array of pore volumes. - /// \param[in, out] tof Array of time-of-flight values. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[out] tof Array of time-of-flight values. void TransportModelTracerTof::solveTof(const double* darcyflux, const double* porevolume, + const double* source, std::vector& tof) { 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]; @@ -58,8 +70,13 @@ namespace Opm void TransportModelTracerTof::solveSingleCell(const int cell) { // Compute flux terms. + // Sources have zero tof, and therefore do not contribute + // to upwind_term. Sinks on the other hand, must be added + // to the downwind_flux (note sign change resulting from + // different sign conventions: pos. source is injection, + // pos. flux is outflow). double upwind_term = 0.0; - double downwind_flux = 0.0; + double downwind_flux = std::max(-source_[cell], 0.0); for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { int f = grid_.cell_faces[i]; double flux; diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 1ab7ef3b4..21375d895 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -40,11 +40,13 @@ namespace Opm TransportModelTracerTof(const UnstructuredGrid& grid); /// Solve for time-of-flight at next timestep. - /// \param[in] darcyflux Array of signed face fluxes. - /// \param[in] porevolume Array of pore volumes. - /// \param[in, out] tof Array of time-of-flight values. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[out] tof Array of time-of-flight values. void solveTof(const double* darcyflux, const double* porevolume, + const double* source, std::vector& tof); private: @@ -55,6 +57,7 @@ namespace Opm const UnstructuredGrid& grid_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell + const double* source_; // one volumetric source term per cell double* tof_; };