OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

remove implicit transport code

Browse Source 
moved to opm-upscaling
This commit is contained in:
Arne Morten Kvarving 2018-01-09 16:09:28 +01:00
parent cc20b89e73
commit 529278d2b9
14 changed files with 0 additions and 2535 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
CMakeLists_files.cmake
View File

@ -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

198
opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp
View File

@ -1,198 +0,0 @@
/*===========================================================================
//
// File: CSRMatrixBlockAssembler.hpp
//
// Created: 2011-10-03 12:40:56+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER
#define OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER
#include <opm/core/transport/implicit/JacobianSystem.hpp>
#include <opm/core/linalg/sparse_sys.h>
#include <cassert>
#include <cstddef>
#include <cstdlib>
#include <algorithm>
#include <vector>
namespace Opm {
namespace ImplicitTransportDefault {
template <>
class MatrixZero <struct CSRMatrix> {
public:
static void
zero(struct CSRMatrix& A) {
csrmatrix_zero(&A);
}
};
template <>
class MatrixBlockAssembler<struct CSRMatrix> {
public:
template <class Block>
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 <class Connections>
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<int>(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 <class Connections>
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<int>(*c));
}
::std::sort(sconn_.begin(), sconn_.end());
esconn_.resize(0);
esconn_.reserve(ndof * sconn_.size());
for (::std::vector<int>::iterator
c = sconn_.begin(), e = sconn_.end(); c != e; ++c) {
for (::std::size_t dof = 0; dof < ndof; ++dof) {
esconn_.push_back(static_cast<int>((*c)*ndof + dof));
}
}
}
void
setCSRSize() {
mat_.m = ia_.size() - 1;
mat_.nnz = ja_.size() ;
}
::std::size_t ndof_;
::std::vector<int> ia_;
::std::vector<int> ja_;
::std::vector<double> sa_;
::std::vector<int> sconn_ ;
::std::vector<int> esconn_;
struct CSRMatrix mat_;
};
}
}
#endif /* OPM_CSRMATRIXBLOCKASSEMBLER_HPP_HEADER */

87
opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp
View File

@ -1,87 +0,0 @@
/*===========================================================================
//
// File: CSRMatrixUmfpackSolver.hpp
//
// Created: 2011-10-03 17:27:26+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_CSRMATRIXUMFPACKSOLVER_HPP_HEADER
#define OPM_CSRMATRIXUMFPACKSOLVER_HPP_HEADER
#include <opm/core/linalg/call_umfpack.h>
#include <opm/common/ErrorMacros.hpp>
namespace Opm
{
namespace ImplicitTransportLinAlgSupport
{
class CSRMatrixUmfpackSolver
{
public:
template <class Vector>
void
solve(const struct CSRMatrix* A,
const Vector b,
Vector x)
{
#if HAVE_SUITESPARSE_UMFPACK_H
call_UMFPACK(const_cast<CSRMatrix*>(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 <class Vector>
void
solve(const struct CSRMatrix& A,
const Vector& b,
Vector& x)
{
#if HAVE_SUITESPARSE_UMFPACK_H
call_UMFPACK(const_cast<CSRMatrix*>(&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 */

267
opm/core/transport/implicit/ImplicitAssembly.hpp
View File

@ -1,267 +0,0 @@
/*===========================================================================
//
// File: ImplicitAssembly.hpp
//
// Created: 2011-09-28 10:00:46+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_IMPLICITASSEMBLY_HPP_HEADER
#define OPM_IMPLICITASSEMBLY_HPP_HEADER
#include <cstddef>
#include <algorithm>
#include <vector>
namespace Opm {
template <class Model>
class ImplicitAssembly {
enum { DofPerCell = Model::DofPerCell };
public:
ImplicitAssembly(Model& model)
: model_(model),
nconn_(-1) ,
asm_buffer_()
{}
template <class Grid ,
class JacobianSystem>
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; pf<tsrc->pf; ++pf){
int nconn = 1;
::std::vector<int> 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 <class ReservoirState,
class Grid ,
class SourceTerms ,
class JacobianSystem>
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 <class Grid>
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 <class Grid>
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 <class Grid, class System>
void
createRowStructure(const Grid& g ,
const int c ,
System& sys) const {
int nconn = countConnections(g, c);
::std::vector<int> 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 <class ReservoirState, class Grid>
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 <class Grid, class System>
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 <class Grid, class SourceTerms, class System>
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<double> asm_buffer_;
};
}
#endif /* OPM_IMPLICITASSEMBLY_HPP_HEADER */

