diff --git a/.hgsubstate b/.hgsubstate index e2142a6..abb19d2 100644 --- a/.hgsubstate +++ b/.hgsubstate @@ -1,2 +1,2 @@ -1589895679c9f47cafdbf98ebccd6a9afe6644d2 dune/porsol/opmpressure -b57c13fcc23003a3e04dfece2900f897c2fe5e99 dune/porsol/opmtransport +07ae49ffab2c469150c67856880b120f8e96930b dune/porsol/opmpressure +888d175acc5817a6206c55132427307e9b869598 dune/porsol/opmtransport diff --git a/dune/porsol/euler/EulerUpstreamImplicit.hpp b/dune/porsol/euler/EulerUpstreamImplicit.hpp index e3bdbe7..76a681d 100644 --- a/dune/porsol/euler/EulerUpstreamImplicit.hpp +++ b/dune/porsol/euler/EulerUpstreamImplicit.hpp @@ -51,8 +51,8 @@ along with OpenRS. If not, see . -#include - +//#include +#include #include #include @@ -141,28 +141,39 @@ namespace Dune { // try to make bard ingerfaces - typedef Opm::SimpleFluid2pWrapper TwophaseFluid; + //typedef Opm::SimpleFluid2pWrapper TwophaseFluid; + typedef TwophaseFluidWrapper TwophaseFluid; typedef Opm::SinglePointUpwindTwoPhase TransportModel; // using namespace Opm::ImplicitTransportDefault - typedef Opm::ImplicitTransportDefault::NewtonVectorCollection< ::std::vector > NVecColl; - typedef Opm::ImplicitTransportDefault::JacobianSystem < struct CSRMatrix, NVecColl > JacSys; - + //typedef Opm::ImplicitTransportDefault::NewtonVectorCollection< ::std::vector > NVecColl; + //typedef Opm::ImplicitTransportDefault::JacobianSystem < struct CSRMatrix, NVecColl > JacSys; + //using namespace Opm::ImplicitTransportDefault; + typedef Dune::FieldVector ScalarVectorBlockType; + typedef Dune::FieldMatrix ScalarMatrixBlockType; + + typedef Dune::BlockVector ScalarBlockVector; + typedef Dune::BCRSMatrix ScalarBCRSMatrix; + + typedef Opm::ImplicitTransportDefault::NewtonVectorCollection< ScalarBlockVector > NVecColl; + typedef Opm::ImplicitTransportDefault::JacobianSystem < ScalarBCRSMatrix, NVecColl > JacSys; + typedef Opm::ImplicitTransport TransportSolver; - // should be initialized by param - Opm::ImplicitTransportDetails::NRReport rpt_; - Opm::ImplicitTransportDetails::NRControl ctrl_; - Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver linsolve_; - - TransportModel model_; - TransportSolver tsolver_; - GridAdapter mygrid_; - Opm::SimpleFluid2pWrapper< ReservoirProperties > myfluid_; + // should be initialized by param + + + //TransportModel model_; + //TransportSolver tsolver_; + GridAdapter mygrid_; + ReservoirProperties myrp_; + + //Opm::SimpleFluid2pWrapper< ReservoirProperties > myfluid_; + //TwophaseFluid myfluid_; bool check_sat_; bool clamp_sat_; diff --git a/dune/porsol/euler/EulerUpstreamImplicit_impl.hpp b/dune/porsol/euler/EulerUpstreamImplicit_impl.hpp index 9e07f88..07d15a0 100644 --- a/dune/porsol/euler/EulerUpstreamImplicit_impl.hpp +++ b/dune/porsol/euler/EulerUpstreamImplicit_impl.hpp @@ -51,6 +51,7 @@ //#include #include #include +#include #include #include @@ -73,7 +74,7 @@ namespace Dune //residual_computer_.initObj(g, r, b); //std::array mu = {{ 1.0, 1.0 }}; //std::array rho = {{ 0.0, 0.0 }}; - myfluid_ = Opm::SimpleFluid2pWrapper(r); + //myfluid_ = Opm::SimpleFluid2pWrapper(r); check_sat_ = true; clamp_sat_ = false; //residual_computer_.initObj(g, r, b); @@ -109,12 +110,13 @@ namespace Dune porevol_[i]= mygrid_.cellVolume(i)*r.porosity(i); } - std::array< double, 3 > gravity; - std::vector< double > trans; - trans.resize(mygrid_.numFaces()); - model_ = TransportModel(myfluid_,mygrid_.c_grid(),porevol_,&gravity[0],&trans[0]); + + myrp_= r; + //model_ = TransportModel(myfluid_,mygrid_.c_grid(),porevol_,&gravity[0],&trans[0]); + //myfluid_.init(r); + //model_.init(myfluid_,mygrid_.c_grid(),porevol_,&gravity[0],&trans[0]); //model.init(myfluid_,mygrid_.c_grid(), &porevol_, gravity, &trans); - tsolver_ = TransportSolver(model_); + //tsolver_ = TransportSolver(model_); } @@ -145,22 +147,21 @@ namespace Dune sat[2*i+1] = 1-saturation[i]; } - typedef typename GI::CellIterator CIt; - typedef typename CIt::FaceIterator FIt; - typedef typename FIt::Vector Vector; - int count=0; - grid_t* cgrid = mygrid_.c_grid(); + + //int count=0; + const grid_t* cgrid = mygrid_.c_grid(); int numhf = cgrid->cell_facepos[cgrid->number_of_cells]; //for (int i=0; i < numhf); ++i){ - + std::vector faceflux(numhf); for (int c = 0, i = 0; c < cgrid->number_of_cells; ++c){ for (; i < cgrid->cell_facepos[c + 1]; ++i) { int f= cgrid->cell_faces[i]; double outflux = pressure_sol.outflux(i); double sgn = 2.0*(cgrid->face_cells[2*f + 0] == c) - 1; - faceflux_[f] = sgn * outflux; + faceflux[f] = sgn * outflux; } } + state.faceflux()=faceflux; double dt_transport = time; int nr_transport_steps = 1; @@ -170,12 +171,26 @@ namespace Dune bool finished = false; clock.start(); std::vector< double > saturation_initial(saturation); + //std::array< double, 3 > mygravity; + std::vector< double > trans; + trans.resize(mygrid_.numFaces()); + TwophaseFluid myfluid(myrp_); + TransportModel model(myfluid,*mygrid_.c_grid(),porevol_,&gravity[0],&trans[0]); + TransportSolver tsolver(model); + Opm::ImplicitTransportDetails::NRReport rpt_; + Opm::ImplicitTransportDetails::NRControl ctrl_; + //Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver linsolve_; + TransportLinearSolver linsolve_; + //std::vector totmob(mygrid_.numCells(), 1.0); + //std::vector src (mygrid_.numCells(), 0.0); + //src[0] = 1.0; + // src[grid->number_of_cells - 1] = -1.0; + TransportSource* tsrc = 0;//create_transport_source(0, 2); while (!finished) { try { for (int q = 0; q < nr_transport_steps; ++q) { - //grid_t* pp=mygrid_.c_grid(); - //tsolver_.solve(*pp, NULL, dt_transport, ctrl_, state, linsolve_, rpt_); + tsolver.solve(*mygrid_.c_grid(), tsrc, dt_transport, ctrl_, state, linsolve_, rpt_); } #ifdef VERBOSE std::cout << "Doing " << 1 diff --git a/examples/Makefile.am b/examples/Makefile.am index 3336c77..ad1cd28 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -19,7 +19,8 @@ spe10_test \ spu_2p \ tpfa_compressible_solver_test \ tpfa_solver_test \ -twophase2_test +twophase2_test \ +spu_2p grdecl_to_legacy_test_SOURCES = grdecl_to_legacy_test.cpp diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 4dbd350..ed67253 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -212,6 +212,10 @@ public: TwophaseFluid(const Dune::ReservoirPropertyCapillary<3>& r) : r_(r) {} + void init(const Dune::ReservoirPropertyCapillary<3>& r) + { + r_ = r; + } template & r_; + Dune::ReservoirPropertyCapillary<3> r_; }; typedef Opm::SinglePointUpwindTwoPhase TransportModel;