From 12aaf935358d5732b15df9a975f6dcee3e6a8224 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 17 Apr 2012 17:24:19 +0200 Subject: [PATCH] Updated tutorials. --- tutorials/tutorial1.cpp | 3 +- tutorials/tutorial2.cpp | 4 +- tutorials/tutorial3.cpp | 387 ++++++++++++++++++++++++++++++++-------- 3 files changed, 316 insertions(+), 78 deletions(-) diff --git a/tutorials/tutorial1.cpp b/tutorials/tutorial1.cpp index 567a594c8..080244fdd 100644 --- a/tutorials/tutorial1.cpp +++ b/tutorials/tutorial1.cpp @@ -39,6 +39,7 @@ // ----------------- Main program ----------------- /// \page tutorial1 +/// \section commentedsource1 Commented source code: /// \code int main() { @@ -85,6 +86,6 @@ int main() /// \image html tutorial1.png /// \page tutorial1 -/// \section sourcecode Source code. +/// \section completecode1 Complete source code: /// \include tutorial1.cpp diff --git a/tutorials/tutorial2.cpp b/tutorials/tutorial2.cpp index ce3d12337..c34d96585 100644 --- a/tutorials/tutorial2.cpp +++ b/tutorials/tutorial2.cpp @@ -46,7 +46,7 @@ #include /// \page tutorial2 -/// \section commentedcode Commented code: +/// \section commentedcode2 Commented code: /// \code int main() { @@ -180,5 +180,5 @@ int main() /// \page tutorial2 -/// \section sourcecode Complete source code. +/// \section completecode2 Complete source code: /// \include tutorial2.cpp diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp index 76d2cbdbb..85e8ac236 100644 --- a/tutorials/tutorial3.cpp +++ b/tutorials/tutorial3.cpp @@ -51,83 +51,129 @@ #include /// \page tutorial3 Multiphase flow -/// Multiphase flow +/// The Darcy law gives +/// \f[u_\alpha= -\frac1{\mu_\alpha} K_\alpha\nabla p_\alpha\f] +/// where \f$\mu_\alpha\f$ and \f$K_\alpha\f$ represent the viscosity +/// and the permeability tensor for each phase \f$\alpha\f$. In the two phase +/// case, we have either \f$\alpha=w\f$ or \f$\alpha=o\f$. +/// In this tutorial, we do not take into account capillary pressure so that +/// \f$p=p_w=p_o\f$ and gravity +/// effects. We denote by \f$K\f$ the absolute permeability tensor and each phase +/// permeability is defined through its relative permeability by the expression +/// \f[K_\alpha=k_{r\alpha}K.\f] +/// The phase mobility are defined as +/// \f[\lambda_\alpha=\frac{k_{r\alpha}}{\mu_\alpha}\f] +/// so that the Darcy law may be rewritten as +/// \f[u_\alpha= -\lambda_\alpha K\nabla p.\f] +/// The conservation of mass for each phase writes: +/// \f[\frac{\partial}{\partial t}(\phi\rho_\alpha s_\alpha)+\nabla\cdot(\rho_\alpha u_\alpha)=q_\alpha\f] +/// where \f$s_\alpha\f$ denotes the saturation of the phase \f$\alpha\f$ and \f$q_\alpha\f$ is a source term. Let +/// us consider a two phase flow with oil and water. We assume that the phases are incompressible. Since +/// \f$s_w+s_o=1\f$, we get +/// \f[ \nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o}.\f] +/// where we define +/// \f[u=u_w+u_o.\f] +/// Let the total mobility be equal to +/// \f[\lambda=\lambda_w+\lambda_o\f] +/// Then, we have +/// \f[u=-\lambda K\nabla p.\f] +/// The set of equations +/// \f[\nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o},\quad u=-\lambda K\nabla p.\f] +/// is referred to as the pressure equation. We introduce +/// the fractional flow \f$f_w\f$ +/// as +/// \f[f_w=\frac{\lambda_w}{\lambda_w+\lambda_o}\f] +/// and obtain +/// \f[\phi\frac{\partial s}{\partial t}+\nabla\cdot(f_w u)=\frac{q_w}{\rho_w}\f] +/// which is referred to as the transport equation. The pressure and +/// transport equation are coupled. In this tutorial, we implement a splitting scheme, +/// where, at each time step, we decouple the two equations. We solve first +/// the pressure equation and then update the water saturation by solving +/// the transport equation assuming that \f$u\f$ is constant in time in the time step +/// interval we are considering. +/// \page tutorial3 +/// \section commentedcode3 Commented code: +/// \page tutorial3 +/// \details +/// We define a class which computes mobility, capillary pressure and +/// the minimum and maximum saturation value for each cell. +/// \code class TwophaseFluid { public: - 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]); - } + /// \endcode + /// \page tutorial3 + /// \details Constructor operator. Takes in the fluid properties defined + /// \c props + /// \code + TwophaseFluid(const Opm::IncompPropertiesInterface& props); + /// \endcode - double density(int phase) const - { - return props_.density()[phase]; - } + /// \page tutorial3 + /// \details Density for each phase. + /// \code + double density(int phase) const; + /// \endcode - template - void 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]; - // 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. - dmob[0*2 + 0] /= mu[0]; - dmob[0*2 + 1] /= mu[1]; - dmob[1*2 + 0] /= mu[0]; - dmob[1*2 + 1] /= mu[1]; - } + /// \page tutorial3 + /// \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 - template - void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const - { - double pcow[2]; - double dpcow[4]; - props_.capPress(1, &s[0], &c, pcow, dpcow); - pcap = pcow[0]; - dpcap = dpcow[0]; - } + /// \page tutorial3 + /// \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 - double s_min(int c) const - { - return smin_[2*c + 0]; - } + /// \page tutorial3 + /// \details Returns the minimum permitted saturation. + /// \code + double s_min(int c) const; + /// \endcode - double s_max(int c) const - { - return smax_[2*c + 0]; - } + /// \page tutorial3 + /// \details Returns the maximum permitted saturation + /// \code + double s_max(int c) const; + /// \endcode + /// \page tutorial3 + /// \details Private variables + /// \code private: const Opm::IncompPropertiesInterface& props_; std::vector smin_; std::vector smax_; }; +/// \endcode +/// \page tutorial3 +/// \details We set up the transport model. +/// \code typedef Opm::SinglePointUpwindTwoPhase TransportModel; +/// \endcode +/// \page tutorial3 +/// \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; +typedef NewtonVectorCollection< ::std::vector > NVecColl; +typedef JacobianSystem< struct CSRMatrix, NVecColl > JacSys; template class MaxNorm { @@ -137,7 +183,12 @@ public: return AccumulationNorm ::norm(v); } }; +/// \endcode +/// \page tutorial3 +/// \details +/// We set up the solver. +/// \code typedef Opm::ImplicitTransport TransportSolver; - +/// \endcode +/// \page tutorial3 +/// \details +/// Main function +/// \code int main () { - + /// \endcode + /// \page tutorial3 + /// \details + /// We define the grid. A cartesian grid with 1200 cells. + /// \code int dim = 3; int nx = 20; int ny = 20; @@ -162,41 +221,88 @@ int main () GridManager grid(nx, ny, nz, dx, dy, dz); int num_cells = grid.c_grid()->number_of_cells; int num_faces = grid.c_grid()->number_of_faces; + /// \endcode + /// \page tutorial3 + /// \details + /// We define the properties of the fluid.\n + /// Number of phases. + /// \code int num_phases = 2; using namespace unit; using namespace prefix; + /// \endcode + /// \page tutorial3 + /// \details density vector (one component per phase). + /// \code std::vector rho(2, 1000.); + /// \endcode + /// \page tutorial3 + /// \details viscosity vector (one component per phase). + /// \code std::vector mu(2, 1.*centi*Poise); + /// \endcode + /// \page tutorial3 + /// \details porosity and permeability of the rock. + /// \code double porosity = 0.5; double k = 10*milli*darcy; - + /// \endcode + + /// \page tutorial3 + /// \details We define the relative permeability function. We use a basic fluid + /// 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; + /// \endcode + /// \page tutorial3 + /// \details We construct a basic fluid with the properties we have defined above. + /// Each property is constant and hold for all cells. + /// \code IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu, porosity, k, dim, num_cells); + TwophaseFluid fluid(props); + /// \endcode - - + + /// \page tutorial3 + /// \details Gravity parameters. Here, we set zero gravity. + /// \code const double *grav = 0; std::vector omega; + /// \endcode + /// \page tutorial3 + /// \details We set up the pressure solver. + /// \code LinearSolverUmfpack linsolver; IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver); + /// \endcode + - TwophaseFluid fluid(props); - + /// \page tutorial3 + /// \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 tutorial3 + /// \details We set up the wells. Here, there are no well and we let them empty. + /// \code std::vector empty_wdp; std::vector empty_well_bhp; std::vector empty_well_flux; + /// \endcode + /// \page tutorial3 + /// \details We set up the source term for the transport solver. + /// \code TransportSource* tsrc = create_transport_source(2, 2); double ssrc[] = { 1.0, 0.0 }; double ssink[] = { 0.0, 1.0 }; @@ -208,48 +314,112 @@ int main () append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc); } } - + /// \endcode + + /// \page tutorial3 + /// \details We compute the pore volume + /// \code std::vector porevol; computePorevolume(*grid.c_grid(), props, porevol); + /// \endcode + /// \page tutorial3 + /// \details We set up the transport solver. + /// \code const bool guess_old_solution = true; TransportModel model (fluid, *grid.c_grid(), porevol, grav, guess_old_solution); TransportSolver tsolver(model); + /// \endcode + /// \page tutorial3 + /// \details Time integration parameters + /// \code double dt = 0.1*day; int num_time_steps = 20; + /// \endcode - TwophaseState state; + /// \page tutorial3 + /// \details Control paramaters for the implicit solver. + /// \code ImplicitTransportDetails::NRReport rpt; ImplicitTransportDetails::NRControl ctrl; - std::vector totmob; + /// \endcode + + + /// \page tutorial3 + /// \details We define a vector which contains all cell indexes. We use this + /// vector to set up parameters on the whole domains. std::vector allcells(num_cells); for (int cell = 0; cell < num_cells; ++cell) { allcells[cell] = cell; } - FlowBCManager bcs; - // Linear solver init. + /// \page tutorial3 + /// \details We set up the boundary conditions. Letting bcs empty is equivalent + /// to no flow boundary conditions. + /// \code + FlowBCManager bcs; + /// \endcode + + /// \page tutorial3 + /// \details + /// Linear solver init. + /// \code using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; CSRMatrixUmfpackSolver linsolve; + /// \endcode + /// \page tutorial3 + /// \details + /// We initialise water saturation to minimum everywhere. + /// \code + TwophaseState state; state.init(*grid.c_grid()); - // By default: initialise water saturation to minimum everywhere. state.setWaterSat(allcells, props, TwophaseState::MinSat); + /// \endcode + /// \page tutorial3 + /// \details We introduce a vector which contains the total mobility + /// on all cells. + /// \code + std::vector totmob; + /// \endcode + /// \page tutorial3 + /// \details This string will contain the name of a VTK output vector. + /// \code std::ostringstream vtkfilename; + /// \endcode + + /// \page tutorial3 + /// \details Loop over the time steps. + /// \code for (int i = 0; i < num_time_steps; ++i) { - + /// \endcode + /// \page tutorial3 + /// \details Compute the total mobility. It is needed by the pressure solver + /// \code computeTotalMobility(props, allcells, state.saturation(), totmob); + /// \endcode + /// \page tutorial3 + /// \details Solve the pressure equation + /// \code psolver.solve(totmob, omega, src, empty_wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), empty_well_bhp, empty_well_flux); + /// \endcode + /// \page tutorial3 + /// \details Transport solver + /// \code tsolver.solve(*grid.c_grid(), tsrc, dt, ctrl, state, linsolve, rpt); + /// \endcode + /// \page tutorial3 + /// \details Write the output to file. + /// \code vtkfilename.str(""); vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); @@ -257,12 +427,78 @@ int main () dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); Opm::writeVtkData(*grid.c_grid(), dm, vtkfile); - } } +/// \endcode + +/// \page tutorial3 +/// \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 tutorial3 + /// \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 tutorial3 +/// \section results3 Results. /// /// /// @@ -279,5 +515,6 @@ int main () ///
\image html tutorial3-000.png
/// \page tutorial3 -/// \section sourcecode Complete source code. +/// \details +/// \section completecode3 Complete source code: /// \include tutorial3.cpp