diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 2d58ec918..404452e0e 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -20,10 +20,11 @@ #include #include #include +#include #include #include -#include #include +#include namespace Opm { @@ -403,32 +404,6 @@ namespace Opm - // Initial version: only a constant interpolation. - static void interpolateVelocity(const UnstructuredGrid& grid, - const int cell, - const double* darcyflux, - const double* /*x*/, - double* v) - { - const int dim = grid.dimensions; - std::fill(v, v + dim, 0.0); - const double* cc = grid.cell_centroids + cell*dim; - for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) { - const int face = grid.cell_faces[hface]; - const double* fc = grid.face_centroids + face*dim; - double flux = 0.0; - if (cell == grid.face_cells[2*face]) { - flux = darcyflux[face]; - } else { - flux = -darcyflux[face]; - } - for (int dd = 0; dd < dim; ++dd) { - v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell]; - } - } - } - - // --------------- Methods of TransportModelTracerTofDiscGal --------------- @@ -436,11 +411,19 @@ namespace Opm /// Construct solver. /// \param[in] grid A 2d or 3d grid. - TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid) + /// \param[in] use_cvi If true, use corner point velocity interpolation. + /// Otherwise, use the basic constant interpolation. + TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, + const bool use_cvi) : grid_(grid), coord_(grid.dimensions), velocity_(grid.dimensions) { + if (use_cvi) { + velocity_interpolation_.reset(new VelocityInterpolationECVI(grid)); + } else { + velocity_interpolation_.reset(new VelocityInterpolationConstant(grid)); + } } @@ -482,9 +465,11 @@ namespace Opm tof_coeff_ = &tof_coeff[0]; rhs_.resize(num_basis); jac_.resize(num_basis*num_basis); + orig_jac_.resize(num_basis*num_basis); basis_.resize(num_basis); basis_nb_.resize(num_basis); grad_basis_.resize(num_basis*grid_.dimensions); + velocity_interpolation_->setupFluxes(darcyflux); reorderAndTransport(grid_, darcyflux); } @@ -560,7 +545,7 @@ namespace Opm quad.quadPtCoord(quad_pt, &coord_[0]); DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]); - interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]); + velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { @@ -633,6 +618,7 @@ namespace Opm std::vector piv(num_basis); MAT_SIZE_T ldb = num_basis; MAT_SIZE_T info = 0; + orig_jac_ = jac_; dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); if (info != 0) { // Print the local matrix and rhs. @@ -640,7 +626,7 @@ namespace Opm << " with A = \n"; for (int row = 0; row < n; ++row) { for (int col = 0; col < n; ++col) { - std::cerr << " " << jac_[row + n*col]; + std::cerr << " " << orig_jac_[row + n*col]; } std::cerr << '\n'; } diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 0285e1f02..9a2e3b9c2 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -21,15 +21,18 @@ #define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED #include +#include #include #include #include + struct UnstructuredGrid; namespace Opm { class IncompPropertiesInterface; + class VelocityInterpolationInterface; /// Implements a discontinuous Galerkin solver for /// (single-phase) time-of-flight using reordering. @@ -45,7 +48,10 @@ namespace Opm public: /// Construct solver. /// \param[in] grid A 2d or 3d grid. - TransportModelTracerTofDiscGal(const UnstructuredGrid& grid); + /// \param[in] use_cvi If true, use corner point velocity interpolation. + /// Otherwise, use the basic constant interpolation. + TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, + const bool use_cvi); /// Solve for time-of-flight. @@ -72,7 +78,12 @@ namespace Opm virtual void solveMultiCell(const int num_cells, const int* cells); private: + // Disable copying and assignment. + TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&); + TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&); + const UnstructuredGrid& grid_; + boost::shared_ptr velocity_interpolation_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell @@ -80,6 +91,7 @@ namespace Opm double* tof_coeff_; std::vector rhs_; // single-cell right-hand-side std::vector jac_; // single-cell jacobian + std::vector orig_jac_; // single-cell jacobian (copy) // Below: storage for quantities needed by solveSingleCell(). std::vector coord_; std::vector basis_; diff --git a/opm/core/utility/writeVtkData.hpp b/opm/core/utility/writeVtkData.hpp index 200b2b261..96cdfac22 100644 --- a/opm/core/utility/writeVtkData.hpp +++ b/opm/core/utility/writeVtkData.hpp @@ -26,25 +26,23 @@ #include #include #include +#include struct UnstructuredGrid; namespace Opm { - /// Intended to map strings (giving the output field names) to data. - typedef std::map*> DataMap; - /// Vtk output for cartesian grids. void writeVtkData(const std::tr1::array& dims, - const std::tr1::array& cell_size, - const DataMap& data, - std::ostream& os); + const std::tr1::array& cell_size, + const DataMap& data, + std::ostream& os); /// Vtk output for general grids. void writeVtkData(const UnstructuredGrid& grid, - const DataMap& data, - std::ostream& os); + const DataMap& data, + std::ostream& os); } // namespace Opm #endif // OPM_WRITEVTKDATA_HEADER_INCLUDED diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index d407c0335..8ffb30273 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -761,16 +761,17 @@ namespace Opm wells_->ctrls[self_index_]->current = ~ wells_->ctrls[self_index_]->current; } } - + std::pair WellNode::getWorstOffending(const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase, ProductionSpecification::ControlMode mode) { const int np = phaseUsage().num_phases; const int index = self_index_*np; - return std::make_pair(this, rateByMode(&well_reservoirrates_phase[index], - &well_surfacerates_phase[index], - mode)); + return std::pair(this, + rateByMode(&well_reservoirrates_phase[index], + &well_surfacerates_phase[index], + mode)); } void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode, @@ -782,7 +783,7 @@ namespace Opm && (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) { return; } - if (!wells_->type[self_index_] == INJECTOR) { + if (wells_->type[self_index_] != INJECTOR) { ASSERT(target == 0.0); return; } @@ -859,7 +860,7 @@ namespace Opm std::cout << "Returning" << std::endl; return; } - if (!wells_->type[self_index_] == PRODUCER) { + if (wells_->type[self_index_] != PRODUCER) { ASSERT(target == 0.0); return; }