diff --git a/.hgignore b/.hgignore index f5b993d..de6c326 100644 --- a/.hgignore +++ b/.hgignore @@ -34,3 +34,4 @@ missing *_test dune/porsol/blackoil/test/bo_fluid_pressuredeps dune/porsol/blackoil/test/bo_fluid_p_and_z_deps +examples/spu_2p diff --git a/examples/Makefile.am b/examples/Makefile.am index 8dcf9dc..3336c77 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,23 +1,24 @@ SUBDIRS = check_PROGRAMS = -noinst_PROGRAMS = \ -aniso_implicitcap_test \ -aniso_simulator_test \ -blackoil_sim_test \ -grdecl_to_legacy_test \ -ifsh_solver_test \ -implicitcap_test \ -istl_test \ -known_answer_test \ -make_vtk_test \ -mimetic_aniso_solver_test \ -mimetic_periodic_test \ -mimetic_solver_test \ -simulator_test \ -spe10_test \ -tpfa_compressible_solver_test \ -tpfa_solver_test \ +noinst_PROGRAMS = \ +aniso_implicitcap_test \ +aniso_simulator_test \ +blackoil_sim_test \ +grdecl_to_legacy_test \ +ifsh_solver_test \ +implicitcap_test \ +istl_test \ +known_answer_test \ +make_vtk_test \ +mimetic_aniso_solver_test \ +mimetic_periodic_test \ +mimetic_solver_test \ +simulator_test \ +spe10_test \ +spu_2p \ +tpfa_compressible_solver_test \ +tpfa_solver_test \ twophase2_test @@ -36,6 +37,9 @@ mimetic_periodic_test_SOURCES = mimetic_periodic_test.cpp istl_test_SOURCES = istl_test.cpp +spu_2p_SOURCES = spu_2p.cpp +spu_2p_LDADD = $(LDADD) ../lib/libduneporsol.la + spe10_test_SOURCES = spe10_test.cpp tpfa_solver_test_SOURCES = tpfa_solver_test.cpp diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp new file mode 100644 index 0000000..4dbd350 --- /dev/null +++ b/examples/spu_2p.cpp @@ -0,0 +1,407 @@ +/*=========================================================================== +// +// File: spu_2p.cpp +// +// Created: 2011-10-10 10:35:13+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 + +#include +#include +#include + +#include + +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 PressureLinearSolver { +public: + PressureLinearSolver() + { + Dune::parameter::ParameterGroup params; + + params.insertParameter("linsolver_tolerance", + boost::lexical_cast(5.0e-9)); + + params.insertParameter("linsoler_verbosity", + boost::lexical_cast(1)); + + params.insertParameter("linsolver_type", + boost::lexical_cast(1)); + + ls_.init(params); + } + + void + solve(struct CSRMatrix* A, + const double* b, + double* x) + { + Dune::LinearSolverISTL::LinearSolverResults res = + ls_.solve(A->m, A->nnz, A->ia, A->ja, A->sa, b, x); + } + +private: + Dune::LinearSolverISTL ls_; +}; + +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_); + + PressureLinearSolver 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_; +}; + +class TwophaseFluid { +public: + TwophaseFluid(const Dune::ReservoirPropertyCapillary<3>& r) + : r_(r) + {} + + template + void + mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const { + const double s1 = s[0]; + + r_.phaseMobilities (c, s1, mob ); + r_.phaseMobilitiesDeriv(c, s1, dmob); + } + + double density(int p) const { + if (p == 0) { + return r_.densityFirstPhase(); + } else { + return r_.densitySecondPhase(); + } + } + + +private: + const Dune::ReservoirPropertyCapillary<3>& r_; +}; + +typedef Opm::SinglePointUpwindTwoPhase TransportModel; + +using namespace Opm::ImplicitTransportDefault; + +typedef Dune::FieldVector ScalarVectorBlockType; +typedef Dune::FieldMatrix ScalarMatrixBlockType; + +typedef Dune::BlockVector ScalarBlockVector; +typedef Dune::BCRSMatrix ScalarBCRSMatrix; + +typedef NewtonVectorCollection< ScalarBlockVector > NVecColl; +typedef JacobianSystem < ScalarBCRSMatrix, NVecColl > JacSys; + +class TransportLinearSolver { +public: + void + solve(const ScalarBCRSMatrix& A, + const ScalarBlockVector& b, + ScalarBlockVector& x) { + + Dune::MatrixAdapter opA(A); + + Dune::SeqILU0 precond(A, 1.0); + + int maxit = A.N(); + double tol = 5.0e-7; + int verb = 1; + + Dune::BiCGSTABSolver + solver(opA, precond, tol, maxit, verb); + + ScalarBlockVector bcpy(b); + Dune::InverseOperatorResult res; + solver.apply(x, bcpy, res); + } +}; + +template +class MaxNorm { +public: + static double + norm(const Vector& v) { + return v.infinity_norm(); + } +}; + +typedef Opm::ImplicitTransport TransportSolver; + +void +compute_porevolume(const grid_t* g, + const Rock& rock, + std::vector& porevol) +{ + const ::std::vector& poro = rock.poro(); + + assert (poro.size() == (::std::size_t)(g->number_of_cells)); + + porevol.resize(rock.poro().size()); + + ::std::transform(poro.begin(), poro.end(), + g->cell_volumes, + porevol.begin(), + ::std::multiplies()); +} + +class TransportSource { +public: + TransportSource() : nsrc(0) {} + + int nsrc ; + ::std::vector< int > cell ; + ::std::vector pressure ; + ::std::vector flux ; + ::std::vector saturation; +}; + +template +void +append_transport_source(int c, double p, double v, const Arr& s, + TransportSource& src) +{ + src.cell .push_back(c); + src.pressure .push_back(p); + src.flux .push_back(v); + src.saturation.insert(src.saturation.end(), + s.begin(), s.end()); + ++src.nsrc; +} + +int +main(int argc, char** argv) +{ + Dune::parameter::ParameterGroup param(argc, argv); + + Dune::CpGrid cp_grid; + Dune::ReservoirPropertyCapillary<3> res_prop; + + setupGridAndProps(param, cp_grid, res_prop); + + res_prop.init(cp_grid.size(0), 1, 1); + res_prop.setViscosities(1.0, 1.0); + res_prop.setDensities (0.0, 0.0); + + GridAdapter grid; + grid.init(cp_grid); + + Rock rock(grid.c_grid()->number_of_cells, grid.c_grid()->dimensions); + + rock.perm_homogeneous(1); + rock.poro_homogeneous(1); + + PressureSolver psolver(grid.c_grid(), rock); + + std::vector totmob(grid.c_grid()->number_of_cells, 1.0); + std::vector src (grid.c_grid()->number_of_cells, 0.0); + + src[0] = 1.0; + src[grid.c_grid()->number_of_cells - 1] = -1.0; + + ReservoirState<> state(grid.c_grid()); + + psolver.solve(grid.c_grid(), totmob, src, state); + + TransportSource tsrc; + ::std::vector ssrc (2, 0.0); ssrc[0] = 1.0; + ::std::vector ssink(2, 0.0); + append_transport_source(0, 1.0, src[0] , ssrc , tsrc); + append_transport_source(grid.c_grid()->number_of_cells - 1, + 1.0, src.back(), ssink, tsrc); + + Opm::ImplicitTransportDetails::NRReport rpt; + Opm::ImplicitTransportDetails::NRControl ctrl; + + std::vector porevol; + compute_porevolume(grid.c_grid(), rock, porevol); + + TwophaseFluid fluid (res_prop); + TransportModel model (fluid, *grid.c_grid(), porevol); + TransportSolver tsolver(model); + + double dt = 1e2; + ctrl.max_it = 20 ; + + TransportLinearSolver linsolve; + tsolver.solve(*grid.c_grid(), &tsrc, dt, ctrl, state, linsolve, rpt); + + std::cerr << "Number of linear solves: " << rpt.nit << '\n' + << "Process converged: " << (rpt.flag > 0) << '\n' + << "Convergence flag: " << rpt.flag << '\n' + << "Final residual norm: " << rpt.norm_res << '\n' + << "Final increment norm: " << rpt.norm_dx << '\n'; + + ::std::ofstream sfile("saturation-00.txt"); + sfile.setf(::std::ios::showpos | ::std::ios::scientific); + sfile.precision(15); + ::std::copy(state.saturation().begin(), + state.saturation().end (), + ::std::ostream_iterator(sfile, "\n")); +}