diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2779f49b..b8ae581f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -69,8 +69,6 @@ list (APPEND MAIN_SOURCE_FILES opm/core/simulator/TwophaseState.cpp opm/core/simulator/SimulatorReport.cpp opm/core/transport/TransportSolverTwophaseInterface.cpp - opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp - opm/core/transport/implicit/transport_source.c opm/core/transport/reorder/ReorderSolverInterface.cpp opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp @@ -251,17 +249,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/simulator/initStateEquil_impl.hpp opm/core/simulator/initState_impl.hpp opm/core/transport/TransportSolverTwophaseInterface.hpp - opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp - opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp - opm/core/transport/implicit/ImplicitAssembly.hpp - opm/core/transport/implicit/ImplicitTransport.hpp - opm/core/transport/implicit/JacobianSystem.hpp - opm/core/transport/implicit/NormSupport.hpp - opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp - opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp - opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp - opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp - opm/core/transport/implicit/transport_source.h opm/core/transport/reorder/ReorderSolverInterface.hpp opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp diff --git a/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp b/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp deleted file mode 100644 index 24dd04ec..00000000 --- a/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp +++ /dev/null @@ -1,198 +0,0 @@ -/*=========================================================================== -// -// File: CSRMatrixBlockAssembler.hpp -// -// Created: 2011-10-03 12:40:56+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER -#define OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER - -#include - -#include - -#include -#include -#include -#include -#include - - -namespace Opm { - namespace ImplicitTransportDefault { - - template <> - class MatrixZero { - public: - static void - zero(struct CSRMatrix& A) { - csrmatrix_zero(&A); - } - }; - - template <> - class MatrixBlockAssembler { - public: - template - void - assembleBlock(::std::size_t ndof, - ::std::size_t i , - ::std::size_t j , - const Block& b ) { - - assert (ndof > 0); - assert (ndof == ndof_); - - const ::std::size_t start = ia_[i*ndof + 0]; - const ::std::size_t off = - csrmatrix_elm_index(i * ndof, j * ndof, &mat_) - start; - - for (::std::size_t row = 0; row < ndof; ++row) { - const ::std::size_t J = ia_[i*ndof + row] + off; - - for (::std::size_t col = 0; col < ndof; ++col) { - sa_[J + col] += b[col*ndof + row]; - } - } - } - - template - void - createBlockRow(::std::size_t i , - const Connections& conn, - ::std::size_t ndof) { - - assert (ndof > 0); - assert (ndof == ndof_); - assert (i == (ia_.size() - 1) / ndof_); (void) i; - - expandSortConn(conn, ndof); - const int nconn = static_cast(esconn_.size()); - - for (::std::size_t dof = 0; dof < ndof; ++dof) { - ja_.insert(ja_.end(), esconn_.begin(), esconn_.end()); - ia_.push_back(ia_.back() + nconn); - } - - sa_.insert(sa_.end(), nconn * ndof, double(0.0)); - - finalizeStructure(); - } - - void - finalizeStructure() { - construct(); - setCSRSize(); - } - - void - setSize(size_t ndof, size_t m, size_t n, size_t nnz = 0) { - (void) n; - - clear(); - - allocate(ndof, m, nnz); - - ia_.push_back(0); - ndof_ = ndof; - } - - struct CSRMatrix& matrix() { return mat_; } - const struct CSRMatrix& matrix() const { return mat_; } - - private: - void - allocate(::std::size_t ndof, ::std::size_t m, ::std::size_t nnz) { - ia_.reserve(1 + ( m * ndof)); - ja_.reserve(0 + (nnz * ndof)); - sa_.reserve(0 + (nnz * ndof)); - } - - void - clear() { - ia_.resize(0); - ja_.resize(0); - sa_.resize(0); - } - - void - construct() { - mat_.ia = &ia_[0]; - mat_.ja = &ja_[0]; - mat_.sa = &sa_[0]; - } - - template - void - expandSortConn(const Connections& conn, ::std::size_t ndof) { - sconn_.resize(0); - sconn_.reserve(conn.size()); - - for (typename Connections::const_iterator - c = conn.begin(), e = conn.end(); c != e; ++c) { - sconn_.push_back(static_cast(*c)); - } - - ::std::sort(sconn_.begin(), sconn_.end()); - - esconn_.resize(0); - esconn_.reserve(ndof * sconn_.size()); - - for (::std::vector::iterator - c = sconn_.begin(), e = sconn_.end(); c != e; ++c) { - for (::std::size_t dof = 0; dof < ndof; ++dof) { - esconn_.push_back(static_cast((*c)*ndof + dof)); - } - } - } - - void - setCSRSize() { - mat_.m = ia_.size() - 1; - mat_.nnz = ja_.size() ; - } - - ::std::size_t ndof_; - - ::std::vector ia_; - ::std::vector ja_; - ::std::vector sa_; - - ::std::vector sconn_ ; - ::std::vector esconn_; - - struct CSRMatrix mat_; - }; - } -} - -#endif /* OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER */ diff --git a/opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp b/opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp deleted file mode 100644 index 518ea222..00000000 --- a/opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp +++ /dev/null @@ -1,87 +0,0 @@ -/*=========================================================================== -// -// File: CSRMatrixUmfpackSolver.hpp -// -// Created: 2011-10-03 17:27:26+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_CSRMATRIXUMFPACKSOLVER_HPP_HEADER -#define OPM_CSRMATRIXUMFPACKSOLVER_HPP_HEADER - -#include -#include - -namespace Opm -{ - namespace ImplicitTransportLinAlgSupport - { - - class CSRMatrixUmfpackSolver - { - public: - - - template - void - solve(const struct CSRMatrix* A, - const Vector b, - Vector x) - { -#if HAVE_SUITESPARSE_UMFPACK_H - call_UMFPACK(const_cast(A), b, x); -#else - OPM_THROW(std::runtime_error, "Cannot use implicit transport solver without UMFPACK. " - "Reconfigure opm-core with SuiteSparse/UMFPACK support and recompile."); -#endif - } - - - template - void - solve(const struct CSRMatrix& A, - const Vector& b, - Vector& x) - { -#if HAVE_SUITESPARSE_UMFPACK_H - call_UMFPACK(const_cast(&A), &b[0], &x[0]); -#else - OPM_THROW(std::runtime_error, "Cannot use implicit transport solver without UMFPACK. " - "Reconfigure opm-core with SuiteSparse/UMFPACK support and recompile."); -#endif - } - - - }; // class CSRMatrixUmfpackSolver - - } // namespace ImplicitTransportLinAlgSupport -} // namespace Opm - -#endif /* OPM_CSRMATRIXUMFPACKSOLVER_HPP_HEADER */ diff --git a/opm/core/transport/implicit/ImplicitAssembly.hpp b/opm/core/transport/implicit/ImplicitAssembly.hpp deleted file mode 100644 index 705af102..00000000 --- a/opm/core/transport/implicit/ImplicitAssembly.hpp +++ /dev/null @@ -1,267 +0,0 @@ -/*=========================================================================== -// -// File: ImplicitAssembly.hpp -// -// Created: 2011-09-28 10:00:46+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_IMPLICITASSEMBLY_HPP_HEADER -#define OPM_IMPLICITASSEMBLY_HPP_HEADER - -#include -#include -#include - -namespace Opm { - template - class ImplicitAssembly { - enum { DofPerCell = Model::DofPerCell }; - - public: - ImplicitAssembly(Model& model) - : model_(model), - nconn_(-1) , - asm_buffer_() - {} - - template - void - createSystem(const Grid& g ,//const SourceTerms* tsrc, - JacobianSystem& sys) const { - - ::std::size_t m = g.number_of_cells; - ::std::size_t nnz = g.number_of_cells - + countConnections(g); -// + tsrc->pf; - - sys.setSize(DofPerCell, m, nnz); - - for (int c = 0; c < g.number_of_cells; ++c) { - this->createRowStructure(g, c, sys); - } - /* is done by modifiying grid - // add elements for periodic boundary - for (int pf = 0; pfpf; ++pf){ - int nconn = 1; - ::std::vector connections(1); - connections[0]=(tsrc->periodic_cells[2*pf]); - sys.matasm().createBlockRow(tsrc->periodic_cells[2*pf+1], connections, DofPerCell); - connections[0]=(tsrc->periodic_cells[2*pf+1]); - sys.matasm().createBlockRow(tsrc->periodic_cells[2*pf], connections, DofPerCell); - } - */ - sys.matasm().finalizeStructure(); - } - - template - void - assemble(const ReservoirState& state, - const Grid& g , - const SourceTerms* src , - const double dt , - JacobianSystem& sys ) { - - for (int c = 0; c < g.number_of_cells; ++c) { - this->computeCellContrib(state, g, dt, c); - this->assembleCellContrib(g, c, sys); - } - - if (src != 0) { - this->assembleSourceContrib(g, src, dt, sys); - } - } - - private: - template - int - countConnections(const Grid& g, int c) const { - int i, f, c1, c2, n; - - for (i = g.cell_facepos[c + 0], n = 0; - i < g.cell_facepos[c + 1]; ++i) { - f = g.cell_faces[i]; - c1 = g.face_cells[2*f + 0]; - c2 = g.face_cells[2*f + 1]; - - n += (c1 >= 0) && (c2 >= 0); - } - - return n; - } - - template - int - countConnections(const Grid& g) const { - int n = 0; - - for (int c = 0; c < g.number_of_cells; ++c) { - n += countConnections(g, c); - } - - return n; - } - - template - void - createRowStructure(const Grid& g , - const int c , - System& sys) const { - - int nconn = countConnections(g, c); - - ::std::vector connections; - connections.reserve (nconn + 1); - connections.push_back(c); - - for (int i = g.cell_facepos[c + 0]; - i < g.cell_facepos[c + 1]; ++i) { - int f = g.cell_faces[i]; - int c1 = g.face_cells[2*f + 0]; - int c2 = g.face_cells[2*f + 1]; - - if ((c1 >= 0) && (c2 >= 0)) { - connections.push_back((c1 == c) ? c2 : c1); - } - } - sys.matasm().createBlockRow(c, connections, DofPerCell); - - } - - template - int - computeCellContrib(const ReservoirState& state, - const Grid& g , - const double dt , - const int c ) { - const int ndof = DofPerCell; - const int ndof2 = ndof * ndof; - - nconn_ = countConnections(g, c); - - asm_buffer_.resize((2*nconn_ + 1)*ndof2 + (nconn_ + 2)*ndof); - std::fill(asm_buffer_.begin(), asm_buffer_.end(), 0.0); - - double* F = &asm_buffer_[(2*nconn_ + 1) * ndof2]; - double* J1 = &asm_buffer_[(0*nconn_ + 1) * ndof2]; - double* J2 = J1 + (1*nconn_ + 0) * ndof2 ; - - model_.initResidual(c, F); - F += ndof; - - for (int i = g.cell_facepos[c + 0]; - i < g.cell_facepos[c + 1]; ++i) { - int f = g.cell_faces[i]; - int c1 = g.face_cells[2*f + 0]; - int c2 = g.face_cells[2*f + 1]; - - if ((c1 >= 0) && (c2 >= 0)) { - model_.fluxConnection(state, g, dt, c, f, J1, J2, F); - J1 += ndof2; J2 += ndof2; F += ndof; - } - } - - model_.accumulation(g, c, &asm_buffer_[0], F); - - return 1; - } - - template - void - assembleCellContrib(const Grid& g , - const int c , - System& sys) const { - const int ndof = DofPerCell; - const int ndof2 = ndof * ndof; - - const double* J1 = &asm_buffer_[0]; - const double* J2 = J1 + ((1*nconn_ + 1) * ndof2); - - // Assemble contributions from accumulation term - sys.matasm().assembleBlock(ndof, c, c, J1); J1 += ndof2; - - // Assemble connection contributions. - for (int i = g.cell_facepos[c + 0]; - i < g.cell_facepos[c + 1]; ++i) { - int f = g.cell_faces[i]; - int c1 = g.face_cells[2*f + 0]; - int c2 = g.face_cells[2*f + 1]; - - c2 = (c1 == c) ? c2 : c1; - - if (c2 >= 0) { - sys.matasm().assembleBlock(ndof, c, c , J1); - sys.matasm().assembleBlock(ndof, c, c2, J2); - - J1 += ndof2; - J2 += ndof2; - } - } - // Assemble residual - const double* F = &asm_buffer_[(2*nconn_ + 1) * ndof2]; - for (int conn = 0; conn < nconn_ + 2; ++conn, F += ndof) { - sys.vector().assembleBlock(ndof, c, F); - } - } - - template - void - assembleSourceContrib(const Grid& g, - const SourceTerms* src, - const double dt, - System& sys) { - const int ndof = DofPerCell; - const int ndof2 = ndof * ndof; - - for (int i = 0; i < src->nsrc; ++i) { - std::fill_n(asm_buffer_.begin(), ndof2 + ndof, 0.0); - - double *J = &asm_buffer_[0]; - double *F = J + ndof2; - - model_.sourceTerms(g, src, i, dt, J, F); - - const int c = src->cell[i]; - - sys.matasm().assembleBlock(ndof, c, c, J); - sys.vector().assembleBlock(ndof, c, F); - } - } - - Model& model_ ; - int nconn_ ; - std::vector asm_buffer_; - }; -} -#endif /* OPM_IMPLICITASSEMBLY_HPP_HEADER */ diff --git a/opm/core/transport/implicit/ImplicitTransport.hpp b/opm/core/transport/implicit/ImplicitTransport.hpp deleted file mode 100644 index b775575d..00000000 --- a/opm/core/transport/implicit/ImplicitTransport.hpp +++ /dev/null @@ -1,234 +0,0 @@ -/*=========================================================================== -// -// File: ImplicitTransport.hpp -// -// Created: 2011-09-29 10:38:42+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_IMPLICITTRANSPORT_HPP_HEADER -#define OPM_IMPLICITTRANSPORT_HPP_HEADER - -#include - -#include - -namespace Opm { - namespace ImplicitTransportDetails { - struct NRControl { - NRControl() - : max_it(1), - atol(1.0e-6), - rtol(5.0e-7), - dxtol(1.0e-8), - max_it_ls(5), - verbosity(0) - {} - - int max_it; - double atol ; - double rtol ; - double dxtol ; - int max_it_ls; - int verbosity; - }; - - struct NRReport { - int nit; - int flag; - double norm_res; - double norm_dx; - }; - - template - Ostream& operator<<(Ostream& os, const NRReport& rpt) - { - os << "Number of linear solves: " << rpt.nit << '\n' - << "Process converged: " << (rpt.flag > 0) << '\n' - << "Convergence flag: " << rpt.flag; - switch (rpt.flag) { - case 1: - os << " (absolute tolerance satisfied)\n"; - break; - case 2: - os << " (relative tolerance satisfied)\n"; - break; - case 3: - os << " (solution change tolerance satisfied)\n"; - break; - case -1: - os << " (failed to converge)\n"; - break; - default: - os << " (unknown convergence flag)\n"; - } - os << "Final residual norm: " << rpt.norm_res << '\n' - << "Final increment norm: " << rpt.norm_dx << '\n'; - return os; - } - } - - template class VNorm , - template class VNeg , - template class VZero , - template class MZero , - template class VAsgn > - class ImplicitTransport { - public: - ImplicitTransport(Model& model) - : model_(model), - asm_ (model) - {} - - template - void solve(const Grid& g , - const SourceTerms* src , - const double dt , - const ImplicitTransportDetails::NRControl& ctrl , - ReservoirState& state , - LinearSolver& linsolve, - ImplicitTransportDetails::NRReport& rpt ) { - bool init; - typedef typename JacobianSystem::vector_type vector_type; - typedef typename JacobianSystem::matrix_type matrix_type; - - asm_.createSystem(g, sys_); - model_.initStep(state, g, sys_); - init = model_.initIteration(state, g, sys_); - - MZero::zero(sys_.writableMatrix()); - VZero::zero(sys_.vector().writableResidual()); - - asm_.assemble(state, g, src, dt, sys_); - - const double nrm_res0 = - VNorm::norm(sys_.vector().residual()); - - rpt.norm_res = nrm_res0; - rpt.norm_dx = -1.0; - rpt.nit = 0; - - bool done = false; //rpt.norm_res < ctrl.atol; - - while (! done) { - VZero::zero(sys_.vector().writableIncrement()); - - linsolve.solve(sys_.matrix(), - sys_.vector().residual(), - sys_.vector().writableIncrement()); - - VNeg::negate(sys_.vector().writableIncrement()); - - //model_.finishIteration(state, g, sys_.vector()); - - rpt.norm_dx = - VNorm::norm(sys_.vector().increment()); - - // Begin line search - double residual=VNorm::norm(sys_.vector().residual()); - int lin_it=0; - bool finished=rpt.norm_res::assign(alpha, dx_old, - sys_.vector().writableIncrement()); - VAsgn::assign(x_old, - sys_.vector().writableSolution()); - - sys_.vector().addIncrement(); - init = model_.initIteration(state, g, sys_); - if (init) { - MZero::zero(sys_.writableMatrix()); - VZero::zero(sys_.vector().writableResidual()); - asm_.assemble(state, g, src, dt, sys_); - residual = VNorm::norm(sys_.vector().residual()); - if (ctrl.verbosity > 1){ - std::cout << "Line search iteration " << std::scientific << lin_it - << " norm :" << residual << " alpha " << alpha << '\n'; - } - }else{ - if (ctrl.verbosity > 1){ - std::cout << "Line search iteration " << std::scientific << lin_it - << " Value out of range, continue search. alpha " << alpha << '\n'; - } - residual = 1e99; - } - finished=(residual < rpt.norm_res) || (lin_it> ctrl.max_it_ls); - lin_it +=1; - }// End line search - rpt.norm_res = - VNorm::norm(sys_.vector().residual()); - - rpt.nit++; - if (ctrl.verbosity > 0){ - std::cout << "Iteration " << std::scientific << rpt.nit - << " norm :" << rpt.norm_res << " alpha " << alpha << std::endl; - } - done = (rpt.norm_res < ctrl.atol) || - (rpt.norm_res < ctrl.rtol * nrm_res0) || - (rpt.norm_dx < ctrl.dxtol) || - (lin_it > ctrl.max_it_ls) || - (rpt.nit == ctrl.max_it); - } - - model_.finishStep(g, sys_.vector().solution(), state); - - if (rpt.norm_res < ctrl.atol) { rpt.flag = 1; } - else if (rpt.norm_res < ctrl.rtol * nrm_res0) { rpt.flag = 2; } - else if (rpt.norm_dx < ctrl.dxtol) { rpt.flag = 3; } - else { rpt.flag = -1; } - } - - private: - ImplicitTransport (const ImplicitTransport&); - ImplicitTransport& operator=(const ImplicitTransport&); - -#if 0 - using Model::initStep; - using Model::initIteration; - using Model::finishIteration; - using Model::finishStep; -#endif - - Model& model_; - ImplicitAssembly asm_; - JacobianSystem sys_; - }; -} -#endif /* OPM_IMPLICITTRANSPORT_HPP_HEADER */ diff --git a/opm/core/transport/implicit/JacobianSystem.hpp b/opm/core/transport/implicit/JacobianSystem.hpp deleted file mode 100644 index 400bfb2a..00000000 --- a/opm/core/transport/implicit/JacobianSystem.hpp +++ /dev/null @@ -1,256 +0,0 @@ -/*=========================================================================== -// -// File: JacobianSystem.hpp -// -// Created: 2011-09-30 19:23:31+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_JACOBIANSYSTEM_HPP_HEADER -#define OPM_JACOBIANSYSTEM_HPP_HEADER - -#include -#include -#include - -#include -//#include -#include -#include - -namespace Opm { - namespace ImplicitTransportDefault { - template - class VectorAdder { - public: - // y += x - static void - add(const BaseVec& x, BaseVec& y) { - typedef typename BaseVec::value_type VT; - ::std::transform(x.begin(), x.end(), - y.begin(), - y.begin(), - ::std::plus()); - } - }; - - template - class VectorNegater { - public: - // x *= -1 - static void - negate(BaseVec& x) { - typedef typename BaseVec::value_type VT; - ::std::transform(x.begin(), x.end(), - x.begin(), - ::std::negate()); - } - }; - - template - class VectorZero { - public: - static void - zero(BaseVec& x) { - typedef typename BaseVec::value_type VT; - - ::std::fill(x.begin(), x.end(), VT(0.0)); - } - }; - - template - class VectorAssign { - public: - // y <- x - static void - assign(const BaseVec& x, BaseVec& y) { - ::std::copy(x.begin(), x.end(), y.begin()); - } - - // y <- a*x - template - static void - assign(const Scalar& a, const BaseVec& x, BaseVec& y) { - typedef typename BaseVec::value_type VT; - ::std::transform(x.begin(), x.end(), - y.begin(), - ::std::bind2nd(::std::multiplies(), a)); - } - }; - - template - class MatrixZero; - - template - class VectorBlockAssembler { - public: - template - static void - assemble(::std::size_t ndof, - ::std::size_t i , - const Block& b , - BaseVec& vec ) { - - for (::std::size_t d = 0; d < ndof; ++d) { - vec[i*ndof + d] += b[d]; - } - } - }; - - template - class VectorSizeSetter { - public: - VectorSizeSetter(BaseVec& v) : v_(v) {} - - void - setSize(::std::size_t ndof, ::std::size_t m) { - v_.resize(ndof * m); - } - - private: - BaseVec& v_; - }; - - template class VSzSetter = VectorSizeSetter , - template class VAdd = VectorAdder , - template class VBlkAsm = VectorBlockAssembler> - class NewtonVectorCollection { - enum { Residual = 0, Increment = 1, Solution = 2 }; - - public: - void - setSize(::std::size_t ndof, ::std::size_t m) { -#if 0 - typedef typename ::std::array::iterator VAI; - - for (VAI i = vcoll_.begin(), e = vcoll_.end(); i != e; ++i) { - VSzSetter(*i).setSize(ndof, m); - } -#endif - for (::std::size_t i = 0; i < sizeof (vcoll_) / sizeof (vcoll_[0]); ++i) { - VSzSetter(vcoll_[i]).setSize(ndof, m); - } - - ndof_ = ndof; - } - - void - addIncrement() { - VAdd::add(vcoll_[ Increment ], vcoll_[ Solution ]); - } - - template - void - assembleBlock(::std::size_t ndof, - ::std::size_t i , - const Block& b ) { - - assert (ndof_ > 0 ); - assert (ndof == ndof_); - - VBlkAsm::assemble(ndof, i, b, vcoll_[ Residual ]); - } - - typedef BaseVec vector_type; - - const vector_type& increment() const { return vcoll_[ Increment ]; } - const vector_type& residual () const { return vcoll_[ Residual ]; } - const vector_type& solution () const { return vcoll_[ Solution ]; } - - // Write access for Newton solver purposes - vector_type& writableIncrement() { return vcoll_[ Increment ]; } - vector_type& writableResidual () { return vcoll_[ Residual ]; } - vector_type& writableSolution () { return vcoll_[ Solution ]; } - - private: - ::std::size_t ndof_ ; - BaseVec vcoll_[3]; - //::std::array vcoll_; - }; - - template - class MatrixBlockAssembler - /* { - public: - template - void - assembleBlock(size_t n, size_t i, size j, const Block& b); - - template - void - createBlockRow(size_t i, const Connections& conn, size_t n); - - void - finalizeStructure(); - - void - setSize(size_t ndof, size_t m, size_t n, size_t nnz = 0); - - Matrix& matrix() ; - const Matrix& matrix() const; - } */; - - template - class JacobianSystem { - public: - JacobianSystem() {} - - typedef Matrix matrix_type; - typedef MatrixBlockAssembler assembler_type; - typedef NVecCollection vector_collection_type; - typedef typename NVecCollection::vector_type vector_type; - - assembler_type& matasm() { return mba_ ; } - NVecCollection& vector() { return sysvec_ ; } - const matrix_type& matrix() const { return mba_.matrix(); } - matrix_type& writableMatrix() { return mba_.matrix(); } - - void - setSize(::std::size_t ndof, - ::std::size_t m , - ::std::size_t nnz = 0) { - - mba_ .setSize(ndof, m, m, nnz); - sysvec_.setSize(ndof, m ); - } - - private: - JacobianSystem (const JacobianSystem&); - JacobianSystem& operator=(const JacobianSystem&); - - MatrixBlockAssembler mba_ ; // Coefficient matrix - NVecCollection sysvec_; // Residual - }; - } -} - -#endif /* OPM_JACOBIANSYSTEM_HPP_HEADER */ diff --git a/opm/core/transport/implicit/NormSupport.hpp b/opm/core/transport/implicit/NormSupport.hpp deleted file mode 100644 index 4bc29b8a..00000000 --- a/opm/core/transport/implicit/NormSupport.hpp +++ /dev/null @@ -1,100 +0,0 @@ -/*=========================================================================== -// -// File: NormSupport.hpp -// -// Created: 2011-10-04 19:37:35+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_NORMSUPPORT_HPP_HEADER -#define OPM_NORMSUPPORT_HPP_HEADER - -#include -#include -#include -#include - -namespace Opm { - namespace ImplicitTransportDefault { - template - class MaxAbs : public ::std::binary_function { - public: - double - operator()(double x, const T& y) { - return std::max(std::abs(x), std::abs(y)); - } - - static double - postprocess(double nrm_inf) { return nrm_inf; } - }; - - template - class SumAbs : public ::std::binary_function { - public: - double - operator()(double x, const T& y) { - return std::abs(x) + std::abs(y); - } - - static double - postprocess(double nrm_1) { return nrm_1; } - }; - - template - class Euclid : public ::std::binary_function { - public: - double - operator()(double x, const T& y) { - const double ay = ::std::abs(y); - - return std::abs(x) + ay*ay; - } - - static double - postprocess(double nrm2) { return ::std::sqrt(nrm2); } - }; - - template class NormImpl> - class AccumulationNorm { - public: - static double - norm(const Vector& v) { - typedef typename Vector::value_type VT; - - double nrm = ::std::accumulate(v.begin(), v.end(), VT(0.0), - NormImpl()); - - return NormImpl::postprocess(nrm); - } - }; - } -} - -#endif /* OPM_NORMSUPPORT_HPP_HEADER */ diff --git a/opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp b/opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp deleted file mode 100644 index 3bf8fcf9..00000000 --- a/opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp +++ /dev/null @@ -1,70 +0,0 @@ -/*=========================================================================== -// -// 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/implicit/SimpleFluid2pWrappingProps_impl.hpp b/opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp deleted file mode 100644 index fa8a5c3b..00000000 --- a/opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp +++ /dev/null @@ -1,104 +0,0 @@ -/*=========================================================================== -// -// 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) { - OPM_THROW(std::runtime_error, "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/implicit/SinglePointUpwindTwoPhase.hpp b/opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp deleted file mode 100644 index 741331a0..00000000 --- a/opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp +++ /dev/null @@ -1,725 +0,0 @@ -/*=========================================================================== -// -// File: SinglePointUpwindTwoPhase.hpp -// -// Created: 2011-09-28 14:21:34+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -/** - * \file - * Numerical model and support classes needed to model transport of two - * incompressible fluid phases. Intended for the ImplicitTransport system. - */ - -#ifndef OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER -#define OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER - -#include -#include - -#include -#include -#include - -namespace Opm { - namespace spu_2p { - /** - * Internal class to manage the direct and derived quantities needed to - * formulate the fluid transport system. - * - * Note: This class elides off-diagonal elements of the phase mobility - * Jacobian, \f$(\partial_{s_\beta} \lambda_\alpha)_{\alpha\beta}\f$. - * These elements are assumed to be strictly equal to zero. In other - * words, the relative permeability of phase \f$\alpha\f$ is assumed to - * depend only on the saturation of phase \f$\alpha\f$. This convention - * allows storing only the two diagonals of the mobility Jacobian per - * grid cell. - * - * The static gravity term is the scalar value - * \f[ - * \Delta G_i = \mathsf{T}_f\, \vec{g}\cdot(\Bar{x}_f - \Bar{x}_c) - * \f] - * in which @c i is the half face index corresponding to the cell-face - * pair (f,c) and \f$\mathsf{T}_f\f$ is the absolute - * (bacground) two-point transmissibility of face @c f. - * - * The fluid transport problem is formulated in terms of saturation - * changes, \f$\Delta s\f$, per cell. These changes are the primary - * degrees of freedom in this model. - * - * Capillary pressures are defined by the fluid model, but usually - * correspond to \f$p_w - p_n\f$ (e.g., \f$p_\mathit{oil} - - * p_\mathit{water}\f$). - */ - class ModelParameterStorage { - public: - /** - * Constructor. - * - * @param[in] nc Total number of grid cells. - * @param[in] totconn Total number of connections, accumulated per - * cell (``half faces''). - */ - ModelParameterStorage(int nc, int totconn) - : drho_(0.0), mob_(0), dmob_(0), - porevol_(0), dg_(0), ds_(0), pc_(0), dpc_(0), trans_(0), - data_() - { - std::size_t alloc_sz; - - alloc_sz = 2 * nc; // mob_ - alloc_sz += 2 * nc; // dmob_ - alloc_sz += 1 * nc; // porevol_ - alloc_sz += 1 * totconn; // dg_ - alloc_sz += 1 * nc; // ds_ - alloc_sz += 1 * nc; // pc_ - alloc_sz += 1 * nc; // dpc_ - alloc_sz += 1 * totconn; // dtrans - data_.resize(alloc_sz); - - mob_ = &data_[0]; - dmob_ = mob_ + (2 * nc ); - porevol_ = dmob_ + (2 * nc ); - dg_ = porevol_ + (1 * nc ); - ds_ = dg_ + (1 * totconn); - pc_ = ds_ + (1 * nc ); - dpc_ = pc_ + (1 * nc ); - trans_ = dpc_ + (1 * nc ); - } - - /** - * Modifiable density difference. - * @return Reference to modifiable internal representation of fluid - * phase density difference. - */ - double& drho () { return drho_ ; } - /** - * Read-only density difference. - * @return Read-only value of current fluid phase difference value. - */ - double drho () const { return drho_ ; } - - /** - * Phase mobility in cell. - * @param[in] c Cell. - * @return Read-write reference to two consecutive phase mobilities - * in cell @c c. - */ - double* mob (int c) { return mob_ + (2*c + 0); } - /** - * Phase mobility in cell. - * @param[in] c Cell. - * @return Read-only reference to two consecutive phase mobilities - * in cell @c c. - */ - const double* mob (int c) const { return mob_ + (2*c + 0); } - - /** - * Diagonal elements of phase mobility derivative (Jacobian). - * - * @param[in] c Cell. - * @return Read-write reference to diagonal elements of phase - * mobility Jacobian in cell @c c. - */ - double* dmob (int c) { return dmob_ + (2*c + 0); } - /** - * Diagonal elements of phase mobility derivative (Jacobian). - * - * @param[in] c Cell. - * @return Read-only reference to two consecutive diagonal elements - * of phase mobility Jacobian in cell @c c. - */ - const double* dmob (int c) const { return dmob_ + (2*c + 0); } - - /** - * Retrieve pore volumes for all cells. - * @return Modifiable vector of pore volumes for all cells. - */ - double* porevol() { return porevol_ ; } - /** - * Pore volume of single cell. - * @param[in] c Cell. - * @return Pore volume of cell @c c. - */ - double porevol(int c) const { return porevol_[c] ; } - - /** - * Static gravity term associated to single half face. - * - * @param[in] i Half face index corresponding to particular - * cell-face pair. - * @return Read-write reference to static gravity term of single - * half face. - */ - double& dg(int i) { return dg_[i] ; } - /** - * Static gravity term associated to single half face. - * @param[in] i Half face index corresponding to particular - * cell-face pair. - * @return Read-only reference to static gravity term of single half - * face. - */ - double dg(int i) const { return dg_[i] ; } - - /** - * Saturation change in particular cell. - * - * @param[in] c - * @return Read-write reference to saturation change (scalar) in - * cell @c c. - */ - double& ds(int c) { return ds_[c] ; } - /** - * Saturation change in particular cell. - * - * @param[in] c - * @return Read-only reference to saturation change (scalar) in cell - * @c c. - */ - double ds(int c) const { return ds_[c] ; } - - /** - * Capillary pressure in particular cell. - * - * @param[in] c Cell. - * @return Read-write reference to capillary pressure in cell @c c. - */ - double& pc(int c) { return pc_[c] ; } - /** - * Capillary pressure in particular cell. - * - * @param[in] c Cell - * @return Read-only reference to capillary pressure in cell @c c. - */ - double pc(int c) const { return pc_[c] ; } - - /** - * Derivative of capillary pressure with respect to saturation. - * - * @param[in] c Cell - * @return Read-write reference to capillary pressure derivative - * with respect to primary saturation in cell @c c. - */ - double& dpc(int c) { return dpc_[c] ; } - /** - * Derivative of capillary pressure with respect to saturation. - * - * @param[in] c Cell - * @return Read-only reference to capillary pressure derivative with - * respect to primary saturation in cell @c c. - */ - double dpc(int c) const { return dpc_[c] ; } - - /** - * Background (absolute) face transmissibility of particular face. - * - * @param[in] f Face - * @return Read-write reference to background face transmissibility - * of face @c f. - */ - double& trans(int f) { return trans_[f] ; } - /** - * Background (absolute) face transmissibility of particular face. - * - * @param[in] f Face - * @return Read-only reference to bacground face transmissibility of - * face @c f. - */ - double trans(int f) const { return trans_[f] ; } - - private: - double drho_ ; /**< Fluid phase density difference */ - double *mob_ ; /**< Fluid phase mobility in all cells */ - double *dmob_ ; /**< Derivative of phase mobility in all cells */ - double *porevol_; /**< Pore volume in all cells */ - double *dg_ ; /**< Static gravity term on all half faces */ - double *ds_ ; /**< Saturation change in all cells */ - double *pc_ ; /**< Capillary pressure in all cells */ - double *dpc_ ; /**< Derivative of cap. pressure in all cells */ - double *trans_ ; /**< Absolute transmissibility on all faces */ - - /** - * Data storage from which individual quantities are managed. - */ - std::vector data_; - }; - } - - - template - class SinglePointUpwindTwoPhase { - public: - template - SinglePointUpwindTwoPhase(const TwophaseFluid& fluid , - const Grid& g , - const std::vector& porevol , - const double* grav = 0, - const bool guess_previous = true) - : fluid_ (fluid) , - gravity_(grav) , - f2hf_ (2 * g.number_of_faces, -1) , - store_ (g.number_of_cells, - g.cell_facepos[ g.number_of_cells ]), - init_step_use_previous_sol_(guess_previous), - sat_tol_(1e-5) - { - - if (gravity_) { - store_.drho() = fluid_.density(0) - fluid_.density(1); - } - - for (int c = 0, i = 0; c < g.number_of_cells; ++c) { - for (; i < g.cell_facepos[c + 1]; ++i) { - const int f = g.cell_faces[i]; - const int p = 1 - (g.face_cells[2*f + 0] == c); - f2hf_[2*f + p] = i; - } - } - - std::copy(porevol.begin(), porevol.end(), store_.porevol()); - } - - void - makefhfQPeriodic(const std::vector& p_faces , - const std::vector& hf_faces, - const std::vector& nb_faces) - { - if (p_faces.empty()) { - return; - } - assert (p_faces.size() == hf_faces.size()); - assert (hf_faces.size() == nb_faces.size()); - - std::vector nbhf(hf_faces.size()); - - for (std::vector::size_type i = 0; i < p_faces.size(); ++i) { - const int nbf = nb_faces[i]; - - assert (2*std::vector::size_type(nbf) + 1 < f2hf_.size()); - assert ((f2hf_[2*nbf + 0] < 0) ^ (f2hf_[2*nbf + 1] < 0)); - - const int p = (f2hf_[2*nbf + 0] < 0) ? 1 : 0; // "Self" - nbhf[ i ] = f2hf_[2*nbf + p]; - } - - for (std::vector::size_type i = 0; i < p_faces.size(); ++i) { - const int f = p_faces [i]; - const int hf = hf_faces[i]; - - assert (0 <= f); - assert (0 <= hf); - assert (2*std::vector::size_type(f) + 1 < f2hf_.size()); - - assert ((f2hf_[2*f + 0] < 0 ) ^ (f2hf_[2*f + 1] < 0 )); - assert ((f2hf_[2*f + 0] == hf) ^ (f2hf_[2*f + 1] == hf)); - - const int p = (f2hf_[2*f + 0] == hf) ? 1 : 0; // "Other" - - f2hf_[2*f + p] = nbhf[ i ]; - } - } - - // ----------------------------------------------------------------- - // System assembly innards - // ----------------------------------------------------------------- - - enum { DofPerCell = 1 }; - - void - initResidual(const int c, double* F) const { - (void) c; // Suppress 'unused' warning - *F = 0.0; - } - - template - void - fluxConnection(const ReservoirState& state, - const Grid& g , - const double dt , - const int c , - const int f , - double* J1 , - double* J2 , - double* F ) const { - - const int *n = g.face_cells + (2 * f); - double dflux = state.faceflux()[f]; - double gflux = gravityFlux(f); - double pcflux,dpcflux[2]; - capFlux(f,n, pcflux, dpcflux); - gflux += pcflux; - - int pix[2]; - double m[2], dm[2]; - upwindMobility(dflux, gflux, n, pix, m, dm); - - assert (! ((m[0] < 0) || (m[1] < 0))); - - double mt = m[0] + m[1]; - assert (mt > 0); - - double sgn = 2.0*(n[0] == c) - 1.0; - dflux *= sgn; - gflux *= sgn; - - - double f1 = m[0] / mt; - const double v1 = dflux + m[1]*gflux; - - // Assemble residual contributions - *F += dt * f1 * v1; - - // Assemble Jacobian (J1 <-> c, J2 <-> other) - double *J[2]; - if (n[0] == c) { - J[0] = J1; J[1] = J2; - // sign is positive - J1[0*2 + 0] += sgn*dt * f1 * dpcflux[0] * m[1]; - J2[0*2 + 0] += sgn*dt * f1 * dpcflux[1] * m[1]; - } else { - J[0] = J2; J[1] = J1; - // sign is negative - J1[0*2 + 0] += sgn*dt * f1 * dpcflux[1] * m[1]; - J2[0*2 + 0] += sgn*dt * f1 * dpcflux[0] * m[1]; - } - - // dF/dm_1 \cdot dm_1/ds - *J[ pix[0] ] += dt * (1 - f1) / mt * v1 * dm[0]; - - /* dF/dm_2 \cdot dm_2/ds */ - *J[ pix[1] ] -= dt * f1 / mt * v1 * dm[1]; - *J[ pix[1] ] += dt * f1 * gflux * dm[1]; - } - - template - void - accumulation(const Grid& g, - const int c, - double* J, - double* F) const { - (void) g; - - const double pv = store_.porevol(c); - - *J += pv; - *F += pv * store_.ds(c); - } - - template - void - sourceTerms(const Grid& g , - const SourceTerms* src, - const int i , - const double dt , - double* J , - double* F ) const { - - (void) g; - - double dflux = -src->flux[i]; // ->flux[] is rate of *inflow* - - if (dflux < 0) { - // src -> cell, affects residual only. - *F += dt * dflux * src->saturation[2*i + 0]; - } else { - // cell -> src - const int c = src->cell[i]; - const double* m = store_.mob (c); - const double* dm = store_.dmob(c); - - const double mt = m[0] + m[1]; - - assert (! ((m[0] < 0) || (m[1] < 0))); - assert (mt > 0); - - const double f = m[0] / mt; - const double df = ((1 - f)*dm[0] - f*dm[1]) / mt; - - *F += dt * dflux * f; - *J += dt * dflux * df; - } - } - template - void - initGravityTrans(const Grid& g , - const std::vector & htrans) { - - assert (htrans.size() == - static_cast::size_type>(g.cell_facepos[ g.number_of_cells ])); - - for (int f = 0; f < g.number_of_faces; ++f) { - store_.trans(f) = 0.0; - } - - for (int c = 0, i = 0; c < g.number_of_cells; ++c) { - for (; i < g.cell_facepos[c + 1]; ++i) { - int f = g.cell_faces[i]; - - assert (htrans[i] > 0.0); - - store_.trans(f) += 1.0 / htrans[i]; - } - } - - for (int f = 0; f < g.number_of_faces; ++f) { - store_.trans(f) = 1.0 / store_.trans(f); - } - - if (gravity_) { - this->computeStaticGravity(g); - } - } - - // ----------------------------------------------------------------- - // Newton control - // ----------------------------------------------------------------- - - template - void - initStep(const ReservoirState& state, - const Grid& g , - JacobianSystem& sys ) { - - (void) state; // Suppress 'unused' warning. - - typename JacobianSystem::vector_type& x = - sys.vector().writableSolution(); - - assert (x.size() == (::std::size_t) (g.number_of_cells)); - - if (init_step_use_previous_sol_) { - std::fill(x.begin(), x.end(), 0.0); - } else { - const std::vector& s = state.saturation(); - for (int c = 0, nc = g.number_of_cells; c < nc; ++c) { - // Impose s=0.5 at next time level as an NR initial value. - x[c] = 0.5 - s[2*c + 0]; - } - } - } - - template - bool - initIteration(const ReservoirState& state, - const Grid& g , - JacobianSystem& sys) { - - double s[2], mob[2], dmob[2 * 2], pc, dpc; - - const typename JacobianSystem::vector_type& x = - sys.vector().solution(); - const ::std::vector& sat = state.saturation(); - - bool in_range = true; - for (int c = 0; c < g.number_of_cells; ++c) { - store_.ds(c) = x[c]; // Store sat-change for accumulation(). - - s[0] = sat[c*2 + 0] + x[c]; - - double s_min = fluid_.s_min(c); - double s_max = fluid_.s_max(c); - - if ( s[0] < (s_min - sat_tol_) || s[0] > (s_max + sat_tol_) ) { - // if (s[0] < s_min){ - // std::cout << "Warning: s out of range, s-s_min = " << s_min-s[0] << std::endl; - // } - // if (s[0] > s_max){ - // std::cout << "Warning: s out of range, s-s_max = " << s[0]-s_max << std::endl; - // } - in_range = false; //line search fails - } - s[0] = std::max(s_min, s[0]); - s[0] = std::min(s_max, s[0]); - s[1] = 1 - s[0]; - - fluid_.mobility(c, s, mob, dmob); - fluid_.pc(c, s, pc, dpc); - - store_.mob (c)[0] = mob [0]; - store_.mob (c)[1] = mob [1]; - store_.dmob(c)[0] = dmob[0*2 + 0]; - store_.dmob(c)[1] = -dmob[1*2 + 1]; - store_.pc(c) = pc; - store_.dpc(c) = dpc; - } - if (!in_range) { -#ifdef VERBOSE - std::cout << "Warning: initIteration() - s was clamped in some cells.\n"; -#endif - } - return in_range; - } - - template - void - finishIteration(const ReservoirState& state, - const Grid& g , - NewtonIterate& it ) { - // Nothing to do at end of iteration in this model. - (void) state; (void) g; (void) it; - } - - template - void - finishStep(const Grid& g , - const SolutionVector& x , - ReservoirState& state) { - - double *s = &state.saturation()[0*2 + 0]; - - for (int c = 0; c < g.number_of_cells; ++c, s += 2) { - s[0] += x[c]; - double s_min = fluid_.s_min(c); - double s_max = fluid_.s_max(c); - -#if 0 - assert(s[0] >= s_min - sat_tol_); - assert(s[0] <= s_max + sat_tol_); -#endif - - s[0] = std::max(s_min, s[0]); - s[0] = std::min(s_max, s[0]); - s[1] = 1.0 - s[0]; - } - } - - private: - void - upwindMobility(const double dflux, - const double gflux, - const int* n , - int* pix , - double* m , - double* dm ) const { - bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) || - ( (! (dflux > 0)) && (! (gflux > 0)) ); - - if (equal_sign) { - - if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; } - else { pix[0] = 1; } - - m[0] = store_.mob(n[ pix[0] ]) [ 0 ]; - - if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; } - else { pix[1] = 1; } - - m[1] = store_.mob(n[ pix[1] ]) [ 1 ]; - - } else { - - if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; } - else { pix[1] = 1; } - - m[1] = store_.mob(n[ pix[1] ]) [ 1 ]; - - if (dflux + m[1]*gflux > 0) { pix[0] = 0; } - else { pix[0] = 1; } - - m[0] = store_.mob(n[ pix[0] ]) [ 0 ]; - } - - dm[0] = store_.dmob(n[ pix[0] ]) [ 0 ]; - dm[1] = store_.dmob(n[ pix[1] ]) [ 1 ]; - } - - template - void - computeStaticGravity(const Grid& g) { - - const int d = g.dimensions; - - for (int c = 0, i = 0; c < g.number_of_cells; ++c) { - const double* cc = g.cell_centroids + (c * d); - - for (; i < g.cell_facepos[c + 1]; ++i) { - const int f = g.cell_faces[i]; - const double* fc = g.face_centroids + (f * d); - - double dg = 0.0; - for (int j = 0; j < d; ++j) { - dg += gravity_[j] * (fc[j] - cc[j]); - } - - store_.dg(i) = store_.trans(f) * dg; - } - } - } - - double - gravityFlux(const int f) const { - double gflux; - - if (gravity_) { - int i1 = f2hf_[2*f + 0]; - int i2 = f2hf_[2*f + 1]; - - assert ((i1 >= 0) && (i2 >= 0)); - - gflux = store_.dg(i1) - store_.dg(i2); - gflux *= store_.drho(); - } else { - gflux = 0.0; - } - - return gflux; - } - void - capFlux(const int f,const int* n,double& pcflux, double* dpcflux) const { - //double capflux; - int i1 = n[0]; - int i2 = n[1]; - assert ((i1 >= 0) && (i2 >= 0)); - //double sgn=-1.0; - pcflux = store_.trans(f)*(store_.pc(i2) - store_.pc(i1)); - dpcflux[0] = -store_.trans(f)*store_.dpc(i1); - dpcflux[1] = store_.trans(f)*store_.dpc(i2); - } - - TwophaseFluid fluid_ ; - const double* gravity_; - std::vector f2hf_ ; - spu_2p::ModelParameterStorage store_ ; - bool init_step_use_previous_sol_; - double sat_tol_; - }; -} -#endif /* OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER */ diff --git a/opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp b/opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp deleted file mode 100644 index 84cb6943..00000000 --- a/opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp +++ /dev/null @@ -1,103 +0,0 @@ -/*=========================================================================== -// -// 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 . -*/ - -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif // HAVE_CONFIG_H - -#include -#include -#include - -#include - -namespace Opm -{ - - TransportSolverTwophaseImplicit::TransportSolverTwophaseImplicit( - const UnstructuredGrid& grid, - const Opm::IncompPropertiesInterface& props, - const std::vector& porevol, - const double* gravity, - const std::vector& half_trans, - const ParameterGroup& param) - : fluid_(props), - model_(fluid_, grid, porevol, gravity, param.getDefault("guess_old_solution", false)), - tsolver_(model_), - grid_(grid), - props_(props) - { - ctrl_.max_it = param.getDefault("max_it", 20); - ctrl_.verbosity = param.getDefault("verbosity", 0); - ctrl_.max_it_ls = param.getDefault("max_it_ls", 5); - model_.initGravityTrans(grid_, half_trans); - tsrc_ = create_transport_source(2, 2); - initial_porevolume_cell0_ = porevol[0]; - } - - TransportSolverTwophaseImplicit::~TransportSolverTwophaseImplicit() - { - destroy_transport_source(tsrc_); - } - - /// Solve for saturation at next timestep. - /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. For interpretation see Opm::computeTransportSource(). - /// \param[in] dt Time step. - /// \param[in, out] state Reservoir state. Calling solve() will read state.faceflux() and - /// read and write state.saturation(). - void TransportSolverTwophaseImplicit::solve(const double* porevolume, - const double* source, - const double dt, - TwophaseState& state) - { - // A very crude check for constant porosity (i.e. no rock-compressibility). - if (porevolume[0] != initial_porevolume_cell0_) { - OPM_THROW(std::runtime_error, "Detected changed pore volumes, but solver cannot handle rock compressibility."); - } - double ssrc[] = { 1.0, 0.0 }; - double dummy[] = { 0.0, 0.0 }; - clear_transport_source(tsrc_); - const int num_phases = 2; - for (int cell = 0; cell < grid_.number_of_cells; ++cell) { - int success = 1; - if (source[cell] > 0.0) { - success = append_transport_source(cell, num_phases, state.pressure()[cell], source[cell], ssrc, dummy, tsrc_); - } else if (source[cell] < 0.0) { - success = append_transport_source(cell, num_phases, state.pressure()[cell], source[cell], dummy, dummy, tsrc_); - } - if (!success) { - OPM_THROW(std::runtime_error, "Failed building TransportSource struct."); - } - } - Opm::ImplicitTransportDetails::NRReport rpt; - tsolver_.solve(grid_, tsrc_, dt, ctrl_, state, linsolver_, rpt); - std::cout << rpt; - } - -} // namespace Opm diff --git a/opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp b/opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp deleted file mode 100644 index f54ea12b..00000000 --- a/opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp +++ /dev/null @@ -1,119 +0,0 @@ -/* - 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 OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED -#define OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - - - -namespace Opm -{ - - // Implicit transport solver using Newton-Raphson. - class TransportSolverTwophaseImplicit : public TransportSolverTwophaseInterface - { - public: - /// Construct solver. - /// \param[in] grid A 2d or 3d grid. - /// \param[in] props Rock and fluid properties. - /// \param[in] porevol Pore volumes - /// \param[in] gravity Gravity vector (null for no gravity). - /// \param[in] half_trans Half-transmissibilities (one-sided) - /// \param[in] maxit Maximum number of non-linear iterations used. - TransportSolverTwophaseImplicit(const UnstructuredGrid& grid, - const Opm::IncompPropertiesInterface& props, - const std::vector& porevol, - const double* gravity, - const std::vector& half_trans, - const ParameterGroup& param); - - virtual ~TransportSolverTwophaseImplicit(); - - /// Solve for saturation at next timestep. - /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. For interpretation see Opm::computeTransportSource(). - /// \param[in] dt Time step. - /// \param[in, out] state Reservoir state. Calling solve() will read state.faceflux() and - /// read and write state.saturation(). - virtual void solve(const double* porevolume, - const double* source, - const double dt, - TwophaseState& state); - - private: - // Disallow copying and assignment. - TransportSolverTwophaseImplicit(const TransportSolverTwophaseImplicit&); - TransportSolverTwophaseImplicit& operator=(const TransportSolverTwophaseImplicit&); - - // Defining types for the underlying transport solver. - typedef Opm::SimpleFluid2pWrappingProps TwophaseFluid; - typedef Opm::SinglePointUpwindTwoPhase TransportModel; - 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; - - // Data members. - Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver linsolver_; - Opm::SimpleFluid2pWrappingProps fluid_; - SinglePointUpwindTwoPhase model_; - TransportSolver tsolver_; - const UnstructuredGrid& grid_; - Opm::ImplicitTransportDetails::NRControl ctrl_; - const Opm::IncompPropertiesInterface& props_; - TransportSource* tsrc_; - double initial_porevolume_cell0_; - }; - -} // namespace Opm - -#endif // OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED diff --git a/opm/core/transport/implicit/transport_source.c b/opm/core/transport/implicit/transport_source.c deleted file mode 100644 index 8daf198f..00000000 --- a/opm/core/transport/implicit/transport_source.c +++ /dev/null @@ -1,180 +0,0 @@ -/*=========================================================================== -// -// File: transport_source.c -// -// Created: 2011-10-05 19:58:38+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include - -#include - - -/* ---------------------------------------------------------------------- */ -static int -expand_source_tables(int alloc, struct TransportSource *src) -/* ---------------------------------------------------------------------- */ -{ - int *c; - double *p, *v, *s, *z; - - c = realloc(src->cell , alloc * 1 * sizeof *c); - p = realloc(src->pressure , alloc * 1 * sizeof *p); - v = realloc(src->flux , alloc * 1 * sizeof *v); - s = realloc(src->saturation, alloc * src->nphase * sizeof *s); - z = realloc(src->surfvolume, alloc * src->nphase * sizeof *z); - - if ((c == NULL) || - (p == NULL) || (v == NULL) || - (s == NULL) || (z == NULL)) { - - free(z); free(s); free(v); free(p); free(c); - - alloc = 0; - } else { - src->cell = c; src->cpty = alloc; - src->pressure = p; src->flux = v; - src->saturation = s; src->surfvolume = z; - } - - return alloc; -} - - -/* ====================================================================== - * Public methods below separator. - * ====================================================================== */ - - -/* ---------------------------------------------------------------------- */ -struct TransportSource * -create_transport_source(int nsrc, int nphase) -/* ---------------------------------------------------------------------- */ -{ - int status; - struct TransportSource *src; - - if (nphase <= 0) { - src = NULL; - } else { - src = malloc(1 * sizeof *src); - - if (src != NULL) { - src->nsrc = src->cpty = 0; - - src->nphase = nphase; - src->cell = NULL; - src->pressure = NULL; src->flux = NULL; - src->saturation = NULL; src->surfvolume = NULL; - - if (nsrc > 0) { - status = expand_source_tables(nsrc, src); - } else { - status = 1; - } - - if (status <= 0) { - destroy_transport_source(src); - src = NULL; - } - } - } - - return src; -} - - -/* ---------------------------------------------------------------------- */ -void -destroy_transport_source(struct TransportSource *src) -/* ---------------------------------------------------------------------- */ -{ - if (src != NULL) { - free(src->surfvolume); free(src->saturation); - free(src->flux) ; free(src->pressure) ; - free(src->cell) ; - } - - free(src); -} - - -/* ---------------------------------------------------------------------- */ -int -append_transport_source(int c, - int nphase, - double p, - double v, - const double *sat, - const double *z, - struct TransportSource *src) -/* ---------------------------------------------------------------------- */ -{ - int alloc, status; - - if (nphase != src->nphase) { - return -1; - } - - if (src->nsrc == src->cpty) { - alloc = (src->cpty > 0) ? 2 * src->cpty : 1; - status = expand_source_tables(alloc, src); - } else { - status = 1; - } - - if (status > 0) { - src->cell [ src->nsrc ] = c; - src->pressure[ src->nsrc ] = p; - src->flux [ src->nsrc ] = v; - - memcpy(src->saturation + (src->nsrc * nphase), sat, - nphase * sizeof *src->saturation); - - memcpy(src->surfvolume + (src->nsrc * nphase), z , - nphase * sizeof *src->surfvolume); - - src->nsrc += 1; - } - - return status > 0; -} - - -/* ---------------------------------------------------------------------- */ -void -clear_transport_source(struct TransportSource *src) -/* ---------------------------------------------------------------------- */ -{ - src->nsrc = 0; -} diff --git a/opm/core/transport/implicit/transport_source.h b/opm/core/transport/implicit/transport_source.h deleted file mode 100644 index b52de23f..00000000 --- a/opm/core/transport/implicit/transport_source.h +++ /dev/null @@ -1,79 +0,0 @@ -/*=========================================================================== -// -// File: transport_source.h -// -// Created: 2011-10-05 19:58:53+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_TRANSPORT_SOURCE_H_HEADER -#define OPM_TRANSPORT_SOURCE_H_HEADER - -#ifdef __cplusplus -extern "C" { -#endif - -struct TransportSource { - int nsrc; - int cpty; - - int nphase; - - int *cell; - double *pressure; - double *flux; - double *saturation; - double *surfvolume; -}; - - -struct TransportSource * -create_transport_source(int nsrc, int nphase); - -void -destroy_transport_source(struct TransportSource *src); - -int -append_transport_source(int c, - int nphase, - double p, - double v, - const double *sat, - const double *z, - struct TransportSource *src); - -void -clear_transport_source(struct TransportSource *src); - -#ifdef __cplusplus -} -#endif - -#endif /* OPM_TRANSPORT_SOURCE_H_HEADER */