diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index f1158598..9bc67bfa 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -608,17 +608,23 @@ main(int argc, char** argv) std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; ptime += pt; + // Process transport sources (to include bdy terms). + if (use_reorder) { + Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, reorder_src); + } else { + clear_transport_source(tsrc); + for (int cell = 0; cell < num_cells; ++cell) { + if (src[cell] > 0.0) { + append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc); + } else if (src[cell] < 0.0) { + append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc); + } + } + } + // Solve transport transport_timer.start(); if (use_reorder) { - // We must treat reorder_src here, - // if we are to handle anything but simple water - // injection, since it is expected to be - // equal to total outflow (if negative) - // and water inflow (if positive). - // Also, for anything but noflow boundaries, - // boundary flows must be accumulated into - // source term following the same convention. Opm::toWaterSat(state.saturation(), reorder_sat); reorder_model.solve(&state.faceflux()[0], &reorder_src[0], stepsize, &reorder_sat[0]); Opm::toBothSat(reorder_sat, state.saturation()); diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 11c3b0df..53a4104a 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -141,6 +141,48 @@ namespace Opm } } + /// Compute two-phase transport source terms from face fluxes, + /// and pressure equation source terms. This puts boundary flows + /// into the source terms for the transport equation. + /// \param[in] grid The grid used. + /// \param[in] src Pressure eq. source terms. The sign convention is: + /// (+) positive total inflow (positive velocity divergence) + /// (-) negative total outflow + /// \param[in] faceflux Signed face fluxes, typically the result from a flow solver. + /// \param[in] inflow_frac Fraction of inflow that consists of first phase. + /// Example: if only water is injected, inflow_frac == 1.0. + /// Note: it is not possible (with this method) to use different fractions + /// for different inflow sources, be they source terms of boundary flows. + /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign: + /// (+) positive inflow of first phase (water) + /// (-) negative total outflow of both phases + void computeTransportSource(const UnstructuredGrid& grid, + const std::vector& src, + const std::vector& faceflux, + const double inflow_frac, + std::vector& transport_src) + { + int nc = grid.number_of_cells; + transport_src.resize(nc); + for (int c = 0; c < nc; ++c) { + transport_src[c] = 0.0; + transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; + for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* f2c = &grid.face_cells[2*f]; + double bdy_influx = 0.0; + if (f2c[0] == c && f2c[1] == -1) { + bdy_influx = -faceflux[f]; + } else if (f2c[0] == -1 && f2c[1] == c) { + bdy_influx = faceflux[f]; + } + if (bdy_influx != 0.0) { + transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; + } + } + } + } + /// @brief Estimates a scalar cell velocity from face fluxes. /// @param[in] grid a grid /// @param[in] face_flux signed per-face fluxes diff --git a/opm/core/utility/miscUtilities.hpp b/opm/core/utility/miscUtilities.hpp index dc061683..22be2383 100644 --- a/opm/core/utility/miscUtilities.hpp +++ b/opm/core/utility/miscUtilities.hpp @@ -73,6 +73,27 @@ namespace Opm std::vector& totmob, std::vector& omega); + /// Compute two-phase transport source terms from face fluxes, + /// and pressure equation source terms. This puts boundary flows + /// into the source terms for the transport equation. + /// \param[in] grid The grid used. + /// \param[in] src Pressure eq. source terms. The sign convention is: + /// (+) positive total inflow (positive velocity divergence) + /// (-) negative total outflow + /// \param[in] faceflux Signed face fluxes, typically the result from a flow solver. + /// \param[in] inflow_frac Fraction of inflow that consists of first phase. + /// Example: if only water is injected, inflow_frac == 1.0. + /// Note: it is not possible (with this method) to use different fractions + /// for different inflow sources, be they source terms of boundary flows. + /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign: + /// (+) positive inflow of first phase (water) + /// (-) negative total outflow of both phases + void computeTransportSource(const UnstructuredGrid& grid, + const std::vector& src, + const std::vector& faceflux, + const double inflow_frac, + std::vector& transport_src); + /// @brief Estimates a scalar cell velocity from face fluxes. /// @param[in] grid a grid /// @param[in] face_flux signed per-face fluxes @@ -91,6 +112,7 @@ namespace Opm void toBothSat(const std::vector& sw, std::vector& sboth); + } // namespace Opm #endif // OPM_MISCUTILITIES_HEADER_INCLUDED