From a48b261a3cd88e3623445488e67a41a80aeba6b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 10 Feb 2012 10:48:18 +0100 Subject: [PATCH] TransportModel* classes are now expected to have a custom solve() method. More: - Using new solve() method in spu_2p. - solve() implemented in terms of protected superclass method reorderAndTransport(). - Removed unused code being replaced by solve(). --- Makefile.am | 8 +-- examples/spu_2p.cpp | 24 ++++----- .../reorder/TransportModelInterface.cpp | 49 +++++++++++++++++++ .../reorder/TransportModelInterface.hpp | 12 +++++ .../reorder/TransportModelTwophase.cpp | 47 ++++++++++-------- .../reorder/TransportModelTwophase.hpp | 20 ++++---- .../transport/reorder/twophasetransport.cpp | 46 ----------------- .../transport/reorder/twophasetransport.hpp | 24 --------- 8 files changed, 116 insertions(+), 114 deletions(-) create mode 100644 opm/core/transport/reorder/TransportModelInterface.cpp delete mode 100644 opm/core/transport/reorder/twophasetransport.cpp delete mode 100644 opm/core/transport/reorder/twophasetransport.hpp diff --git a/Makefile.am b/Makefile.am index 192ebe9e..15305532 100644 --- a/Makefile.am +++ b/Makefile.am @@ -78,11 +78,11 @@ opm/core/pressure/mimetic/hybsys.c \ opm/core/transport/spu_explicit.c \ opm/core/transport/spu_implicit.c \ opm/core/transport/transport_source.c \ +opm/core/transport/reorder/TransportModelInterface.cpp \ opm/core/transport/reorder/TransportModelTwophase.cpp \ opm/core/transport/reorder/reordersequence.c \ opm/core/transport/reorder/nlsolvers.c \ -opm/core/transport/reorder/tarjan.c \ -opm/core/transport/reorder/twophasetransport.cpp +opm/core/transport/reorder/tarjan.c nobase_include_HEADERS = \ @@ -188,8 +188,8 @@ opm/core/transport/reorder/TransportModelInterface.hpp \ opm/core/transport/reorder/TransportModelTwophase.hpp \ opm/core/transport/reorder/nlsolvers.h \ opm/core/transport/reorder/reordersequence.h \ -opm/core/transport/reorder/tarjan.h \ -opm/core/transport/reorder/twophasetransport.hpp +opm/core/transport/reorder/tarjan.h + if UMFPACK libopmcore_la_SOURCES += opm/core/linalg/call_umfpack.c diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 91dd7b64..e5b92f1b 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -59,7 +59,7 @@ #include #include -#include +#include #include #include @@ -753,9 +753,14 @@ main(int argc, char** argv) TwophaseFluid fluid(*props); // Solvers init. + // Pressure solver. PressureSolver psolver(grid->c_grid(), *props); + // Non-reordering solver. TransportModel model (fluid, *grid->c_grid(), porevol, 0, guess_old_solution); TransportSolver tsolver(model); + // Reordering solver. + Opm::TransportModelTwophase reorder_model(*grid->c_grid(), &porevol[0], *props); + // State-related and source-related variables init. std::vector totmob; @@ -763,10 +768,13 @@ main(int argc, char** argv) // We need a separate reorder_sat, because the reorder // code expects a scalar sw, not both sw and so. std::vector reorder_sat(grid->c_grid()->number_of_cells); - double flow_per_sec = 0.1*tot_porevol/Opm::unit::day; + // double flow_per_sec = 0.1*tot_porevol/Opm::unit::day; + double flow_per_sec = 0.1*porevol[0]/Opm::unit::day; std::vector src (grid->c_grid()->number_of_cells, 0.0); - src[0] = flow_per_sec; - src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec; + // src[0] = flow_per_sec; + // src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec; + std::fill(src.begin(), src.begin() + src.size()/2, flow_per_sec); + std::fill(src.begin() + src.size()/2, src.end(), -flow_per_sec); TransportSource* tsrc = create_transport_source(2, 2); double ssrc[] = { 1.0, 0.0 }; double ssink[] = { 0.0, 1.0 }; @@ -840,13 +848,7 @@ main(int argc, char** argv) // boundary flows must be accumulated into // source term following the same convention. transport_timer.start(); - Opm::reorderTransportTwophase(&porevol[0], - &reorder_src[0], - stepsize, - grid->c_grid(), - props.get(), - &state.faceflux()[0], - &reorder_sat[0]); + reorder_model.solve(&state.faceflux()[0], &reorder_src[0], stepsize, &reorder_sat[0]); transport_timer.stop(); double tt = transport_timer.secsSinceStart(); std::cout << "Transport solver took: " << tt << " seconds." << std::endl; diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp new file mode 100644 index 00000000..16ce9c8d --- /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 895ca661..6ad13aac 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 fba4e69d..80dd11b3 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 41846509..f1b88cd7 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 diff --git a/opm/core/transport/reorder/twophasetransport.cpp b/opm/core/transport/reorder/twophasetransport.cpp deleted file mode 100644 index ff069f58..00000000 --- a/opm/core/transport/reorder/twophasetransport.cpp +++ /dev/null @@ -1,46 +0,0 @@ -/* Copyright 2011 (c) Jostein R. Natvig */ - -#include -#include -#include -#include - -#include - - - -void Opm::reorderTransportTwophase(const double *porevolume, - const double *source, - const double dt, - const UnstructuredGrid *grid, - const IncompPropertiesInterface* props, - const double *darcyflux, - double *saturation) -{ - // Set up transport model. - TransportModelTwophase tmodel(grid, props, darcyflux, - porevolume, source, dt, saturation); - - // 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 - tmodel.solveSingleCell(sequence[i]); - } -} - -/* Local Variables: */ -/* c-basic-offset:4 */ -/* End: */ diff --git a/opm/core/transport/reorder/twophasetransport.hpp b/opm/core/transport/reorder/twophasetransport.hpp deleted file mode 100644 index 610437bc..00000000 --- a/opm/core/transport/reorder/twophasetransport.hpp +++ /dev/null @@ -1,24 +0,0 @@ -/* Copyright 2011 (c) Jostein R. Natvig */ - -#ifndef TWOPHASETRANSPORT_HPP_INCLUDED -#define TWOPHASETRANSPORT_HPP_INCLUDED - - -struct UnstructuredGrid; - -namespace Opm -{ - - class IncompPropertiesInterface; - - void reorderTransportTwophase(const double *porevolume, - const double *source, - const double dt, - const UnstructuredGrid *grid, - const IncompPropertiesInterface* props, - const double *darcyflux, - double *saturation); -} // namespace Opm - - -#endif /* TWOPHASETRANSPORT_HPP_INCLUDED */