diff --git a/examples/Makefile.am b/examples/Makefile.am new file mode 100644 index 00000000..01b8a5e0 --- /dev/null +++ b/examples/Makefile.am @@ -0,0 +1,29 @@ +# Additional definitions to ensure proper include paths + +AM_CPPFLAGS = \ +-I$(top_srcdir)/src + +AM_LDFLAGS = \ +$(LAPACK_LIBS) $(BLAS_LIBS) + +#----------------------------------------------------------------------- + +# Declare example programs +noinst_PROGRAMS = + +# "spu_2p" depends on availability of OPMPressure and UMFPACK +if UMFPACK +if OPMPRESSURE_SPARSE_SYS +noinst_PROGRAMS += spu_2p +endif +endif + +#----------------------------------------------------------------------- + +# Example program sources + +spu_2p_SOURCES = \ +spu_2p.cpp \ +call_umfpack.c \ +cart_grid.c \ +transport_source.c diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp new file mode 100644 index 00000000..ae1ead97 --- /dev/null +++ b/examples/spu_2p.cpp @@ -0,0 +1,289 @@ +/*=========================================================================== +// +// File: spu_2p.cpp +// +// Created: 2011-10-05 10:29:01+0200 +// +// Authors: Ingeborg S. Ligaarden +// Jostein R. Natvig +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2011 SINTEF ICT, Applied Mathematics. + Copyright 2011 Statoil ASA. + + 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 +#include + +#include +#include +#include + +#include + +#include + +#include + +#include +#include +#include +#include + +#include + +#include +#include + +template +Ostream& +operator<<(Ostream& os, const Collection& c) +{ + typedef typename Collection::value_type VT; + + os << "[ "; + std::copy(c.begin(), c.end(), ::std::ostream_iterator(os, " ")); + os << "]"; + + return os; +} + +class Rock { +public: + Rock(::std::size_t nc, ::std::size_t dim) + : dim_ (dim ), + perm_(nc * dim * dim), + poro_(nc ) {} + + const ::std::vector& perm() const { return perm_; } + const ::std::vector& poro() const { return poro_; } + + void + perm_homogeneous(double k) { + setVector(0.0, perm_); + + const ::std::size_t d2 = dim_ * dim_; + + for (::std::size_t c = 0, nc = poro_.size(); c < nc; ++c) { + for (::std::size_t i = 0; i < dim_; ++i) { + perm_[c*d2 + i*(dim_ + 1)] = k; + } + } + } + + void + poro_homogeneous(double phi) { + setVector(phi, poro_); + } + +private: + void + setVector(double x, ::std::vector& v) { + ::std::fill(v.begin(), v.end(), x); + } + + ::std::size_t dim_ ; + ::std::vector perm_; + ::std::vector poro_; +}; + +template +class ReservoirState { +public: + ReservoirState(const grid_t* g) + : press_ (g->number_of_cells), + fpress_(g->number_of_faces), + flux_ (g->number_of_faces), + sat_ (np * g->number_of_cells) + {} + + ::std::vector& pressure () { return press_ ; } + ::std::vector& facepressure() { return fpress_; } + ::std::vector& faceflux () { return flux_ ; } + ::std::vector& saturation () { return sat_ ; } + + const ::std::vector& faceflux () const { return flux_; } + const ::std::vector& saturation () const { return sat_ ; } + +private: + ::std::vector press_ ; + ::std::vector fpress_; + ::std::vector flux_ ; + ::std::vector sat_ ; +}; + +class PressureSolver { +public: + PressureSolver(grid_t* g, const Rock& rock) + : htrans_(g->cell_facepos[ g->number_of_cells ]), + trans_ (g->number_of_faces), + gpress_(g->cell_facepos[ g->number_of_cells ]) + { + tpfa_htrans_compute(g, &rock.perm()[0], &htrans_[0]); + + h_ = ifs_tpfa_construct(g); + } + + ~PressureSolver() { + ifs_tpfa_destroy(h_); + } + + template + void + solve(grid_t* g , + const ::std::vector& totmob, + const ::std::vector& src , + State& state ) { + + tpfa_eff_trans_compute(g, &totmob[0], &htrans_[0], &trans_[0]); + + // No gravity + ::std::fill(gpress_.begin(), gpress_.end(), double(0.0)); + + ifs_tpfa_assemble(g, &trans_[0], &src[0], &gpress_[0], h_); + + using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; + + CSRMatrixUmfpackSolver linsolve; + linsolve.solve(h_->A, h_->b, h_->x); + + ifs_tpfa_press_flux(g, &trans_[0], h_, + &state.pressure()[0], + &state.faceflux()[0]); + } + +private: + ::std::vector htrans_; + ::std::vector trans_ ; + ::std::vector gpress_; + + struct ifs_tpfa_data* h_; +}; + + +typedef Opm::SimpleFluid2p<> TwophaseFluid; +typedef Opm::SinglePointUpwindTwoPhase TransportModel; + +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); + } +}; + +typedef Opm::ImplicitTransport TransportSolver; + +void +compute_porevolume(const grid_t* g, + const Rock& rock, + std::vector& porevol) +{ + const ::std::vector& poro = rock.poro(); + +#if 0 + assert (poro.size() == static_cast<::std::size_t>(g->number_of_cells)); +#endif + + porevol.resize(rock.poro().size()); + + ::std::transform(poro.begin(), poro.end(), + g->cell_volumes, + porevol.begin(), + ::std::multiplies()); +} + +int +main() +{ + grid_t* grid = create_cart_grid(2, 1, 1); + + Rock rock(grid->number_of_cells, grid->dimensions); + + rock.perm_homogeneous(1); + rock.poro_homogeneous(1); + + PressureSolver psolver(grid, rock); + + std::vector totmob(grid->number_of_cells, 1.0); + std::vector src (grid->number_of_cells, 0.0); + + src[0] = 1.0; + src[grid->number_of_cells - 1] = -1.0; + + ReservoirState<> state(grid); + + psolver.solve(grid, totmob, src, state); + + 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 }; + + append_transport_source(0, 2, 0, src[0], ssrc, zdummy, tsrc); + append_transport_source(grid->number_of_cells - 1, 2, 0, + src.back(), ssink, zdummy, tsrc); + + Opm::ImplicitTransportDetails::NRReport rpt; + Opm::ImplicitTransportDetails::NRControl ctrl; + ctrl.max_it = 1; + + using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; + CSRMatrixUmfpackSolver linsolve; + + std::array mu = {{ 1.0, 1.0 }}; + std::array rho = {{ 0.0, 0.0 }}; + TwophaseFluid fluid(mu, rho); + + std::vector porevol; + compute_porevolume(grid, rock, porevol); + + TransportModel model (fluid, *grid, porevol); + TransportSolver tsolver(model); + + double dt = 1; + tsolver.solve(*grid, tsrc, dt, ctrl, state, linsolve, rpt); + + std::cerr << state.saturation() << '\n'; + + destroy_transport_source(tsrc); + destroy_cart_grid(grid); +}