diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp new file mode 100644 index 000000000..16ce9c8d0 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelInterface.cpp @@ -0,0 +1,49 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include + +#include +#include + + +void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux) +{ + // Compute sequence of single-cell problems + std::vector sequence(grid.number_of_cells); + std::vector components(grid.number_of_cells + 1); + int ncomponents; + compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents); + + // Assume all strong components are single-cell domains. + assert(ncomponents == grid.number_of_cells); + for (int i = 0; i < grid.number_of_cells; ++i) { +#ifdef MATLAB_MEX_FILE + if (interrupt_signal) { + mexPrintf("Reorder loop interrupted by user: %d of %d " + "cells finished.\n", i, grid.number_of_cells); + break; + } +#endif + solveSingleCell(sequence[i]); + } +} diff --git a/opm/core/transport/reorder/TransportModelInterface.hpp b/opm/core/transport/reorder/TransportModelInterface.hpp index 895ca6618..6ad13aac6 100644 --- a/opm/core/transport/reorder/TransportModelInterface.hpp +++ b/opm/core/transport/reorder/TransportModelInterface.hpp @@ -20,16 +20,28 @@ #ifndef OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED #define OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED +class UnstructuredGrid; + namespace Opm { + /// Interface for reordering transport models. + /// A transport model must provide the solveSingleCell() + /// method, and is expected to implement a solve() method + /// that will have an interface geared to the model's needs. + /// (The solve() method is therefore not virtual in this class). + /// The reorderAndTransport() method is provided as an + /// aid to implementing solve() in subclasses. class TransportModelInterface { public: virtual ~TransportModelInterface() {} virtual void solveSingleCell(int cell) = 0; + protected: + void reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux); }; + } // namespace Opm #endif // OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index fba4e69db..80dd11b37 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -27,28 +27,35 @@ namespace Opm { - TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid* grid, - const Opm::IncompPropertiesInterface* props, - const double* darcyflux, + TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid& grid, const double* porevolume, - const double* source, - const double dt, - double* saturation) + const Opm::IncompPropertiesInterface& props) : grid_(grid), - props_(props), - darcyflux_(darcyflux), porevolume_(porevolume), - source_(source), - dt_(dt), - saturation_(saturation), - fractionalflow_(grid->number_of_cells, 0.0) + props_(props), + darcyflux_(0), + source_(0), + dt_(0.0), + saturation_(0), + fractionalflow_(grid.number_of_cells, -1.0) { - if (props->numPhases() != 2) { + if (props.numPhases() != 2) { THROW("Property object must have 2 phases"); } - visc_ = props->viscosity(); + visc_ = props.viscosity(); } + void TransportModelTwophase::solve(const double* darcyflux, + const double* source, + const double dt, + double* saturation) + { + darcyflux_ = darcyflux; + source_ = source; + dt_ = dt; + saturation_ = saturation; + reorderAndTransport(grid_, darcyflux); + } // Residual function r(s) for a single-cell implicit Euler transport // @@ -73,17 +80,17 @@ namespace Opm outflux = tm.source_[cell] <= 0 ? -tm.source_[cell] : 0.0; dtpv = tm.dt_/tm.porevolume_[cell]; - for (int i = tm.grid_->cell_facepos[cell]; i < tm.grid_->cell_facepos[cell+1]; ++i) { - int f = tm.grid_->cell_faces[i]; + for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { + int f = tm.grid_.cell_faces[i]; double flux; int other; // Compute cell flux - if (cell == tm.grid_->face_cells[2*f]) { + if (cell == tm.grid_.face_cells[2*f]) { flux = tm.darcyflux_[f]; - other = tm.grid_->face_cells[2*f+1]; + other = tm.grid_.face_cells[2*f+1]; } else { flux =-tm.darcyflux_[f]; - other = tm.grid_->face_cells[2*f]; + other = tm.grid_.face_cells[2*f]; } // Add flux to influx or outflux, if interiour. if (other != -1) { @@ -118,7 +125,7 @@ namespace Opm { double sat[2] = { s, 1.0 - s }; double mob[2]; - props_->relperm(1, sat, &cell, mob, 0); + props_.relperm(1, sat, &cell, mob, 0); mob[0] /= visc_[0]; mob[1] /= visc_[1]; return mob[0]/(mob[0] + mob[1]); diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index 418465092..f1b88cd77 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelTwophase.hpp @@ -33,21 +33,23 @@ namespace Opm class TransportModelTwophase : public TransportModelInterface { public: - TransportModelTwophase(const UnstructuredGrid* grid, - const Opm::IncompPropertiesInterface* props, - const double* darcyflux, + TransportModelTwophase(const UnstructuredGrid& grid, const double* porevolume, - const double* source, - const double dt, - double* saturation); + const Opm::IncompPropertiesInterface& props); + + void solve(const double* darcyflux, + const double* source, + const double dt, + double* saturation); virtual void solveSingleCell(int cell); private: - const UnstructuredGrid* grid_; - const IncompPropertiesInterface* props_; - const double* darcyflux_; // one flux per grid face + const UnstructuredGrid& grid_; const double* porevolume_; // one volume per cell + const IncompPropertiesInterface& props_; + + const double* darcyflux_; // one flux per grid face const double* source_; // one source per cell double dt_; double* saturation_; // one per cell