Updated tutorial4 to use the reorder solver.

tutorials/tutorial4.cpp

@@ -35,14 +35,7 @@
 #include <opm/core/pressure/FlowBCManager.hpp>
 #include <opm/core/fluid/IncompPropertiesBasic.hpp>
 
-#include <opm/core/transport/transport_source.h>
+#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
-#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
-#include <opm/core/transport/NormSupport.hpp>
-#include <opm/core/transport/ImplicitAssembly.hpp>
-#include <opm/core/transport/ImplicitTransport.hpp>
-#include <opm/core/transport/JacobianSystem.hpp>
-#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
-#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
 
 #include <opm/core/simulator/TwophaseState.hpp>
 
@@ -52,113 +45,6 @@
 #include <opm/core/wells/WellCollection.hpp>
 
 /// \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 <class Sat, class Mob, class DMob>  
-    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 <class Sat, class Pcap, class DPcap>  
-    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<double> smin_;
-    std::vector<double> smax_;
-};
-/// \endcode
-
-/// \page tutorial4
-/// \details We set up the transport model.
-/// \code
-typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> 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<double> >  NVecColl;
-typedef JacobianSystem< struct CSRMatrix, NVecColl > JacSys;
-
-template <class Vector>
-class MaxNorm {
-public:
-    static double
-    norm(const Vector& v) {
-        return AccumulationNorm <Vector, MaxAbs>::norm(v);
-    }
-};
-/// \endcode
-
-/// \page tutorial4
-/// \details
-/// We set up the solver.
-/// \code
-typedef Opm::ImplicitTransport<TransportModel,
-        JacSys        ,
-        MaxNorm       ,
-        VectorNegater ,
-        VectorZero    ,
-        MatrixZero    ,
-        VectorAssign  > 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);
+    GridManager grid_manager(nx, ny, nz, dx, dy, dz);
-    int num_cells = grid.c_grid()->number_of_cells;
+    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;
+    SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
-    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<double> 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);
+    std::vector<double> src(num_cells, 0.0);
-    double ssrc[]   = { 1.0, 0.0 };
+    src[0] = 1.;
-    double ssink[]  = { 0.0, 1.0 };
+    src[num_cells-1] = -1.;
-    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);
-        }
-    }
     /// \endcode
 
     /// \page tutorial4
     /// \details We compute the pore volume
     /// \code
     std::vector<double> 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;
+    const double tolerance = 1e-9;
-    TransportModel  model  (fluid, *grid.c_grid(), porevol, grav, guess_old_solution);
+    const int max_iterations = 30;
-    TransportSolver tsolver(model);
+    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<int> 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 <class Sat,
-          class Mob,
-          class DMob>
-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 <class Sat,
-          class Pcap,
-          class DPcap>
-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