diff --git a/tutorials/tutorial4.cpp b/tutorials/tutorial4.cpp index 7bdbd1b93..9358da13a 100644 --- a/tutorials/tutorial4.cpp +++ b/tutorials/tutorial4.cpp @@ -35,14 +35,7 @@ #include #include -#include -#include -#include -#include -#include -#include -#include -#include +#include #include @@ -52,113 +45,6 @@ #include /// \page tutorial4 Well controls -/// -/// \page tutorial4 -/// \section commentedcode4 Commented code: -/// \page tutorial4 -/// \details -/// We define a class which computes mobility, capillary pressure and -/// the minimum and maximum saturation value for each cell. -/// \code -class TwophaseFluid -{ -public: - /// \endcode - /// \page tutorial4 - /// \details Constructor operator. Takes in the fluid properties defined - /// \c props - /// \code - TwophaseFluid(const Opm::IncompPropertiesInterface& props); - /// \endcode - - /// \page tutorial4 - /// \details Density for each phase. - /// \code - double density(int phase) const; - /// \endcode - - /// \page tutorial4 - /// \details Computes the mobility and its derivative with respect to saturation - /// for a given cell \c c and saturation \c sat. The template classes \c Sat, - /// \c Mob, \c DMob are typically arrays. By using templates, we do not have to - /// investigate how these array objects are implemented - /// (as long as they have an \c operator[] method). - /// \code - template - void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const; - /// \endcode - - /// \page tutorial4 - /// \details Computes the capillary pressure and its derivative with respect - /// to saturation - /// for a given cell \c c and saturation \c sat. - /// \code - template - void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const; - /// \endcode - - /// \page tutorial4 - /// \details Returns the minimum permitted saturation. - /// \code - double s_min(int c) const; - /// \endcode - - /// \page tutorial4 - /// \details Returns the maximum permitted saturation - /// \code - double s_max(int c) const; - /// \endcode - - /// \page tutorial4 - /// \details Private variables - /// \code -private: - const Opm::IncompPropertiesInterface& props_; - std::vector smin_; - std::vector smax_; -}; -/// \endcode - -/// \page tutorial4 -/// \details We set up the transport model. -/// \code -typedef Opm::SinglePointUpwindTwoPhase TransportModel; -/// \endcode - -/// \page tutorial4 -/// \details -/// The transport equation is nonlinear. We use an implicit transport solver -/// which implements a Newton-Raphson solver. -/// We define the format of the objects -/// which will be used by the solver. -/// \code -using namespace Opm::ImplicitTransportDefault; -typedef NewtonVectorCollection< ::std::vector > NVecColl; -typedef JacobianSystem< struct CSRMatrix, NVecColl > JacSys; - -template -class MaxNorm { -public: - static double - norm(const Vector& v) { - return AccumulationNorm ::norm(v); - } -}; -/// \endcode - -/// \page tutorial4 -/// \details -/// We set up the solver. -/// \code -typedef Opm::ImplicitTransport TransportSolver; -/// \endcode - /// \page tutorial4 /// \details @@ -179,8 +65,9 @@ int main () double dy = 10.; double dz = 10.; using namespace Opm; - GridManager grid(nx, ny, nz, dx, dy, dz); - int num_cells = grid.c_grid()->number_of_cells; + GridManager grid_manager(nx, ny, nz, dx, dy, dz); + const UnstructuredGrid& grid = *grid_manager.c_grid(); + int num_cells = grid.number_of_cells; /// \endcode /// \page tutorial4 @@ -214,8 +101,7 @@ int main () /// description and set this function to be linear. For more realistic fluid, the /// saturation function is given by the data. /// \code - SaturationPropsBasic::RelPermFunc rel_perm_func; - rel_perm_func = SaturationPropsBasic::Linear; + SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear; /// \endcode /// \page tutorial4 @@ -224,7 +110,6 @@ int main () /// \code IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu, porosity, k, dim, num_cells); - TwophaseFluid fluid(props); /// \endcode @@ -236,15 +121,6 @@ int main () /// \endcode - - - /// \page tutorial4 - /// \details We set up the source term - /// \code - std::vector src(num_cells, 0.0); - src[0] = 1.; - src[num_cells-1] = -1.; - /// \endcode /// \page tutorial4 /// \details We set up necessary information for the wells @@ -259,34 +135,27 @@ int main () /// \endcode /// \page tutorial4 - /// \details We set up the source term for the transport solver. + /// \details We set up the source term. Positive numbers indicate that the cell is a source, + /// while negative numbers indicate a sink. /// \code - TransportSource* tsrc = create_transport_source(2, 2); - double ssrc[] = { 1.0, 0.0 }; - double ssink[] = { 0.0, 1.0 }; - double zdummy[] = { 0.0, 0.0 }; - 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); - } - } + std::vector src(num_cells, 0.0); + src[0] = 1.; + src[num_cells-1] = -1.; /// \endcode /// \page tutorial4 /// \details We compute the pore volume /// \code std::vector porevol; - Opm::computePorevolume(*grid.c_grid(), props.porosity(), porevol); + Opm::computePorevolume(grid, props.porosity(), porevol); /// \endcode /// \page tutorial4 - /// \details We set up the transport solver. + /// \details Set up the transport solver. This is a reordering implicit Euler transport solver. /// \code - const bool guess_old_solution = true; - TransportModel model (fluid, *grid.c_grid(), porevol, grav, guess_old_solution); - TransportSolver tsolver(model); + const double tolerance = 1e-9; + const int max_iterations = 30; + Opm::TransportModelTwophase transport_solver(grid, props, tolerance, max_iterations); /// \endcode /// \page tutorial4 @@ -296,14 +165,6 @@ int main () int num_time_steps = 20; /// \endcode - /// \page tutorial4 - /// \details Control paramaters for the implicit solver. - /// \code - ImplicitTransportDetails::NRReport rpt; - ImplicitTransportDetails::NRControl ctrl; - /// \endcode - - /// \page tutorial4 /// \details We define a vector which contains all cell indexes. We use this @@ -321,21 +182,13 @@ int main () FlowBCManager bcs; /// \endcode - /// \page tutorial4 - /// \details - /// Linear solver init. - /// \code - using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; - CSRMatrixUmfpackSolver linsolve; - /// \endcode - /// \page tutorial4 /// \details /// We set up a two-phase state object, and /// initialise water saturation to minimum everywhere. /// \code TwophaseState state; - state.init(*grid.c_grid(), 2); + state.init(grid, 2); state.setFirstSat(allcells, props, TwophaseState::MinSat); /// \endcode @@ -450,7 +303,7 @@ int main () /// last argument. /// \code LinearSolverUmfpack linsolver; - IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver, wells); + IncompTpfa psolver(grid, props.permeability(), grav, linsolver, wells); /// \endcode @@ -469,7 +322,7 @@ int main () /// \page tutorial4 /// \details In order to use the well controls, we need to generate the WDP for each well. /// \code - Opm::computeWDP(*wells, *grid.c_grid(), state.saturation(), props.density(), gravity, true, wdp); + Opm::computeWDP(*wells, grid, state.saturation(), props.density(), gravity, true, wdp); /// \endcode /// \page tutorial4 @@ -511,7 +364,7 @@ int main () /// \page tutorial4 /// \details Transport solver /// \code - tsolver.solve(*grid.c_grid(), tsrc, dt, ctrl, state, linsolve, rpt); + transport_solver.solve(&state.faceflux()[0], &porevol[0], &src[0], dt, state.saturation()); /// \endcode /// \page tutorial4 @@ -523,75 +376,11 @@ int main () Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - Opm::writeVtkData(*grid.c_grid(), dm, vtkfile); + Opm::writeVtkData(grid, dm, vtkfile); } } /// \endcode -/// \page tutorial4 -/// \details Implementation of the TwophaseFluid class. -/// \code -TwophaseFluid::TwophaseFluid(const Opm::IncompPropertiesInterface& props) - : props_(props), - smin_(props.numCells()*props.numPhases()), - smax_(props.numCells()*props.numPhases()) -{ - const int num_cells = props.numCells(); - std::vector cells(num_cells); - for (int c = 0; c < num_cells; ++c) { - cells[c] = c; - } - props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]); - } - -double TwophaseFluid::density(int phase) const -{ - return props_.density()[phase]; -} - -template -void TwophaseFluid::mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const -{ - props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]); - const double* mu = props_.viscosity(); - mob[0] /= mu[0]; - mob[1] /= mu[1]; - /// \endcode - /// \page tutorial4 - /// \details We - /// recall that we use Fortran ordering for kr derivatives, - /// therefore dmob[i*2 + j] is row j and column i of the - /// matrix. - /// Each row corresponds to a kr function, so which mu to - /// divide by also depends on the row, j. - /// \code - dmob[0*2 + 0] /= mu[0]; - dmob[0*2 + 1] /= mu[1]; - dmob[1*2 + 0] /= mu[0]; - dmob[1*2 + 1] /= mu[1]; -} - -template -void TwophaseFluid::pc(int /*c */, const Sat& /* s*/, Pcap& pcap, DPcap& dpcap) const -{ - pcap = 0.; - dpcap = 0.; -} - -double TwophaseFluid::s_min(int c) const -{ - return smin_[2*c + 0]; -} - -double TwophaseFluid::s_max(int c) const -{ - return smax_[2*c + 0]; -} -/// \endcode /// \page tutorial4