234
opm/core/transport/implicit/ImplicitTransport.hpp
View File

@ -1,234 +0,0 @@
/*===========================================================================
//
// File: ImplicitTransport.hpp
//
// Created: 2011-09-29 10:38:42+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_IMPLICITTRANSPORT_HPP_HEADER
#define OPM_IMPLICITTRANSPORT_HPP_HEADER
#include <opm/core/transport/implicit/ImplicitAssembly.hpp>
#include <iostream>
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 <class Ostream>
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 Model ,
class JacobianSystem ,
template <class> class VNorm ,
template <class> class VNeg ,
template <class> class VZero ,
template <class> class MZero ,
template <class> class VAsgn >
class ImplicitTransport {
public:
ImplicitTransport(Model& model)
: model_(model),
asm_ (model)
{}
template <class Grid ,
class SourceTerms ,
class ReservoirState,
class LinearSolver >
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<matrix_type>::zero(sys_.writableMatrix());
VZero<vector_type>::zero(sys_.vector().writableResidual());
asm_.assemble(state, g, src, dt, sys_);
const double nrm_res0 =
VNorm<vector_type>::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<vector_type>::zero(sys_.vector().writableIncrement());
linsolve.solve(sys_.matrix(),
sys_.vector().residual(),
sys_.vector().writableIncrement());
VNeg<vector_type>::negate(sys_.vector().writableIncrement());
//model_.finishIteration(state, g, sys_.vector());
rpt.norm_dx =
VNorm<vector_type>::norm(sys_.vector().increment());
// Begin line search
double residual=VNorm<vector_type>::norm(sys_.vector().residual());
int lin_it=0;
bool finished=rpt.norm_res<ctrl.atol;
double alpha=2.0;
// store old solution and increment before line search
vector_type dx_old(sys_.vector().increment());
vector_type x_old(sys_.vector().solution());
while(! finished){
alpha/=2.0;
VAsgn<vector_type>::assign(alpha, dx_old,
sys_.vector().writableIncrement());
VAsgn<vector_type>::assign(x_old,
sys_.vector().writableSolution());
sys_.vector().addIncrement();
init = model_.initIteration(state, g, sys_);
if (init) {
MZero<matrix_type>::zero(sys_.writableMatrix());
VZero<vector_type>::zero(sys_.vector().writableResidual());
asm_.assemble(state, g, src, dt, sys_);
residual = VNorm<vector_type>::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<vector_type>::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<Model> asm_;
JacobianSystem sys_;
};
}
#endif /* OPM_IMPLICITTRANSPORT_HPP_HEADER */

256
opm/core/transport/implicit/JacobianSystem.hpp
View File

