diff --git a/examples/sim_2p_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp index ebfa6de4..b47edd7c 100644 --- a/examples/sim_2p_incomp_reorder.cpp +++ b/examples/sim_2p_incomp_reorder.cpp @@ -43,7 +43,7 @@ #include #include -#include +#include #include #include @@ -200,7 +200,7 @@ main(int argc, char** argv) if (!use_deck) { // Simple simulation without a deck. WellsManager wells; // no wells. - SimulatorIncompTwophase simulator(param, + SimulatorIncompTwophaseReorder simulator(param, *grid->c_grid(), *props, rock_comp->isActive() ? rock_comp.get() : 0, @@ -255,7 +255,7 @@ main(int argc, char** argv) } // Create and run simulator. - SimulatorIncompTwophase simulator(param, + SimulatorIncompTwophaseReorder simulator(param, *grid->c_grid(), *props, rock_comp->isActive() ? rock_comp.get() : 0, diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 5a546045..3a131518 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -74,7 +74,7 @@ #include #include -#include +#include #include #include @@ -455,7 +455,7 @@ main(int argc, char** argv) // Reordering solver. const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); const int nl_maxiter = param.getDefault("nl_maxiter", 30); - Opm::TransportModelTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); + Opm::TransportSolverTwophaseReorder reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); if (use_gauss_seidel_gravity) { reorder_model.initGravity(grav); } diff --git a/opm/core/GridManager.cpp b/opm/core/GridManager.cpp index 76fa6fce..14ff6af3 100644 --- a/opm/core/GridManager.cpp +++ b/opm/core/GridManager.cpp @@ -63,13 +63,19 @@ namespace Opm /// Construct a 2d cartesian grid with cells of unit size. GridManager::GridManager(int nx, int ny) { - ug_ = create_grid_cart2d(nx, ny); + ug_ = create_grid_cart2d(nx, ny, 1.0, 1.0); if (!ug_) { THROW("Failed to construct grid."); } } - + GridManager::GridManager(int nx, int ny,double dx, double dy) + { + ug_ = create_grid_cart2d(nx, ny, dx, dy); + if (!ug_) { + THROW("Failed to construct grid."); + } + } /// Construct a 3d cartesian grid with cells of unit size. diff --git a/opm/core/GridManager.hpp b/opm/core/GridManager.hpp index 36846f2c..38fb09d7 100644 --- a/opm/core/GridManager.hpp +++ b/opm/core/GridManager.hpp @@ -47,6 +47,9 @@ namespace Opm /// Construct a 2d cartesian grid with cells of unit size. GridManager(int nx, int ny); + /// Construct a 2d cartesian grid with cells of size [dx, dy]. + GridManager(int nx, int ny, double dx, double dy); + /// Construct a 3d cartesian grid with cells of unit size. GridManager(int nx, int ny, int nz); diff --git a/opm/core/grid/cart_grid.c b/opm/core/grid/cart_grid.c index 989bbfe6..8adb2bd1 100644 --- a/opm/core/grid/cart_grid.c +++ b/opm/core/grid/cart_grid.c @@ -91,7 +91,7 @@ static void fill_cart_geometry_2d(struct UnstructuredGrid *G, const double *y); struct UnstructuredGrid* -create_grid_cart2d(int nx, int ny) +create_grid_cart2d(int nx, int ny, double dx, double dy) { int i; double *x, *y; @@ -104,8 +104,8 @@ create_grid_cart2d(int nx, int ny) G = NULL; } else { - for (i = 0; i < nx + 1; i++) { x[i] = i; } - for (i = 0; i < ny + 1; i++) { y[i] = i; } + for (i = 0; i < nx + 1; i++) { x[i] = i*dx; } + for (i = 0; i < ny + 1; i++) { y[i] = i*dy; } G = create_grid_tensor2d(nx, ny, x, y); } diff --git a/opm/core/grid/cart_grid.h b/opm/core/grid/cart_grid.h index 2ab9c673..527ef656 100644 --- a/opm/core/grid/cart_grid.h +++ b/opm/core/grid/cart_grid.h @@ -44,7 +44,7 @@ extern "C" { #endif -struct UnstructuredGrid; + /** * Form geometrically Cartesian grid in two space dimensions with unit-sized @@ -52,12 +52,14 @@ struct UnstructuredGrid; * * @param[in] nx Number of cells in @c x direction. * @param[in] ny Number of cells in @c y direction. + * @param[in] dx Length, in meters, of each cell's @c x extent. + * @param[in] dy Length, in meters, of each cell's @c y extent. * * @return Fully formed grid structure containing valid geometric primitives. * Must be destroyed using function destroy_grid(). */ struct UnstructuredGrid * -create_grid_cart2d(int nx, int ny); +create_grid_cart2d(int nx, int ny,double dx, double dy); /** diff --git a/opm/core/grid/cpgpreprocess/geometry.c b/opm/core/grid/cpgpreprocess/geometry.c index 769d5028..be833a46 100644 --- a/opm/core/grid/cpgpreprocess/geometry.c +++ b/opm/core/grid/cpgpreprocess/geometry.c @@ -2,6 +2,7 @@ * Copyright 2010 (c) SINTEF ICT, Applied Mathematics. * Jostein R. Natvig */ +#include #include #include #include "geometry.h" @@ -47,11 +48,20 @@ compute_face_geometry_3d(double *coords, int nfaces, double cface[3] = {0}; double n[3] = {0}; - const double twothirds = 0.666666666666666666666666666667; + double twothirds = 0.666666666666666666666666666667; + double a; + int num_face_nodes; + double area; + /*#pragma omp parallel for */ + + /*#pragma omp parallel for shared(fnormals,fcentroids,fareas)*/ +#pragma omp parallel for default(none) \ + private(f,x,u,v,w,i,k,node,cface,n,a,num_face_nodes,area) \ + shared(fnormals,fcentroids,fareas \ + ,coords, nfaces, nodepos, facenodes) \ + firstprivate(ndims, twothirds) for (f=0; f0)) + /* if(!(a>0)) { fprintf(stderr, "Internal error in compute_face_geometry."); } - + */ /* face normal */ for (i=0; i 0.0)) { - fprintf(stderr, - "Internal error in mex_compute_geometry(%*d): " - "negative volume\n", ndigits, c); - } - for (i=0; i* ptr = new SaturationPropsFromDeck(); @@ -41,9 +44,13 @@ namespace Opm BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, - const parameter::ParameterGroup& param) + const parameter::ParameterGroup& param, + bool init_rock) { - rock_.init(deck, grid); + if(init_rock){ + rock_.init(deck, grid); + } + const int pvt_samples = param.getDefault("pvt_tab_size", 200); pvt_.init(deck, pvt_samples); diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index aa5b590f..da8732a5 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -45,7 +45,7 @@ namespace Opm /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); + const UnstructuredGrid& grid, bool init_rock=true ); /// Initialize from deck, grid and parameters. /// \param[in] deck Deck input parser @@ -60,7 +60,8 @@ namespace Opm /// be done, and the input fluid data used directly for linear interpolation. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, - const parameter::ParameterGroup& param); + const parameter::ParameterGroup& param, + bool init_rock=true); /// Destructor. virtual ~BlackoilPropertiesFromDeck(); diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp index 0de37d55..132d29be 100644 --- a/opm/core/simulator/SimulatorCompressibleTwophase.cpp +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -46,7 +46,7 @@ #include #include #include -#include +#include #include #include @@ -101,7 +101,7 @@ namespace Opm const double* gravity_; // Solvers CompressibleTpfa psolver_; - TransportModelCompressibleTwophase tsolver_; + TransportSolverCompressibleTwophaseReorder tsolver_; // Needed by column-based gravity segregation solver. std::vector< std::vector > columns_; // Misc. data diff --git a/opm/core/simulator/SimulatorIncompTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp index 2302d47b..810df899 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -45,8 +45,9 @@ #include #include #include -#include - +//#include +//#include +#include #include #include #include @@ -98,7 +99,9 @@ namespace Opm const FlowBoundaryConditions* bcs_; // Solvers IncompTpfa psolver_; - TransportModelTwophase tsolver_; + // this should maybe be packed in a separate file + boost::scoped_ptr tsolver_; + //ImpliciteTwoPhaseTransportSolver tsolver_; // Needed by column-based gravity segregation solver. std::vector< std::vector > columns_; // Misc. data @@ -316,7 +319,7 @@ namespace Opm WellsManager& wells_manager, const std::vector& src, const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, + LinearSolverInterface& linsolverp, const double* gravity) : grid_(grid), props_(props), @@ -325,15 +328,49 @@ namespace Opm wells_(wells_manager.c_wells()), src_(src), bcs_(bcs), - psolver_(grid, props, rock_comp_props, linsolver, + psolver_(grid, props, rock_comp_props, linsolverp, param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_maxiter", 10), - gravity, wells_manager.c_wells(), src, bcs), - tsolver_(grid, props, - param.getDefault("nl_tolerance", 1e-9), - param.getDefault("nl_maxiter", 30)) + gravity, wells_manager.c_wells(), src, bcs) { + const bool use_reorder = param.getDefault("use_reorder", true); + if(use_reorder){ + /* + tsolver_.reset(new Opm::TransportModelTwoPhase(grid, props, rock_comp_props, linsolver, + param.getDefault("nl_pressure_residual_tolerance", 0.0), + param.getDefault("nl_pressure_change_tolerance", 1.0), + param.getDefault("nl_pressure_maxiter", 10), + gravity, wells_manager.c_wells(), src, bcs)); + */ + }else{ + //Opm::ImplicitTransportDetails::NRReport rpt; + Opm::ImplicitTransportDetails::NRControl ctrl; + ctrl.max_it = param.getDefault("max_it", 20); + ctrl.verbosity = param.getDefault("verbosity", 0); + ctrl.max_it_ls = param.getDefault("max_it_ls", 5); + const bool guess_old_solution = param.getDefault("guess_old_solution", false); + Opm::SimpleFluid2pWrappingProps fluid(props); + std::vector porevol; + //if (rock_comp->isActive()) { + // computePorevolume(grid, props->porosity(), rock_comp, state.pressure(), porevol); + //} else { + computePorevolume(grid, props.porosity(), porevol); + //} + SinglePointUpwindTwoPhase + model(fluid, grid, porevol, gravity, guess_old_solution); + model.initGravityTrans(grid_, psolver_.getHalfTrans()); + tsolver_.reset(new Opm::ImplicitTwoPhaseTransportSolver( + wells_manager, + *rock_comp_props, + ctrl, + model, + grid, + props, + param)); + + } + // For output. output_ = param.getDefault("output", true); if (output_) { @@ -356,12 +393,7 @@ namespace Opm // Transport related init. num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); - use_segregation_split_ = param.getDefault("use_segregation_split", false); - if (gravity != 0 && use_segregation_split_){ - tsolver_.initGravity(gravity); - extractColumn(grid_, columns_); - } - + use_segregation_split_ = param.getDefault("use_segregation_split", false); // Misc init. const int num_cells = grid.number_of_cells; allcells_.resize(num_cells); @@ -427,9 +459,6 @@ namespace Opm outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputVectorMatlab(std::string("reorder_it"), - tsolver_.getReorderIterations(), - timer.currentStepNum(), output_dir_); } SimulatorReport sreport; @@ -523,8 +552,14 @@ namespace Opm double injected[2] = { 0.0 }; double produced[2] = { 0.0 }; for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { - tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], - stepsize, state.saturation()); + //tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], + // stepsize, state.saturation()); + tsolver_->solve(&transport_src[0], + &porevol[0], + stepsize, + state, + well_state); + double substep_injected[2] = { 0.0 }; double substep_produced[2] = { 0.0 }; Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, @@ -533,17 +568,18 @@ namespace Opm injected[1] += substep_injected[1]; produced[0] += substep_produced[0]; produced[1] += substep_produced[1]; - if (use_segregation_split_) { - tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation()); - } + + watercut.push(timer.currentTime() + timer.currentStepLength(), produced[0]/(produced[0] + produced[1]), tot_produced[0]/tot_porevol_init); + if (wells_) { wellreport.push(props_, *wells_, state.saturation(), timer.currentTime() + timer.currentStepLength(), well_state.bhp(), well_state.perfRates()); } + } transport_timer.stop(); double tt = transport_timer.secsSinceStart(); @@ -571,9 +607,6 @@ namespace Opm outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputVectorMatlab(std::string("reorder_it"), - tsolver_.getReorderIterations(), - timer.currentStepNum(), output_dir_); outputWaterCut(watercut, output_dir_); if (wells_) { outputWellReport(wellreport, output_dir_); diff --git a/opm/core/simulator/SimulatorIncompTwophaseReorder.cpp b/opm/core/simulator/SimulatorIncompTwophaseReorder.cpp new file mode 100644 index 00000000..322ec37d --- /dev/null +++ b/opm/core/simulator/SimulatorIncompTwophaseReorder.cpp @@ -0,0 +1,594 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + + +namespace Opm +{ + + class SimulatorIncompTwophaseReorder::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity); + + SimulatorReport run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state); + + private: + // Data. + // Parameters for output. + bool output_; + bool output_vtk_; + std::string output_dir_; + int output_interval_; + // Parameters for well control + bool check_well_controls_; + int max_well_control_iterations_; + // Parameters for transport solver. + int num_transport_substeps_; + bool use_segregation_split_; + // Observed objects. + const UnstructuredGrid& grid_; + const IncompPropertiesInterface& props_; + const RockCompressibility* rock_comp_props_; + WellsManager& wells_manager_; + const Wells* wells_; + const std::vector& src_; + const FlowBoundaryConditions* bcs_; + // Solvers + IncompTpfa psolver_; + TransportSolverTwophaseReorder tsolver_; + // Needed by column-based gravity segregation solver. + std::vector< std::vector > columns_; + // Misc. data + std::vector allcells_; + }; + + + + + SimulatorIncompTwophaseReorder::SimulatorIncompTwophaseReorder(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity)); + } + + + + + + SimulatorReport SimulatorIncompTwophaseReorder::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { + return pimpl_->run(timer, state, well_state); + } + + static void reportVolumes(std::ostream &os, double satvol[2], double tot_porevol_init, + double tot_injected[2], double tot_produced[2], + double injected[2], double produced[2], + double init_satvol[2]) + { + std::cout.precision(5); + const int width = 18; + os << "\nVolume balance report (all numbers relative to total pore volume).\n"; + os << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; + os << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init << std::endl; + os << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init << std::endl; + os << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; + os << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; + os << " In-place + prod - inj: " + << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init + << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl; + os << " Init - now - pr + inj: " + << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init + << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init + << std::endl; + os.precision(8); + } + + static void outputStateVtk(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + THROW("Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + } + + static void outputVectorMatlab(const std::string& name, + const std::vector& vec, + const int step, + const std::string& output_dir) + { + std::ostringstream fname; + fname << output_dir << "/" << name; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + std::copy(vec.begin(), vec.end(), std::ostream_iterator(file, "\n")); + } + + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + watercut.write(os); + } + + + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + wellreport.write(os); + } + + + static bool allNeumannBCs(const FlowBoundaryConditions* bcs) + { + if (bcs == NULL) { + return true; + } else { + return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE) + == bcs->type + bcs->nbc; + } + } + + + static bool allRateWells(const Wells* wells) + { + if (wells == NULL) { + return true; + } + const int nw = wells->number_of_wells; + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells->ctrls[w]; + if (wc->current >= 0) { + if (wc->type[wc->current] == BHP) { + return false; + } + } + } + return true; + } + + + + + + SimulatorIncompTwophaseReorder::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + rock_comp_props_(rock_comp_props), + wells_manager_(wells_manager), + wells_(wells_manager.c_wells()), + src_(src), + bcs_(bcs), + psolver_(grid, props, rock_comp_props, linsolver, + param.getDefault("nl_pressure_residual_tolerance", 0.0), + param.getDefault("nl_pressure_change_tolerance", 1.0), + param.getDefault("nl_pressure_maxiter", 10), + gravity, wells_manager.c_wells(), src, bcs), + tsolver_(grid, props, + param.getDefault("nl_tolerance", 1e-9), + param.getDefault("nl_maxiter", 30)) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + + // Well control related init. + check_well_controls_ = param.getDefault("check_well_controls", false); + max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); + + // Transport related init. + num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); + use_segregation_split_ = param.getDefault("use_segregation_split", false); + if (gravity != 0 && use_segregation_split_){ + tsolver_.initGravity(gravity); + extractColumn(grid_, columns_); + } + + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + SimulatorReport SimulatorIncompTwophaseReorder::Impl::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { + std::vector transport_src; + + // Initialisation. + std::vector porevol; + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } else { + computePorevolume(grid_, props_.porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + std::vector initial_porevol = porevol; + + // Main simulation loop. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + double init_satvol[2] = { 0.0 }; + double satvol[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol); + std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init + << " " << init_satvol[1]/tot_porevol_init << std::endl; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); + Opm::WellReport wellreport; + std::vector fractional_flows; + std::vector well_resflows_phase; + if (wells_) { + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); + } + std::fstream tstep_os; + if (output_) { + std::string filename = output_dir_ + "/step_timing.param"; + tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); + } + for (; !timer.done(); ++timer) { + // Report timestep and (optionally) write state to disk. + step_timer.start(); + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputVectorMatlab(std::string("reorder_it"), + tsolver_.getReorderIterations(), + timer.currentStepNum(), output_dir_); + } + + SimulatorReport sreport; + + // Solve pressure equation. + if (check_well_controls_) { + computeFractionalFlow(props_, allcells_, state.saturation(), fractional_flows); + wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); + } + bool well_control_passed = !check_well_controls_; + int well_control_iteration = 0; + do { + // Run solver. + pressure_timer.start(); + std::vector initial_pressure = state.pressure(); + psolver_.solve(timer.currentStepLength(), state, well_state); + + // Renormalize pressure if rock is incompressible, and + // there are no pressure conditions (bcs or wells). + // It is deemed sufficient for now to renormalize + // using geometric volume instead of pore volume. + if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive()) + && allNeumannBCs(bcs_) && allRateWells(wells_)) { + // Compute average pressures of previous and last + // step, and total volume. + double av_prev_press = 0.0; + double av_press = 0.0; + double tot_vol = 0.0; + const int num_cells = grid_.number_of_cells; + for (int cell = 0; cell < num_cells; ++cell) { + av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell]; + av_press += state.pressure()[cell]*grid_.cell_volumes[cell]; + tot_vol += grid_.cell_volumes[cell]; + } + // Renormalization constant + const double ren_const = (av_prev_press - av_press)/tot_vol; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += ren_const; + } + const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells; + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += ren_const; + } + } + + // Stop timer and report. + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + sreport.pressure_time = pt; + + // Optionally, check if well controls are satisfied. + if (check_well_controls_) { + Opm::computePhaseFlowRatesPerWell(*wells_, + well_state.perfRates(), + fractional_flows, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { + THROW("Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); + } + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; + } + } + } while (!well_control_passed); + + // Update pore volumes if rock is compressible. + if (rock_comp_props_ && rock_comp_props_->isActive()) { + initial_porevol = porevol; + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } + + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, + wells_, well_state.perfRates(), transport_src); + + // Solve transport. + transport_timer.start(); + double stepsize = timer.currentStepLength(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], + stepsize, state.saturation()); + double substep_injected[2] = { 0.0 }; + double substep_produced[2] = { 0.0 }; + Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, + substep_injected, substep_produced); + injected[0] += substep_injected[0]; + injected[1] += substep_injected[1]; + produced[0] += substep_produced[0]; + produced[1] += substep_produced[1]; + if (use_segregation_split_) { + tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation()); + } + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } + } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + sreport.transport_time = tt; + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.saturation(), satvol); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + reportVolumes(std::cout,satvol, tot_porevol_init, + tot_injected, tot_produced, + injected, produced, + init_satvol); + sreport.total_time = step_timer.secsSinceStart(); + if (output_) { + sreport.reportParam(tstep_os); + } + } + + if (output_) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputVectorMatlab(std::string("reorder_it"), + tsolver_.getReorderIterations(), + timer.currentStepNum(), output_dir_); + outputWaterCut(watercut, output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } + tstep_os.close(); + } + + total_timer.stop(); + + SimulatorReport report; + report.pressure_time = ptime; + report.transport_time = ttime; + report.total_time = total_timer.secsSinceStart(); + return report; + } + + +} // namespace Opm diff --git a/opm/core/simulator/SimulatorIncompTwophaseReorder.hpp b/opm/core/simulator/SimulatorIncompTwophaseReorder.hpp new file mode 100644 index 00000000..a9616f02 --- /dev/null +++ b/opm/core/simulator/SimulatorIncompTwophaseReorder.hpp @@ -0,0 +1,99 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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 . +*/ + +#ifndef OPM_SIMULATORINCOMPTWOPHASEREORDER_HEADER_INCLUDED +#define OPM_SIMULATORINCOMPTWOPHASEREORDER_HEADER_INCLUDED + +#include +#include + +struct UnstructuredGrid; +struct Wells; +struct FlowBoundaryConditions; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class IncompPropertiesInterface; + class RockCompressibility; + class WellsManager; + class LinearSolverInterface; + class SimulatorTimer; + class TwophaseState; + class WellState; + struct SimulatorReport; + + /// Class collecting all necessary components for a two-phase simulation. + class SimulatorIncompTwophaseReorder + { + public: + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] well_manager well manager, may manage no (null) wells + /// \param[in] src source terms + /// \param[in] bcs boundary conditions, treat as all noflow if null + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorIncompTwophaseReorder(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity); + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + /// \return simulation report, with timing data + SimulatorReport run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state); + + private: + class Impl; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + boost::shared_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORINCOMPTWOPHASEREORDER_HEADER_INCLUDED diff --git a/opm/core/transport/ImplicitTwoPhaseTransportSolver.cpp b/opm/core/transport/ImplicitTwoPhaseTransportSolver.cpp new file mode 100644 index 00000000..3a8ee384 --- /dev/null +++ b/opm/core/transport/ImplicitTwoPhaseTransportSolver.cpp @@ -0,0 +1,86 @@ +/*=========================================================================== +// +// File: ImpliciteTwoPhaseTransportSolver.cpp +// +// Author: hnil +// +// Created: 9 Nov 2012 +//==========================================================================*/ +/* + 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 +namespace Opm{ + + ImplicitTwoPhaseTransportSolver::ImplicitTwoPhaseTransportSolver( + const Opm::WellsManager& wells, + const Opm::RockCompressibility& rock_comp, + const ImplicitTransportDetails::NRControl& ctrl, + SinglePointUpwindTwoPhase& model, + const UnstructuredGrid& grid, + const Opm::IncompPropertiesInterface& props, + const parameter::ParameterGroup& param) + : tsolver_(model), + grid_(grid), + ctrl_(ctrl), + props_(props), + rock_comp_(rock_comp), + wells_(wells) + { + tsrc_ = create_transport_source(2, 2); + //linsolver_(param), + //src_(num_cells, 0.0); + } + + void ImplicitTwoPhaseTransportSolver::solve(const double* porevolume, + const double* source, + const double dt, + TwophaseState& state, + WellState& well_state) + { + std::vector porevol; + if (rock_comp_.isActive()) { + computePorevolume(grid_, props_.porosity(), rock_comp_, state.pressure(), porevol); + } + std::vector src(grid_.number_of_cells, 0.0); + //Opm::wellsToSrc(*wells->c_wells(), num_cells, src); + Opm::computeTransportSource(grid_, src, state.faceflux(), 1.0, + wells_.c_wells(), well_state.perfRates(), src); + double ssrc[] = { 1.0, 0.0 }; + double ssink[] = { 0.0, 1.0 }; + double zdummy[] = { 0.0, 0.0 }; + for (int cell = 0; cell < grid_.number_of_cells; ++cell) { + clear_transport_source(tsrc_); + 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_); + } + } + // Boundary conditions. + Opm::ImplicitTransportDetails::NRReport rpt; + tsolver_.solve(grid_, tsrc_, dt, ctrl_, state, linsolver_, rpt); + std::cout << rpt; + + } +} diff --git a/opm/core/transport/ImplicitTwoPhaseTransportSolver.hpp b/opm/core/transport/ImplicitTwoPhaseTransportSolver.hpp new file mode 100644 index 00000000..2d3d7710 --- /dev/null +++ b/opm/core/transport/ImplicitTwoPhaseTransportSolver.hpp @@ -0,0 +1,129 @@ +/*=========================================================================== +// +// File: ImpliciteTwoPhaseTransportSolver.hpp +// +// Author: hnil +// +// Created: 9 Nov 2012 +//==========================================================================*/ +/* + 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 . +*/ + +#ifndef IMPLICITETWOPHASETRANSPORTSOLVER_HPP +#define IMPLICITTWOPHASETRANSPORTSOLVER_HPP +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +namespace Opm{ + // implicite transprot solver + class ImplicitTwoPhaseTransportSolver : public TwoPhaseTransportSolver + { + public: + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + /// \param[in] props Rock and fluid properties. + /// \param[in] tol Tolerance used in the solver. + /// \param[in] maxit Maximum number of non-linear iterations used. + ImplicitTwoPhaseTransportSolver( + const Opm::WellsManager& wells, + const Opm::RockCompressibility& rock_comp, + const ImplicitTransportDetails::NRControl& ctrl, + SinglePointUpwindTwoPhase& model, + const UnstructuredGrid& grid, + const Opm::IncompPropertiesInterface& props, + const parameter::ParameterGroup& param); + + ~ImplicitTwoPhaseTransportSolver(){ + destroy_transport_source(tsrc_); + } + + /// Solve for saturation at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[in] dt Time step. + /// \param[in, out] saturation Phase saturations. + void solve(const double* porevolume, + const double* source, + const double dt, + Opm::TwophaseState& state, + Opm::WellState& well_state); + + private: + ImplicitTwoPhaseTransportSolver(const ImplicitTwoPhaseTransportSolver&); + ImplicitTwoPhaseTransportSolver& operator=(const ImplicitTwoPhaseTransportSolver&); + typedef Opm::SimpleFluid2pWrappingProps TwophaseFluid; + typedef Opm::SinglePointUpwindTwoPhase TransportModel; + + //using namespace ImplicitTransportDefault; + + typedef ImplicitTransportDefault::NewtonVectorCollection< ::std::vector > NVecColl; + typedef ImplicitTransportDefault::JacobianSystem < struct CSRMatrix, NVecColl > JacSys; + + template + class MaxNorm { + public: + static double + norm(const Vector& v) { + return ImplicitTransportDefault::AccumulationNorm ::norm(v); + } + }; + + typedef Opm::ImplicitTransport TransportSolver; + //Opm::LinearSolverFactory + Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver linsolver_; + TransportSolver tsolver_; + const UnstructuredGrid& grid_; + const Opm::ImplicitTransportDetails::NRControl& ctrl_; + const Opm::IncompPropertiesInterface& props_; + const Opm::RockCompressibility& rock_comp_; + const Opm::WellsManager& wells_; + + TransportSource* tsrc_; + + }; +} +#endif // IMPLICITETWOPHASETRANSPORTSOLVER_HPP diff --git a/opm/core/transport/SimpleFluid2pWrappingProps.hpp b/opm/core/transport/SimpleFluid2pWrappingProps.hpp new file mode 100644 index 00000000..0ace1294 --- /dev/null +++ b/opm/core/transport/SimpleFluid2pWrappingProps.hpp @@ -0,0 +1,70 @@ +/*=========================================================================== +// +// File: SimpleFluid2pWrappingProps.hpp +// +// Author: hnil +// +// Created: 15 Nov 2012 +//==========================================================================*/ +/* + 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 . +*/ + + +#ifndef SIMPLEFLUID2PWRAPPINGPROPS_HPP +#define SIMPLEFLUID2PWRAPPINGPROPS_HPP + +#include +#include + +namespace Opm{ + class SimpleFluid2pWrappingProps + { + public: + SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props); + + double density(int phase) const; + + + template + void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const; + + + template + void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const; + + double s_min(int c) const; + + + double s_max(int c) const; + + + private: + const Opm::IncompPropertiesInterface& props_; + std::vector smin_; + std::vector smax_; + }; +} + +#include +#endif // SIMPLEFLUID2PWRAPPINGPROPS_HPP diff --git a/opm/core/transport/SimpleFluid2pWrappingProps_impl.hpp b/opm/core/transport/SimpleFluid2pWrappingProps_impl.hpp new file mode 100644 index 00000000..c1aba4ad --- /dev/null +++ b/opm/core/transport/SimpleFluid2pWrappingProps_impl.hpp @@ -0,0 +1,104 @@ +/*=========================================================================== +// +// File: SimpleFluid2pWrappingProps.cpp +// +// Author: hnil +// +// Created: 15 Nov 2012 +//==========================================================================*/ +/* + 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 . +*/ +#ifndef SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP +#define SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP + +#include +#include +#include +namespace Opm{ + + inline SimpleFluid2pWrappingProps::SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props) + : props_(props), + smin_(props.numCells()*props.numPhases()), + smax_(props.numCells()*props.numPhases()) + { + if (props.numPhases() != 2) { + THROW("SimpleFluid2pWrapper requires 2 phases."); + } + 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]); + } + + inline double SimpleFluid2pWrappingProps::density(int phase) const + { + return props_.density()[phase]; + } + + template + void SimpleFluid2pWrappingProps::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]; + } + + template + void SimpleFluid2pWrappingProps::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]; + ASSERT(pcow[1] == 0.0); + dpcap = dpcow[0]; + ASSERT(dpcow[1] == 0.0); + ASSERT(dpcow[2] == 0.0); + ASSERT(dpcow[3] == 0.0); + } + + inline double SimpleFluid2pWrappingProps::s_min(int c) const + { + return smin_[2*c + 0]; + } + + inline double SimpleFluid2pWrappingProps::s_max(int c) const + { + return smax_[2*c + 0]; + } + +} +#endif // SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP diff --git a/opm/core/transport/TwoPhaseTransportSolver.cpp b/opm/core/transport/TwoPhaseTransportSolver.cpp new file mode 100644 index 00000000..6515670e --- /dev/null +++ b/opm/core/transport/TwoPhaseTransportSolver.cpp @@ -0,0 +1,32 @@ +/*=========================================================================== +// +// File: TwoPhaseTransportSolver.cpp +// +// Author: hnil +// +// Created: 9 Nov 2012 +//==========================================================================*/ +/* + 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 "TwoPhaseTransportSolver.hpp" + + diff --git a/opm/core/transport/TwoPhaseTransportSolver.hpp b/opm/core/transport/TwoPhaseTransportSolver.hpp new file mode 100644 index 00000000..7e678f06 --- /dev/null +++ b/opm/core/transport/TwoPhaseTransportSolver.hpp @@ -0,0 +1,60 @@ +/*=========================================================================== +// +// File: TwoPhaseTransportSolver.hpp +// +// Author: hnil +// +// Created: 9 Nov 2012 +//==========================================================================*/ +/* + 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 . +*/ + + +#ifndef TWOPHASETRANSPORTSOLVER_HPP +#define TWOPHASETRANSPORTSOLVER_HPP + +#include +#include +namespace Opm +{ + + /// Base class for tranport solvers + class TwoPhaseTransportSolver + { + public: + /// Virtual destructor to enable inheritance. + virtual ~TwoPhaseTransportSolver() {} + + /// Solve for saturation at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[in] dt Time step. + /// \param[in, out] saturation Phase saturations. + virtual void solve(const double* porevolume, + const double* source, + const double dt, + TwophaseState& state, + WellState& wstate) = 0; + }; + +} + +#endif // TWOPHASETRANSPORTSOLVER_HPP diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/ReorderSolverInterface.cpp similarity index 92% rename from opm/core/transport/reorder/TransportModelInterface.cpp rename to opm/core/transport/reorder/ReorderSolverInterface.cpp index 0a542b88..79969aaa 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/opm/core/transport/reorder/ReorderSolverInterface.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include @@ -26,7 +26,7 @@ #include -void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux) +void Opm::ReorderSolverInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux) { // Compute reordered sequence of single-cell problems sequence_.resize(grid.number_of_cells); diff --git a/opm/core/transport/reorder/TransportModelInterface.hpp b/opm/core/transport/reorder/ReorderSolverInterface.hpp similarity index 91% rename from opm/core/transport/reorder/TransportModelInterface.hpp rename to opm/core/transport/reorder/ReorderSolverInterface.hpp index 763ffd2e..a1071ee5 100644 --- a/opm/core/transport/reorder/TransportModelInterface.hpp +++ b/opm/core/transport/reorder/ReorderSolverInterface.hpp @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED -#define OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED +#ifndef OPM_REORDERSOLVERINTERFACE_HEADER_INCLUDED +#define OPM_REORDERSOLVERINTERFACE_HEADER_INCLUDED #include @@ -35,10 +35,10 @@ namespace Opm /// class.) The reorderAndTransport() method is provided as an aid /// to implementing solve() in subclasses, together with the /// sequence() and components() methods for accessing the ordering. - class TransportModelInterface + class ReorderSolverInterface { public: - virtual ~TransportModelInterface() {} + virtual ~ReorderSolverInterface() {} private: virtual void solveSingleCell(const int cell) = 0; virtual void solveMultiCell(const int num_cells, const int* cells) = 0; diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TofDiscGalReorder.cpp similarity index 97% rename from opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp rename to opm/core/transport/reorder/TofDiscGalReorder.cpp index 54a6e585..2ef010a4 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TofDiscGalReorder.cpp @@ -19,7 +19,7 @@ #include #include -#include +#include #include #include #include @@ -34,10 +34,6 @@ namespace Opm { - // --------------- Methods of TransportModelTracerTofDiscGal --------------- - - - /// Construct solver. /// \param[in] grid A 2d or 3d grid. /// \param[in] param Parameters for the solver. @@ -55,8 +51,8 @@ namespace Opm /// computing (unlimited) solution. /// AsSimultaneousPostProcess Apply to each cell independently, using un- /// limited solution in neighbouring cells. - TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, - const parameter::ParameterGroup& param) + TofDiscGalReorder::TofDiscGalReorder(const UnstructuredGrid& grid, + const parameter::ParameterGroup& param) : grid_(grid), use_cvi_(false), use_limiter_(false), @@ -127,7 +123,7 @@ namespace Opm /// cell comes before the K coefficients corresponding /// to the second cell etc. /// K depends on degree and grid dimension. - void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux, + void TofDiscGalReorder::solveTof(const double* darcyflux, const double* porevolume, const double* source, std::vector& tof_coeff) @@ -173,7 +169,7 @@ namespace Opm - void TransportModelTracerTofDiscGal::solveSingleCell(const int cell) + void TofDiscGalReorder::solveSingleCell(const int cell) { // Residual: // For each cell K, basis function b_j (spanning V_h), @@ -396,7 +392,7 @@ namespace Opm - void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells) + void TofDiscGalReorder::solveMultiCell(const int num_cells, const int* cells) { std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; for (int i = 0; i < num_cells; ++i) { diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TofDiscGalReorder.hpp similarity index 92% rename from opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp rename to opm/core/transport/reorder/TofDiscGalReorder.hpp index cf55ca8c..ba6c72dc 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TofDiscGalReorder.hpp @@ -17,10 +17,10 @@ along with OPM. If not, see . */ -#ifndef OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED -#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED +#ifndef OPM_TOFDISCGALREORDER_HEADER_INCLUDED +#define OPM_TOFDISCGALREORDER_HEADER_INCLUDED -#include +#include #include #include #include @@ -45,7 +45,7 @@ namespace Opm /// \tau is specified to be zero on all inflow boundaries. /// The user may specify the polynomial degree of the basis function space /// used, but only degrees 0 and 1 are supported so far. - class TransportModelTracerTofDiscGal : public TransportModelInterface + class TofDiscGalReorder : public ReorderSolverInterface { public: /// Construct solver. @@ -68,8 +68,8 @@ namespace Opm /// computing (unlimited) solution. /// AsSimultaneousPostProcess Apply to each cell independently, using un- /// limited solution in neighbouring cells. - TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, - const parameter::ParameterGroup& param); + TofDiscGalReorder(const UnstructuredGrid& grid, + const parameter::ParameterGroup& param); /// Solve for time-of-flight. @@ -95,8 +95,8 @@ namespace Opm private: // Disable copying and assignment. - TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&); - TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&); + TofDiscGalReorder(const TofDiscGalReorder&); + TofDiscGalReorder& operator=(const TofDiscGalReorder&); // Data members const UnstructuredGrid& grid_; diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TofReorder.cpp similarity index 96% rename from opm/core/transport/reorder/TransportModelTracerTof.cpp rename to opm/core/transport/reorder/TofReorder.cpp index 27562844..ca5183e0 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TofReorder.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include @@ -30,8 +30,8 @@ namespace Opm /// Construct solver. /// \param[in] grid A 2d or 3d grid. - TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid, - const bool use_multidim_upwind) + TofReorder::TofReorder(const UnstructuredGrid& grid, + const bool use_multidim_upwind) : grid_(grid), darcyflux_(0), porevolume_(0), @@ -53,7 +53,7 @@ namespace Opm /// (+) inflow flux, /// (-) outflow flux. /// \param[out] tof Array of time-of-flight values. - void TransportModelTracerTof::solveTof(const double* darcyflux, + void TofReorder::solveTof(const double* darcyflux, const double* porevolume, const double* source, std::vector& tof) @@ -137,7 +137,7 @@ namespace Opm - void TransportModelTracerTof::solveSingleCell(const int cell) + void TofReorder::solveSingleCell(const int cell) { if (use_multidim_upwind_) { solveSingleCellMultidimUpwind(cell); @@ -243,7 +243,7 @@ namespace Opm - void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells) + void TofReorder::solveMultiCell(const int num_cells, const int* cells) { std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; for (int i = 0; i < num_cells; ++i) { diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TofReorder.hpp similarity index 91% rename from opm/core/transport/reorder/TransportModelTracerTof.hpp rename to opm/core/transport/reorder/TofReorder.hpp index 04983981..b1a9bf9c 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TofReorder.hpp @@ -17,10 +17,10 @@ along with OPM. If not, see . */ -#ifndef OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED -#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED +#ifndef OPM_TOFREORDER_HEADER_INCLUDED +#define OPM_TOFREORDER_HEADER_INCLUDED -#include +#include #include #include #include @@ -38,14 +38,14 @@ namespace Opm /// where v is the fluid velocity, \tau is time-of-flight and /// \phi is the porosity. This is a boundary value problem, where /// \tau is specified to be zero on all inflow boundaries. - class TransportModelTracerTof : public TransportModelInterface + class TofReorder : public ReorderSolverInterface { public: /// Construct solver. /// \param[in] grid A 2d or 3d grid. /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. - TransportModelTracerTof(const UnstructuredGrid& grid, - const bool use_multidim_upwind = false); + TofReorder(const UnstructuredGrid& grid, + const bool use_multidim_upwind = false); /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp similarity index 93% rename from opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp rename to opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp index 045a0203..04271b3b 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp @@ -18,7 +18,7 @@ */ -#include +#include #include #include #include @@ -39,7 +39,7 @@ namespace Opm typedef RegulaFalsi RootFinder; - TransportModelCompressibleTwophase::TransportModelCompressibleTwophase( + TransportSolverCompressibleTwophaseReorder::TransportSolverCompressibleTwophaseReorder( const UnstructuredGrid& grid, const Opm::BlackoilPropertiesInterface& props, const double tol, @@ -76,7 +76,7 @@ namespace Opm props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]); } - void TransportModelCompressibleTwophase::solve(const double* darcyflux, + void TransportSolverCompressibleTwophaseReorder::solve(const double* darcyflux, const double* pressure, const double* porevolume0, const double* porevolume, @@ -132,7 +132,7 @@ namespace Opm // We need the formula influx = B_i sum_{j->i} b_j v_{ij} - B_i q_w. // outflux = B_i sum_{i->j} b_i v_{ij} - B_i q = sum_{i->j} v_{ij} - B_i q // Influxes are negative, outfluxes positive. - struct TransportModelCompressibleTwophase::Residual + struct TransportSolverCompressibleTwophaseReorder::Residual { int cell; double B_cell; @@ -142,8 +142,8 @@ namespace Opm // @@@ TODO: figure out change to rock-comp. terms with fluid compr. double comp_term; // Now: used to be: q - sum_j v_ij double dtpv; // dt/pv(i) - const TransportModelCompressibleTwophase& tm; - explicit Residual(const TransportModelCompressibleTwophase& tmodel, int cell_index) + const TransportSolverCompressibleTwophaseReorder& tm; + explicit Residual(const TransportSolverCompressibleTwophaseReorder& tmodel, int cell_index) : tm(tmodel) { cell = cell_index; @@ -187,7 +187,7 @@ namespace Opm }; - void TransportModelCompressibleTwophase::solveSingleCell(const int cell) + void TransportSolverCompressibleTwophaseReorder::solveSingleCell(const int cell) { Residual res(*this, cell); int iters_used; @@ -196,7 +196,7 @@ namespace Opm } - void TransportModelCompressibleTwophase::solveMultiCell(const int num_cells, const int* cells) + void TransportSolverCompressibleTwophaseReorder::solveMultiCell(const int num_cells, const int* cells) { // Experiment: when a cell changes more than the tolerance, // mark all downwind cells as needing updates. After @@ -302,7 +302,7 @@ namespace Opm } - double TransportModelCompressibleTwophase::fracFlow(double s, int cell) const + double TransportSolverCompressibleTwophaseReorder::fracFlow(double s, int cell) const { double sat[2] = { s, 1.0 - s }; double mob[2]; @@ -321,15 +321,15 @@ namespace Opm // [[ incompressible was: r(s) = s - s0 + dt/pv*sum_{j adj i}( gravmod_ij * gf_ij ) ]] // // r(s) = s - B*z0 + dt/pv*( influx + outflux*f(s) ) - struct TransportModelCompressibleTwophase::GravityResidual + struct TransportSolverCompressibleTwophaseReorder::GravityResidual { int cell; int nbcell[2]; double s0; double dtpv; // dt/pv(i) double gf[2]; - const TransportModelCompressibleTwophase& tm; - explicit GravityResidual(const TransportModelCompressibleTwophase& tmodel, + const TransportSolverCompressibleTwophaseReorder& tm; + explicit GravityResidual(const TransportSolverCompressibleTwophaseReorder& tmodel, const std::vector& cells, const int pos, const double* gravflux) // Always oriented towards next in column. Size = colsize - 1. @@ -376,7 +376,7 @@ namespace Opm } }; - void TransportModelCompressibleTwophase::mobility(double s, int cell, double* mob) const + void TransportSolverCompressibleTwophaseReorder::mobility(double s, int cell, double* mob) const { double sat[2] = { s, 1.0 - s }; props_.relperm(1, sat, &cell, mob, 0); @@ -386,7 +386,7 @@ namespace Opm - void TransportModelCompressibleTwophase::initGravity(const double* grav) + void TransportSolverCompressibleTwophaseReorder::initGravity(const double* grav) { // Set up transmissibilities. std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); @@ -403,7 +403,7 @@ namespace Opm - void TransportModelCompressibleTwophase::initGravityDynamic() + void TransportSolverCompressibleTwophaseReorder::initGravityDynamic() { // Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f) // + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ] @@ -436,7 +436,7 @@ namespace Opm - void TransportModelCompressibleTwophase::solveSingleCellGravity(const std::vector& cells, + void TransportSolverCompressibleTwophaseReorder::solveSingleCellGravity(const std::vector& cells, const int pos, const double* gravflux) { @@ -451,7 +451,7 @@ namespace Opm - int TransportModelCompressibleTwophase::solveGravityColumn(const std::vector& cells) + int TransportSolverCompressibleTwophaseReorder::solveGravityColumn(const std::vector& cells) { // Set up column gravflux. const int nc = cells.size(); @@ -504,7 +504,7 @@ namespace Opm - void TransportModelCompressibleTwophase::solveGravity(const std::vector >& columns, + void TransportSolverCompressibleTwophaseReorder::solveGravity(const std::vector >& columns, const double dt, std::vector& saturation, std::vector& surfacevol) diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp similarity index 93% rename from opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp rename to opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp index 23542de8..c1c9f3fb 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp @@ -17,10 +17,10 @@ along with OPM. If not, see . */ -#ifndef OPM_TRANSPORTMODELCOMPRESSIBLETWOPHASE_HEADER_INCLUDED -#define OPM_TRANSPORTMODELCOMPRESSIBLETWOPHASE_HEADER_INCLUDED +#ifndef OPM_TRANSPORTSOLVERCOMPRESSIBLETWOPHASEREORDER_HEADER_INCLUDED +#define OPM_TRANSPORTSOLVERCOMPRESSIBLETWOPHASEREORDER_HEADER_INCLUDED -#include +#include #include struct UnstructuredGrid; @@ -32,7 +32,7 @@ namespace Opm /// Implements a reordering transport solver for compressible, /// non-miscible two-phase flow. - class TransportModelCompressibleTwophase : public TransportModelInterface + class TransportSolverCompressibleTwophaseReorder : public ReorderSolverInterface { public: /// Construct solver. @@ -40,7 +40,7 @@ namespace Opm /// \param[in] props Rock and fluid properties. /// \param[in] tol Tolerance used in the solver. /// \param[in] maxit Maximum number of non-linear iterations used. - TransportModelCompressibleTwophase(const UnstructuredGrid& grid, + TransportSolverCompressibleTwophaseReorder(const UnstructuredGrid& grid, const Opm::BlackoilPropertiesInterface& props, const double tol, const int maxit); diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp similarity index 94% rename from opm/core/transport/reorder/TransportModelTwophase.cpp rename to opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp index 1643dc27..db88d432 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include @@ -40,7 +40,7 @@ namespace Opm typedef RegulaFalsi RootFinder; - TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid& grid, + TransportSolverTwophaseReorder::TransportSolverTwophaseReorder(const UnstructuredGrid& grid, const Opm::IncompPropertiesInterface& props, const double tol, const int maxit) @@ -76,7 +76,7 @@ namespace Opm props.satRange(props.numCells(), &cells[0], &smin_[0], &smax_[0]); } - void TransportModelTwophase::solve(const double* darcyflux, + void TransportSolverTwophaseReorder::solve(const double* darcyflux, const double* porevolume, const double* source, const double dt, @@ -108,7 +108,7 @@ namespace Opm } - const std::vector& TransportModelTwophase::getReorderIterations() const + const std::vector& TransportSolverTwophaseReorder::getReorderIterations() const { return reorder_iterations_; } @@ -120,7 +120,7 @@ namespace Opm // // where influx is water influx, outflux is total outflux. // Influxes are negative, outfluxes positive. - struct TransportModelTwophase::Residual + struct TransportSolverTwophaseReorder::Residual { int cell; double s0; @@ -128,8 +128,8 @@ namespace Opm double outflux; // sum_j max(v_ij, 0) - q double comp_term; // q - sum_j v_ij double dtpv; // dt/pv(i) - const TransportModelTwophase& tm; - explicit Residual(const TransportModelTwophase& tmodel, int cell_index) + const TransportSolverTwophaseReorder& tm; + explicit Residual(const TransportSolverTwophaseReorder& tmodel, int cell_index) : tm(tmodel) { cell = cell_index; @@ -171,7 +171,7 @@ namespace Opm }; - void TransportModelTwophase::solveSingleCell(const int cell) + void TransportSolverTwophaseReorder::solveSingleCell(const int cell) { Residual res(*this, cell); // const double r0 = res(saturation_[cell]); @@ -225,7 +225,7 @@ namespace Opm // } // anon namespace - void TransportModelTwophase::solveMultiCell(const int num_cells, const int* cells) + void TransportSolverTwophaseReorder::solveMultiCell(const int num_cells, const int* cells) { // std::ofstream os("dump"); // std::copy(cells, cells + num_cells, std::ostream_iterator(os, "\n")); @@ -440,7 +440,7 @@ namespace Opm #endif // EXPERIMENT_GAUSS_SEIDEL } - double TransportModelTwophase::fracFlow(double s, int cell) const + double TransportSolverTwophaseReorder::fracFlow(double s, int cell) const { double sat[2] = { s, 1.0 - s }; double mob[2]; @@ -458,15 +458,15 @@ namespace Opm // // r(s) = s - s0 + dt/pv*sum_{j adj i}( gravmod_ij * gf_ij ). // - struct TransportModelTwophase::GravityResidual + struct TransportSolverTwophaseReorder::GravityResidual { int cell; int nbcell[2]; double s0; double dtpv; // dt/pv(i) double gf[2]; - const TransportModelTwophase& tm; - explicit GravityResidual(const TransportModelTwophase& tmodel, + const TransportSolverTwophaseReorder& tm; + explicit GravityResidual(const TransportSolverTwophaseReorder& tmodel, const std::vector& cells, const int pos, const double* gravflux) // Always oriented towards next in column. Size = colsize - 1. @@ -513,7 +513,7 @@ namespace Opm } }; - void TransportModelTwophase::mobility(double s, int cell, double* mob) const + void TransportSolverTwophaseReorder::mobility(double s, int cell, double* mob) const { double sat[2] = { s, 1.0 - s }; props_.relperm(1, sat, &cell, mob, 0); @@ -523,7 +523,7 @@ namespace Opm - void TransportModelTwophase::initGravity(const double* grav) + void TransportSolverTwophaseReorder::initGravity(const double* grav) { // Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j) std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); @@ -547,7 +547,7 @@ namespace Opm - void TransportModelTwophase::solveSingleCellGravity(const std::vector& cells, + void TransportSolverTwophaseReorder::solveSingleCellGravity(const std::vector& cells, const int pos, const double* gravflux) { @@ -564,7 +564,7 @@ namespace Opm - int TransportModelTwophase::solveGravityColumn(const std::vector& cells) + int TransportSolverTwophaseReorder::solveGravityColumn(const std::vector& cells) { // Set up column gravflux. const int nc = cells.size(); @@ -617,7 +617,7 @@ namespace Opm - void TransportModelTwophase::solveGravity(const std::vector >& columns, + void TransportSolverTwophaseReorder::solveGravity(const std::vector >& columns, const double* porevolume, const double dt, std::vector& saturation) diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp similarity index 93% rename from opm/core/transport/reorder/TransportModelTwophase.hpp rename to opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp index f03b6b7a..30a9825d 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp @@ -17,10 +17,10 @@ along with OPM. If not, see . */ -#ifndef OPM_TRANSPORTMODELTWOPHASE_HEADER_INCLUDED -#define OPM_TRANSPORTMODELTWOPHASE_HEADER_INCLUDED +#ifndef OPM_TRANSPORTSOLVERTWOPHASEREORDER_HEADER_INCLUDED +#define OPM_TRANSPORTSOLVERTWOPHASEREORDER_HEADER_INCLUDED -#include +#include #include #include #include @@ -32,7 +32,7 @@ namespace Opm class IncompPropertiesInterface; /// Implements a reordering transport solver for incompressible two-phase flow. - class TransportModelTwophase : public TransportModelInterface + class TransportSolverTwophaseReorder : public ReorderSolverInterface { public: /// Construct solver. @@ -40,7 +40,7 @@ namespace Opm /// \param[in] props Rock and fluid properties. /// \param[in] tol Tolerance used in the solver. /// \param[in] maxit Maximum number of non-linear iterations used. - TransportModelTwophase(const UnstructuredGrid& grid, + TransportSolverTwophaseReorder(const UnstructuredGrid& grid, const Opm::IncompPropertiesInterface& props, const double tol, const int maxit); diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 1bc7fa9c..5b1375b7 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -223,8 +223,7 @@ namespace Opm : w_(0) { } - - + /// Construct from existing wells object. WellsManager::WellsManager(struct Wells* W)