diff --git a/opm/core/grid/FaceQuadrature.hpp b/opm/core/grid/FaceQuadrature.hpp index d36dc65c..505fc705 100644 --- a/opm/core/grid/FaceQuadrature.hpp +++ b/opm/core/grid/FaceQuadrature.hpp @@ -20,6 +20,7 @@ #ifndef OPM_FACEQUADRATURE_HEADER_INCLUDED #define OPM_FACEQUADRATURE_HEADER_INCLUDED +#include #include #include diff --git a/opm/core/props/IncompPropertiesShadow_impl.hpp b/opm/core/props/IncompPropertiesShadow_impl.hpp index 279e71a7..fef3a68a 100644 --- a/opm/core/props/IncompPropertiesShadow_impl.hpp +++ b/opm/core/props/IncompPropertiesShadow_impl.hpp @@ -4,6 +4,8 @@ #error Do not include IncompPropertiesShadow_impl.hpp directly! #endif /* OPM_INCOMPPROPERTIESSHADOW_HEADER_INCLUDED */ +#include + namespace Opm { /** diff --git a/opm/core/simulator/SimulatorIncompTwophase.hpp b/opm/core/simulator/SimulatorIncompTwophase.hpp index 437afe5f..d44cbbfb 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.hpp +++ b/opm/core/simulator/SimulatorIncompTwophase.hpp @@ -39,7 +39,7 @@ namespace Opm class TwophaseState; class WellState; struct SimulatorReport; - class Event; + struct Event; /// Class collecting all necessary components for a two-phase simulation. class SimulatorIncompTwophase diff --git a/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp b/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp index 4e5b4fce..24dd04ec 100644 --- a/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp +++ b/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp @@ -36,14 +36,16 @@ #ifndef OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER #define OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER +#include + +#include + #include #include #include - #include #include -#include namespace Opm { namespace ImplicitTransportDefault { diff --git a/opm/core/transport/implicit/GravityColumnSolver.hpp b/opm/core/transport/implicit/GravityColumnSolver.hpp deleted file mode 100644 index 12132978..00000000 --- a/opm/core/transport/implicit/GravityColumnSolver.hpp +++ /dev/null @@ -1,66 +0,0 @@ -/* - 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_GRAVITYCOLUMNSOLVER_HEADER_INCLUDED -#define OPM_GRAVITYCOLUMNSOLVER_HEADER_INCLUDED - -#include -#include -#include - -namespace Opm -{ - - /// Class for doing gravity segregation (only), - /// on a vertical column of cells. - template - class GravityColumnSolver - { - public: - /// Note: the model will be changed since it stores computed - /// quantities in itself, such as mobilities. - GravityColumnSolver(Model& model, - const UnstructuredGrid& grid, - const double tol, - const int maxit); - - /// \param[in] columns for each column col, columns[col] - /// contains the cells on which to solve the segregation - /// problem. For each column, its cells must be in a single - /// vertical column, connected and ordered - /// (direction doesn't matter). - void solve(const std::vector >& columns, - const double dt, - std::vector& s); - private: - void solveSingleColumn(const std::vector& columns, - const double dt, - std::vector& s, - std::vector& sol_vec); - Model& model_; - const UnstructuredGrid& grid_; - const double tol_; - const int maxit_; -}; - -} // namespace Opm - -#include - -#endif // OPM_GRAVITYCOLUMNSOLVER_HEADER_INCLUDED diff --git a/opm/core/transport/implicit/GravityColumnSolver_impl.hpp b/opm/core/transport/implicit/GravityColumnSolver_impl.hpp deleted file mode 100644 index 8e404d06..00000000 --- a/opm/core/transport/implicit/GravityColumnSolver_impl.hpp +++ /dev/null @@ -1,194 +0,0 @@ -/* - 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 . -*/ - -#include -#include -#include - -#include - -#include - -namespace Opm -{ - - template - GravityColumnSolver::GravityColumnSolver(Model& model, - const UnstructuredGrid& grid, - const double tol, - const int maxit) - : model_(model), grid_(grid), tol_(tol), maxit_(maxit) - { - } - - namespace { - struct ZeroVec - { - double operator[](int) const { return 0.0; } - }; - struct StateWithZeroFlux - { - StateWithZeroFlux(std::vector& s) : sat(s) {} - const ZeroVec& faceflux() const { return zv; } - const std::vector& saturation() const { return sat; } - std::vector& saturation() { return sat; } - ZeroVec zv; - std::vector& sat; - }; - struct Vecs - { - Vecs(int sz) : sol(sz, 0.0) {} - const std::vector& solution() const { return sol; } - std::vector& writableSolution() { return sol; } - std::vector sol; - }; - struct JacSys - { - JacSys(int sz) : v(sz) {} - const Vecs& vector() const { return v; } - Vecs& vector() { return v; } - Vecs v; - typedef std::vector vector_type; - }; - } // anon namespace - - - - /// \param[in] columns for each column columns.second, - /// contains the cells on which to solve the segregation - /// problem. For each column, its cells must be in a single - /// vertical column, and ordered - /// (direction doesn't matter). - template - void GravityColumnSolver::solve(const std::vector >& columns, - const double dt, - std::vector& s) - { - // Initialize model. These things are done for the whole grid! - StateWithZeroFlux state(s); // This holds s by reference. - JacSys sys(grid_.number_of_cells); - std::vector increment(grid_.number_of_cells, 0.0); - model_.initStep(state, grid_, sys); - - int iter = 0; - double max_delta = 1e100; - - while (iter < maxit_) { - model_.initIteration(state, grid_, sys); - // std::map >::const_iterator it; - //for (it = columns.begin(); it != columns.end(); ++it) { - int size = columns.size(); - - #pragma omp parallel for schedule(dynamic) - for(int i = 0; i < size; ++i) { - solveSingleColumn(columns[i], dt, s, increment); - } - - for (int cell = 0; cell < grid_.number_of_cells; ++cell) { - sys.vector().writableSolution()[cell] += increment[cell]; - } - - const double maxelem = *std::max_element(increment.begin(), increment.end()); - const double minelem = *std::min_element(increment.begin(), increment.end()); - max_delta = std::max(maxelem, -minelem); - std::cout << "Gravity column solver iteration " << iter << " max_delta = " << max_delta << std::endl; - if (max_delta < tol_) { - break; - } - ++iter; - } - if (max_delta >= tol_) { - OPM_THROW(std::runtime_error, "Failed to converge!"); - } - // Finalize. - // model_.finishIteration(); // Doesn't do anything in th 2p model. - // finishStep() writes to state, which holds s by reference. - // This will update the entire grid's state... - model_.finishStep(grid_, sys.vector().solution(), state); - } - - - - /// \param[in] column_cells the cells on which to solve the segregation - /// problem. Must be in a single vertical column, - /// and ordered (direction doesn't matter). - template - void GravityColumnSolver::solveSingleColumn(const std::vector& column_cells, - const double dt, - std::vector& s, - std::vector& sol_vec) - { - // This is written only to work with SinglePointUpwindTwoPhase, - // not with arbitrary problem models. - const int col_size = column_cells.size(); - if (col_size == 1) { - sol_vec[column_cells[0]] = 0.0; - return; - } - StateWithZeroFlux state(s); // This holds s by reference. - - // Assemble. - std::vector tridiag_matrix_data(3*col_size - 2, 0.0); - double* DU = &tridiag_matrix_data[0]; - double* D = DU + col_size - 1; - double* DL = D + col_size; - std::vector rhs(col_size, 0.0); - for (int ci = 0; ci < col_size; ++ci) { - double rescontrib, j1contrib, j2contrib; - const int cell = column_cells[ci]; - const int prev_cell = (ci == 0) ? -999 : column_cells[ci - 1]; - const int next_cell = (ci == col_size - 1) ? -999 : column_cells[ci + 1]; - // model_.initResidual(cell, F); - for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) { - const int face = grid_.cell_faces[j]; - const int c1 = grid_.face_cells[2*face + 0]; - const int c2 = grid_.face_cells[2*face + 1]; - if (c1 == prev_cell || c2 == prev_cell || c1 == next_cell || c2 == next_cell) { - j1contrib = j2contrib = rescontrib = 0.0; - model_.fluxConnection(state, grid_, dt, cell, face, &j1contrib, &j2contrib, &rescontrib); - if (c1 == prev_cell || c2 == prev_cell) { - DL[ci-1] += j2contrib; - } else { - assert(c1 == next_cell || c2 == next_cell); - DU[ci] += j2contrib; - } - D[ci] += j1contrib; - rhs[ci] += rescontrib; - } - } - j1contrib = rescontrib = 0.0; - model_.accumulation(grid_, cell, &j1contrib, &rescontrib); - D[ci] += j1contrib; - rhs[ci] += rescontrib; - } - // model_.sourceTerms(); // Not needed - // Solve. - const MAT_SIZE_T num_rhs = 1, colSize = col_size; - MAT_SIZE_T info = 0; - // Solution will be written to rhs. - dgtsv_(&colSize, &num_rhs, DL, D, DU, &rhs[0], &colSize, &info); - if (info != 0) { - OPM_THROW(std::runtime_error, "Lapack reported error in dgtsv: " << info); - } - for (int ci = 0; ci < col_size; ++ci) { - sol_vec[column_cells[ci]] = -rhs[ci]; - } - } - -} // namespace Opm diff --git a/opm/core/utility/Factory.hpp b/opm/core/utility/Factory.hpp index e9ca8f46..84bf803d 100644 --- a/opm/core/utility/Factory.hpp +++ b/opm/core/utility/Factory.hpp @@ -35,6 +35,8 @@ #ifndef OPM_FACTORY_HEADER #define OPM_FACTORY_HEADER +#include + #include #include diff --git a/opm/core/utility/parameters/ParameterGroup_impl.hpp b/opm/core/utility/parameters/ParameterGroup_impl.hpp index 8de64958..c927ca59 100644 --- a/opm/core/utility/parameters/ParameterGroup_impl.hpp +++ b/opm/core/utility/parameters/ParameterGroup_impl.hpp @@ -39,6 +39,7 @@ #include #include +#include #include #include #include