@ -1,256 +0,0 @@
/*===========================================================================
//
// File: JacobianSystem.hpp
//
// Created: 2011-09-30 19:23:31+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_JACOBIANSYSTEM_HPP_HEADER
#define OPM_JACOBIANSYSTEM_HPP_HEADER
#include <cassert>
#include <cmath>
#include <cstddef>
#include <algorithm>
//#include <array>
#include <functional>
#include <numeric>
namespace Opm {
namespace ImplicitTransportDefault {
template <class BaseVec>
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<VT>());
}
};
template <class BaseVec>
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<VT>());
}
};
template <class BaseVec>
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 BaseVec>
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 <class Scalar>
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<VT>(), a));
}
};
template <class Matrix>
class MatrixZero;
template <class BaseVec>
class VectorBlockAssembler {
public:
template <class Block>
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 BaseVec>
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 BaseVec ,
template <class> class VSzSetter = VectorSizeSetter ,
template <class> class VAdd = VectorAdder ,
template <class> 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<BaseVec, 3>::iterator VAI;
for (VAI i = vcoll_.begin(), e = vcoll_.end(); i != e; ++i) {
VSzSetter<BaseVec>(*i).setSize(ndof, m);
}
#endif
for (::std::size_t i = 0; i < sizeof (vcoll_) / sizeof (vcoll_[0]); ++i) {
VSzSetter<BaseVec>(vcoll_[i]).setSize(ndof, m);
}
ndof_ = ndof;
}
void
addIncrement() {
VAdd<BaseVec>::add(vcoll_[ Increment ], vcoll_[ Solution ]);
}
template <class Block>
void
assembleBlock(::std::size_t ndof,
::std::size_t i ,
const Block& b ) {
assert (ndof_ > 0 );
assert (ndof == ndof_);
VBlkAsm<BaseVec>::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<BaseVec, 3> vcoll_;
};
template <class Matrix>
class MatrixBlockAssembler
/* {
public:
template <class Block>
void
assembleBlock(size_t n, size_t i, size j, const Block& b);
template <class Connections>
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 Matrix ,
class NVecCollection>
class JacobianSystem {
public:
JacobianSystem() {}
typedef Matrix matrix_type;
typedef MatrixBlockAssembler<Matrix> 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<Matrix> mba_ ; // Coefficient matrix
NVecCollection sysvec_; // Residual
};
}
}
#endif /* OPM_JACOBIANSYSTEM_HPP_HEADER */

100
opm/core/transport/implicit/NormSupport.hpp
View File

@ -1,100 +0,0 @@
/*===========================================================================
//
// File: NormSupport.hpp
//
// Created: 2011-10-04 19:37:35+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_NORMSUPPORT_HPP_HEADER
#define OPM_NORMSUPPORT_HPP_HEADER
#include <algorithm>
#include <cmath>
#include <functional>
#include <numeric>
namespace Opm {
namespace ImplicitTransportDefault {
template <typename T>
class MaxAbs : public ::std::binary_function <double, T, double> {
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 <typename T>
class SumAbs : public ::std::binary_function <double, T, double> {
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 <typename T>
class Euclid : public ::std::binary_function <double, T, double> {
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 Vector, template <typename> 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<VT>());
return NormImpl<VT>::postprocess(nrm);
}
};
}
}
#endif /* OPM_NORMSUPPORT_HPP_HEADER */

70
opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp
View File

@ -1,70 +0,0 @@
/*===========================================================================
//
// File: SimpleFluid2pWrappingProps.hpp
//
// Author: hnil <hnil@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#ifndef SIMPLEFLUID2PWRAPPINGPROPS_HPP
#define SIMPLEFLUID2PWRAPPINGPROPS_HPP
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <vector>
namespace Opm{
class SimpleFluid2pWrappingProps
{
public:
SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props);
double density(int phase) const;
template <class Sat,
class Mob,
class DMob>
void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const;
template <class Sat,
class Pcap,
class DPcap>
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<double> smin_;
std::vector<double> smax_;
};
}
#include <opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp>
#endif // SIMPLEFLUID2PWRAPPINGPROPS_HPP

104
opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp
View File

@ -1,104 +0,0 @@
/*===========================================================================
//
// File: SimpleFluid2pWrappingProps.cpp
//
// Author: hnil <hnil@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#ifndef SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP
#define SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP
#include <opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp>
#include <cassert>
#include <opm/common/ErrorMacros.hpp>
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<int> 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 <class Sat,
class Mob,
class DMob>
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 <class Sat,
class Pcap,
class DPcap>
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

725
opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp
View File

@ -1,725 +0,0 @@
/*===========================================================================
//
// File: SinglePointUpwindTwoPhase.hpp
//
// Created: 2011-09-28 14:21:34+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
/**
* \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 <cassert>
#include <cstddef>
#include <algorithm>
#include <vector>
#include <iostream>
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 <CODE>(f,c)</CODE> 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<double> data_;
};
}
template <class TwophaseFluid>
class SinglePointUpwindTwoPhase {
public:
template <class Grid>
SinglePointUpwindTwoPhase(const TwophaseFluid& fluid ,
const Grid& g ,
const std::vector<double>& 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<int>& p_faces ,
const std::vector<int>& hf_faces,
const std::vector<int>& nb_faces)
{
if (p_faces.empty()) {
return;
}
assert (p_faces.size() == hf_faces.size());
assert (hf_faces.size() == nb_faces.size());
std::vector<int> nbhf(hf_faces.size());
for (std::vector<int>::size_type i = 0; i < p_faces.size(); ++i) {
const int nbf = nb_faces[i];
assert (2*std::vector<int>::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<int>::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<int>::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 <class ReservoirState,
class Grid >
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 <class Grid>
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 <class Grid ,
class SourceTerms>
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 <class Grid>
void
initGravityTrans(const Grid& g ,
const std::vector<double> & htrans) {
assert (htrans.size() ==
static_cast<std::vector<double>::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 <class ReservoirState,
class Grid ,
class JacobianSystem>
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<double>& 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 <class ReservoirState,
class Grid ,
class JacobianSystem>
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<double>& 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 <class ReservoirState,
class Grid ,
class NewtonIterate >
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 <class Grid ,
class SolutionVector,
class ReservoirState>
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 <class Grid>
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<int> f2hf_ ;
spu_2p::ModelParameterStorage store_ ;
bool init_step_use_previous_sol_;
double sat_tol_;
};
}
#endif /* OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER */

103
opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp
View File

@ -1,103 +0,0 @@
/*===========================================================================
//
// File: ImpliciteTwoPhaseTransportSolver.cpp
//
// Author: hnil <hnil@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <iostream>
namespace Opm
{
TransportSolverTwophaseImplicit::TransportSolverTwophaseImplicit(
const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const std::vector<double>& porevol,
const double* gravity,
const std::vector<double>& 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

119
opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED
#define OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED
#include <opm/core/transport/TransportSolverTwophaseInterface.hpp>
#include <opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp>
#include <opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/transport/implicit/ImplicitTransport.hpp>
#include <opm/core/transport/implicit/transport_source.h>
#include <opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/implicit/NormSupport.hpp>
#include <opm/core/transport/implicit/ImplicitAssembly.hpp>
#include <opm/core/transport/implicit/ImplicitTransport.hpp>
#include <opm/core/transport/implicit/JacobianSystem.hpp>
#include <opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <memory>
#include <vector>
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<double>& porevol,
const double* gravity,
const std::vector<double>& 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<TwophaseFluid> TransportModel;
typedef ImplicitTransportDefault::NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef ImplicitTransportDefault::JacobianSystem < struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public:
static double
norm(const Vector& v) {
return ImplicitTransportDefault::AccumulationNorm <Vector, ImplicitTransportDefault::MaxAbs>::norm(v);
}
};
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
ImplicitTransportDefault::VectorNegater ,
ImplicitTransportDefault::VectorZero ,
ImplicitTransportDefault::MatrixZero ,
ImplicitTransportDefault::VectorAssign > TransportSolver;
// Data members.
Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver linsolver_;
Opm::SimpleFluid2pWrappingProps fluid_;
SinglePointUpwindTwoPhase<Opm::SimpleFluid2pWrappingProps> 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

180
opm/core/transport/implicit/transport_source.c
View File

@ -1,180 +0,0 @@
/*===========================================================================
//
// File: transport_source.c
//
// Created: 2011-10-05 19:58:38+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <stdlib.h>
#include <string.h>
#include <opm/core/transport/implicit/transport_source.h>
/* ---------------------------------------------------------------------- */
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;
}

79
opm/core/transport/implicit/transport_source.h
View File

@ -1,79 +0,0 @@
/*===========================================================================
//
// File: transport_source.h
//
// Created: 2011-10-05 19:58:53+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#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 */