Further reorganising of opm-core.
Deleted some unused code (or moved to opm-porsol), moved all code dealing with time-of-flight to opm/core/tof, moved code for implicit transport solver to opm/core/transport/implicit, spu_[im|ex]plicit.[ch] to opm/core/transport/minimal.
This commit is contained in:
Atgeirr Flø Rasmussen
parent
2405758e2d
commit
987aa5b6fd
@ -1,328 +0,0 @@
|
|||||||
/*
|
|
||||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
|
||||||
|
|
||||||
OPM is free software: you can redistribute it and/or modify
|
|
||||||
it under the terms of the GNU General Public License as published by
|
|
||||||
the Free Software Foundation, either version 3 of the License, or
|
|
||||||
(at your option) any later version.
|
|
||||||
|
|
||||||
OPM is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
You should have received a copy of the GNU General Public License
|
|
||||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
|
|
||||||
#ifndef OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
#define OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
|
|
||||||
#include <opm/core/pressure/fsh.h>
|
|
||||||
#include <opm/core/linalg/sparse_sys.h>
|
|
||||||
#include <opm/core/pressure/mimetic/mimetic.h>
|
|
||||||
#include <opm/core/GridAdapter.hpp>
|
|
||||||
#include <stdexcept>
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Encapsulates the ifsh (= incompressible flow solver hybrid) solver modules.
|
|
||||||
class HybridPressureSolver
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
/// @brief
|
|
||||||
/// Default constructor, does nothing.
|
|
||||||
HybridPressureSolver()
|
|
||||||
: state_(Uninitialized), data_(0)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Destructor.
|
|
||||||
~HybridPressureSolver()
|
|
||||||
{
|
|
||||||
fsh_destroy(data_);
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Initialize the solver's structures for a given grid (at some point also well pattern).
|
|
||||||
/// @tparam Grid This must conform to the SimpleGrid concept.
|
|
||||||
/// @param grid The grid object.
|
|
||||||
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
|
|
||||||
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
|
|
||||||
template <class Grid>
|
|
||||||
void init(const Grid& grid, const double* perm, const double* gravity)
|
|
||||||
{
|
|
||||||
// Build C grid structure.
|
|
||||||
grid_.init(grid);
|
|
||||||
|
|
||||||
// Checking if grids are properly initialized. Move this code to a test somewhere.
|
|
||||||
// GridAdapter grid2;
|
|
||||||
// grid2.init(grid_);
|
|
||||||
// if (grid2 == grid_) {
|
|
||||||
// std::cout << "Grids are equal." << std::endl;
|
|
||||||
// } else {
|
|
||||||
// std::cout << "Grids are NOT equal." << std::endl;
|
|
||||||
// }
|
|
||||||
|
|
||||||
// Build (empty for now) C well structure.
|
|
||||||
well_t* w = 0;
|
|
||||||
|
|
||||||
// Initialize ifsh data.
|
|
||||||
data_ = ifsh_construct(grid_.c_grid(), w);
|
|
||||||
if (!data_) {
|
|
||||||
throw std::runtime_error("Failed to initialize ifsh solver.");
|
|
||||||
}
|
|
||||||
|
|
||||||
// Compute inner products, gravity contributions.
|
|
||||||
int num_cells = grid.numCells();
|
|
||||||
int ngconn = grid_.c_grid()->cell_facepos[num_cells];
|
|
||||||
gpress_.clear();
|
|
||||||
gpress_.resize(ngconn, 0.0);
|
|
||||||
int ngconn2 = data_->sum_ngconn2;
|
|
||||||
Binv_.resize(ngconn2);
|
|
||||||
ncf_.resize(num_cells);
|
|
||||||
typename Grid::Vector grav;
|
|
||||||
std::copy(gravity, gravity + Grid::dimension, &grav[0]);
|
|
||||||
int count = 0;
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
int num_local_faces = grid.numCellFaces(cell);
|
|
||||||
ncf_[cell] = num_local_faces;
|
|
||||||
typename Grid::Vector cc = grid.cellCentroid(cell);
|
|
||||||
for (int local_ix = 0; local_ix < num_local_faces; ++local_ix) {
|
|
||||||
int face = grid.cellFace(cell, local_ix);
|
|
||||||
typename Grid::Vector fc = grid.faceCentroid(face);
|
|
||||||
gpress_[count++] = grav*(fc - cc);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
assert(count == ngconn);
|
|
||||||
|
|
||||||
UnstructuredGrid* g = grid_.c_grid();
|
|
||||||
mim_ip_simple_all(g->number_of_cells, g->dimensions,
|
|
||||||
data_->max_ngconn,
|
|
||||||
g->cell_facepos, g->cell_faces,
|
|
||||||
g->face_cells, g->face_centroids,
|
|
||||||
g->face_normals, g->face_areas,
|
|
||||||
g->cell_centroids, g->cell_volumes,
|
|
||||||
const_cast<double*>(perm), &Binv_[0]);
|
|
||||||
state_ = Initialized;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
enum FlowBCTypes { FBC_UNSET, FBC_PRESSURE, FBC_FLUX };
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Assemble the sparse system.
|
|
||||||
/// You must call init() prior to calling assemble().
|
|
||||||
/// @param sources Source terms, one per cell. Positive numbers
|
|
||||||
/// are sources, negative are sinks.
|
|
||||||
/// @param total_mobilities Scalar total mobilities, one per cell.
|
|
||||||
/// @param omegas Gravity term, one per cell. In a multi-phase
|
|
||||||
/// flow setting this is equal to
|
|
||||||
/// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f]
|
|
||||||
/// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a
|
|
||||||
/// phase density and \f$\lambda_t\f$ is the total mobility.
|
|
||||||
void assemble(const std::vector<double>& sources,
|
|
||||||
const std::vector<double>& total_mobilities,
|
|
||||||
const std::vector<double>& omegas,
|
|
||||||
const std::vector<FlowBCTypes>& bctypes,
|
|
||||||
const std::vector<double>& bcvalues)
|
|
||||||
{
|
|
||||||
if (state_ == Uninitialized) {
|
|
||||||
throw std::runtime_error("Error in HybridPressureSolver::assemble(): You must call init() prior to calling assemble().");
|
|
||||||
}
|
|
||||||
|
|
||||||
// Boundary conditions.
|
|
||||||
|
|
||||||
assert (bctypes.size() ==
|
|
||||||
static_cast<std::vector<FlowBCTypes>::size_type>(grid_.numFaces()));
|
|
||||||
|
|
||||||
FlowBoundaryConditions *bc = gather_boundary_conditions(bctypes, bcvalues);
|
|
||||||
|
|
||||||
// Source terms from user.
|
|
||||||
double* src = const_cast<double*>(&sources[0]); // Ugly? Yes. Safe? I think so.
|
|
||||||
|
|
||||||
// All well related things are zero.
|
|
||||||
well_control_t* wctrl = 0;
|
|
||||||
double* WI = 0;
|
|
||||||
double* wdp = 0;
|
|
||||||
|
|
||||||
// Scale inner products and gravity terms by saturation-dependent factors.
|
|
||||||
UnstructuredGrid* g = grid_.c_grid();
|
|
||||||
Binv_mobilityweighted_.resize(Binv_.size());
|
|
||||||
mim_ip_mobility_update(g->number_of_cells, g->cell_facepos, &total_mobilities[0],
|
|
||||||
&Binv_[0], &Binv_mobilityweighted_[0]);
|
|
||||||
gpress_omegaweighted_.resize(gpress_.size());
|
|
||||||
mim_ip_density_update(g->number_of_cells, g->cell_facepos, &omegas[0],
|
|
||||||
&gpress_[0], &gpress_omegaweighted_[0]);
|
|
||||||
|
|
||||||
|
|
||||||
// Zero the linalg structures.
|
|
||||||
csrmatrix_zero(data_->A);
|
|
||||||
for (std::size_t i = 0; i < data_->A->m; i++) {
|
|
||||||
data_->b[i] = 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
// Assemble the embedded linear system.
|
|
||||||
ifsh_assemble(bc, src, &Binv_mobilityweighted_[0], &gpress_omegaweighted_[0],
|
|
||||||
wctrl, WI, wdp, data_);
|
|
||||||
state_ = Assembled;
|
|
||||||
|
|
||||||
flow_conditions_destroy(bc);
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Encapsulate a sparse linear system in CSR format.
|
|
||||||
struct LinearSystem
|
|
||||||
{
|
|
||||||
int n;
|
|
||||||
int nnz;
|
|
||||||
int* ia;
|
|
||||||
int* ja;
|
|
||||||
double* sa;
|
|
||||||
double* b;
|
|
||||||
double* x;
|
|
||||||
};
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Access the linear system assembled.
|
|
||||||
/// You must call assemble() prior to calling linearSystem().
|
|
||||||
/// @param[out] s The linear system encapsulation to modify.
|
|
||||||
/// After this call, s will point to linear system structures
|
|
||||||
/// that are owned and allocated internally.
|
|
||||||
void linearSystem(LinearSystem& s)
|
|
||||||
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in HybridPressureSolver::linearSystem(): "
|
|
||||||
"You must call assemble() prior to calling linearSystem().");
|
|
||||||
}
|
|
||||||
s.n = data_->A->m;
|
|
||||||
s.nnz = data_->A->nnz;
|
|
||||||
s.ia = data_->A->ia;
|
|
||||||
s.ja = data_->A->ja;
|
|
||||||
s.sa = data_->A->sa;
|
|
||||||
s.b = data_->b;
|
|
||||||
s.x = data_->x;
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Compute cell pressures and face fluxes.
|
|
||||||
/// You must call assemble() (and solve the linear system accessed
|
|
||||||
/// by calling linearSystem()) prior to calling
|
|
||||||
/// computePressuresAndFluxes().
|
|
||||||
/// @param[out] cell_pressures Cell pressure values.
|
|
||||||
/// @param[out] face_areas Face flux values.
|
|
||||||
void computePressuresAndFluxes(std::vector<double>& cell_pressures,
|
|
||||||
std::vector<double>& face_fluxes)
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in HybridPressureSolver::computePressuresAndFluxes(): "
|
|
||||||
"You must call assemble() (and solve the linear system) "
|
|
||||||
"prior to calling computePressuresAndFluxes().");
|
|
||||||
}
|
|
||||||
int num_cells = grid_.c_grid()->number_of_cells;
|
|
||||||
int num_faces = grid_.c_grid()->number_of_faces;
|
|
||||||
cell_pressures.clear();
|
|
||||||
cell_pressures.resize(num_cells, 0.0);
|
|
||||||
face_fluxes.clear();
|
|
||||||
face_fluxes.resize(num_faces, 0.0);
|
|
||||||
fsh_press_flux(grid_.c_grid(), &Binv_mobilityweighted_[0], &gpress_omegaweighted_[0],
|
|
||||||
data_, &cell_pressures[0], &face_fluxes[0], 0, 0);
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Compute cell fluxes from face fluxes.
|
|
||||||
/// You must call assemble() (and solve the linear system accessed
|
|
||||||
/// by calling linearSystem()) prior to calling
|
|
||||||
/// faceFluxToCellFlux().
|
|
||||||
/// @param face_fluxes
|
|
||||||
/// @param face_areas Face flux values (usually output from computePressuresAndFluxes()).
|
|
||||||
/// @param[out] cell_fluxes Cell-wise flux values.
|
|
||||||
/// They are given in cell order, and for each cell there is
|
|
||||||
/// one value for each adjacent face (in the same order as the
|
|
||||||
/// cell-face topology of the grid). Positive values represent
|
|
||||||
/// fluxes out of the cell.
|
|
||||||
void faceFluxToCellFlux(const std::vector<double>& face_fluxes,
|
|
||||||
std::vector<double>& cell_fluxes)
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in HybridPressureSolver::faceFluxToCellFlux(): "
|
|
||||||
"You must call assemble() (and solve the linear system) "
|
|
||||||
"prior to calling faceFluxToCellFlux().");
|
|
||||||
}
|
|
||||||
const UnstructuredGrid& g = *(grid_.c_grid());
|
|
||||||
int num_cells = g.number_of_cells;
|
|
||||||
cell_fluxes.resize(g.cell_facepos[num_cells]);
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
|
|
||||||
int face = g.cell_faces[hface];
|
|
||||||
bool pos = (g.face_cells[2*face] == cell);
|
|
||||||
cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Access the number of connections (faces) per cell. Deprecated, will be removed.
|
|
||||||
const std::vector<int>& numCellFaces()
|
|
||||||
{
|
|
||||||
return ncf_;
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
// Disabling copy and assigment for now.
|
|
||||||
HybridPressureSolver(const HybridPressureSolver&);
|
|
||||||
HybridPressureSolver& operator=(const HybridPressureSolver&);
|
|
||||||
|
|
||||||
enum State { Uninitialized, Initialized, Assembled };
|
|
||||||
State state_;
|
|
||||||
|
|
||||||
// Solver data.
|
|
||||||
fsh_data* data_;
|
|
||||||
// Grid.
|
|
||||||
GridAdapter grid_;
|
|
||||||
// Number of faces per cell.
|
|
||||||
std::vector<int> ncf_;
|
|
||||||
// B^{-1} storage.
|
|
||||||
std::vector<double> Binv_;
|
|
||||||
std::vector<double> Binv_mobilityweighted_;
|
|
||||||
// Gravity contributions.
|
|
||||||
std::vector<double> gpress_;
|
|
||||||
std::vector<double> gpress_omegaweighted_;
|
|
||||||
|
|
||||||
|
|
||||||
FlowBoundaryConditions*
|
|
||||||
gather_boundary_conditions(const std::vector<FlowBCTypes>& bctypes ,
|
|
||||||
const std::vector<double>& bcvalues)
|
|
||||||
{
|
|
||||||
FlowBoundaryConditions* fbc = flow_conditions_construct(0);
|
|
||||||
|
|
||||||
int ok = fbc != 0;
|
|
||||||
std::vector<FlowBCTypes>::size_type i;
|
|
||||||
|
|
||||||
for (i = 0; ok && (i < bctypes.size()); ++i) {
|
|
||||||
if (bctypes[ i ] == FBC_PRESSURE) {
|
|
||||||
ok = flow_conditions_append(BC_PRESSURE,
|
|
||||||
static_cast<int>(i),
|
|
||||||
bcvalues[ i ],
|
|
||||||
fbc);
|
|
||||||
}
|
|
||||||
else if (bctypes[ i ] == FBC_FLUX) {
|
|
||||||
ok = flow_conditions_append(BC_FLUX_TOTVOL,
|
|
||||||
static_cast<int>(i),
|
|
||||||
bcvalues[ i ],
|
|
||||||
fbc);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return fbc;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
}; // class HybridPressureSolver
|
|
||||||
|
|
||||||
|
|
||||||
#endif // OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
@ -1,527 +0,0 @@
|
|||||||
/*
|
|
||||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
|
||||||
|
|
||||||
OPM is free software: you can redistribute it and/or modify
|
|
||||||
it under the terms of the GNU General Public License as published by
|
|
||||||
the Free Software Foundation, either version 3 of the License, or
|
|
||||||
(at your option) any later version.
|
|
||||||
|
|
||||||
OPM is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
You should have received a copy of the GNU General Public License
|
|
||||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
#define OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
|
|
||||||
|
|
||||||
#include <opm/core/pressure/tpfa/cfs_tpfa.h>
|
|
||||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
|
||||||
#include <opm/core/linalg/sparse_sys.h>
|
|
||||||
#include <opm/core/pressure/flow_bc.h>
|
|
||||||
#include <opm/core/pressure/legacy_well.h>
|
|
||||||
#include <opm/core/pressure/tpfa/compr_quant.h>
|
|
||||||
#include <opm/core/GridAdapter.hpp>
|
|
||||||
#include <stdexcept>
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Encapsulates the cfs_tpfa (= compressible flow solver
|
|
||||||
/// two-point flux approximation) solver modules.
|
|
||||||
class TPFACompressiblePressureSolver
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
/// @brief
|
|
||||||
/// Default constructor, does nothing.
|
|
||||||
TPFACompressiblePressureSolver()
|
|
||||||
: state_(Uninitialized), data_(0), bc_(0)
|
|
||||||
{
|
|
||||||
wells_.number_of_wells = 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Destructor.
|
|
||||||
~TPFACompressiblePressureSolver()
|
|
||||||
{
|
|
||||||
flow_conditions_destroy(bc_);
|
|
||||||
cfs_tpfa_destroy(data_);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Initialize the solver's structures for a given grid, for well setup also call initWells().
|
|
||||||
/// @tparam Grid This must conform to the SimpleGrid concept.
|
|
||||||
/// @tparam Wells This must conform to the SimpleWells concept.
|
|
||||||
/// @param grid The grid object.
|
|
||||||
/// @param wells Well specifications.
|
|
||||||
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
|
|
||||||
/// @param perm Porosity by cell.
|
|
||||||
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
|
|
||||||
template <class Grid, class Wells>
|
|
||||||
void init(const Grid& grid, const Wells& wells, const double* perm, const double* porosity,
|
|
||||||
const typename Grid::Vector& gravity)
|
|
||||||
{
|
|
||||||
initWells(wells);
|
|
||||||
init(grid, perm, porosity, gravity);
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Initialize the solver's structures for a given grid, for well setup also call initWells().
|
|
||||||
/// @tparam Grid This must conform to the SimpleGrid concept.
|
|
||||||
/// @param grid The grid object.
|
|
||||||
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
|
|
||||||
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
|
|
||||||
template <class Grid>
|
|
||||||
void init(const Grid& grid, const double* perm, const double* porosity, const typename Grid::Vector& gravity)
|
|
||||||
{
|
|
||||||
// Build C grid structure.
|
|
||||||
grid_.init(grid);
|
|
||||||
|
|
||||||
// Initialize data.
|
|
||||||
int num_phases = 3;
|
|
||||||
well_t* w = 0;
|
|
||||||
if (wells_.number_of_wells != 0) {
|
|
||||||
w = &wells_;
|
|
||||||
}
|
|
||||||
data_ = cfs_tpfa_construct(grid_.c_grid(), w, num_phases);
|
|
||||||
if (!data_) {
|
|
||||||
throw std::runtime_error("Failed to initialize cfs_tpfa solver.");
|
|
||||||
}
|
|
||||||
|
|
||||||
// Compute half-transmissibilities
|
|
||||||
int num_cells = grid.numCells();
|
|
||||||
int ngconn = grid_.c_grid()->cell_facepos[num_cells];
|
|
||||||
ncf_.resize(num_cells);
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
int num_local_faces = grid.numCellFaces(cell);
|
|
||||||
ncf_[cell] = num_local_faces;
|
|
||||||
}
|
|
||||||
htrans_.resize(ngconn);
|
|
||||||
tpfa_htrans_compute(grid_.c_grid(), perm, &htrans_[0]);
|
|
||||||
|
|
||||||
// Compute transmissibilities.
|
|
||||||
trans_.resize(grid_.numFaces());
|
|
||||||
tpfa_trans_compute(grid_.c_grid(), &htrans_[0], &trans_[0]);
|
|
||||||
|
|
||||||
// Compute pore volumes.
|
|
||||||
porevol_.resize(num_cells);
|
|
||||||
for (int i = 0; i < num_cells; ++i) {
|
|
||||||
porevol_[i] = porosity[i]*grid.cellVolume(i);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Set gravity.
|
|
||||||
if (Grid::dimension != 3) {
|
|
||||||
throw std::logic_error("Only 3 dimensions supported currently.");
|
|
||||||
}
|
|
||||||
std::copy(gravity.begin(), gravity.end(), gravity_);
|
|
||||||
|
|
||||||
state_ = Initialized;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// Boundary condition types.
|
|
||||||
enum FlowBCTypes { FBC_UNSET = BC_NOFLOW ,
|
|
||||||
FBC_PRESSURE = BC_PRESSURE ,
|
|
||||||
FBC_FLUX = BC_FLUX_TOTVOL };
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Assemble the sparse system.
|
|
||||||
/// You must call init() prior to calling assemble().
|
|
||||||
/// @param sources Source terms, one per cell. Positive numbers
|
|
||||||
/// are sources, negative are sinks.
|
|
||||||
/// @param total_mobilities Scalar total mobilities, one per cell.
|
|
||||||
/// @param omegas Gravity term, one per cell. In a multi-phase
|
|
||||||
/// flow setting this is equal to
|
|
||||||
/// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f]
|
|
||||||
/// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a
|
|
||||||
/// phase density and \f$\lambda_t\f$ is the total mobility.
|
|
||||||
void assemble(const double* sources,
|
|
||||||
const FlowBCTypes* bctypes,
|
|
||||||
const double* bcvalues,
|
|
||||||
const double dt,
|
|
||||||
const double* totcompr,
|
|
||||||
const double* voldiscr,
|
|
||||||
const double* cellA, // num phases^2 * num cells, fortran ordering!
|
|
||||||
const double* faceA, // num phases^2 * num faces, fortran ordering!
|
|
||||||
const double* wellperfA,
|
|
||||||
const double* phasemobf,
|
|
||||||
const double* phasemobwellperf,
|
|
||||||
const double* cell_pressure,
|
|
||||||
const double* gravcapf,
|
|
||||||
const double* wellperf_gpot,
|
|
||||||
const double* surf_dens)
|
|
||||||
{
|
|
||||||
if (state_ == Uninitialized) {
|
|
||||||
throw std::runtime_error("Error in TPFACompressiblePressureSolver::assemble(): You must call init() prior to calling assemble().");
|
|
||||||
}
|
|
||||||
UnstructuredGrid* g = grid_.c_grid();
|
|
||||||
|
|
||||||
// Boundary conditions.
|
|
||||||
gather_boundary_conditions(bctypes, bcvalues);
|
|
||||||
|
|
||||||
// Source terms from user.
|
|
||||||
double* src = const_cast<double*>(sources); // Ugly? Yes. Safe? I think so.
|
|
||||||
|
|
||||||
// Wells.
|
|
||||||
well_t* wells = NULL;
|
|
||||||
well_control_t* wctrl = NULL;
|
|
||||||
struct completion_data* wcompl = NULL;
|
|
||||||
if (wells_.number_of_wells != 0) {
|
|
||||||
wells = &wells_;
|
|
||||||
wctrl = &wctrl_;
|
|
||||||
wcompl = &wcompl_;
|
|
||||||
// The next objects already have the correct sizes.
|
|
||||||
std::copy(wellperf_gpot, wellperf_gpot + well_gpot_storage_.size(), well_gpot_storage_.begin());
|
|
||||||
std::copy(wellperfA, wellperfA + well_A_storage_.size(), well_A_storage_.begin());
|
|
||||||
std::copy(phasemobwellperf, phasemobwellperf + well_phasemob_storage_.size(), well_phasemob_storage_.begin());
|
|
||||||
}
|
|
||||||
|
|
||||||
// Assemble the embedded linear system.
|
|
||||||
compr_quantities cq = { 3 ,
|
|
||||||
const_cast<double *>(totcompr ) ,
|
|
||||||
const_cast<double *>(voldiscr ) ,
|
|
||||||
const_cast<double *>(cellA ) ,
|
|
||||||
const_cast<double *>(faceA ) ,
|
|
||||||
const_cast<double *>(phasemobf) };
|
|
||||||
|
|
||||||
// Call the assembly routine. After this, linearSystem() may be called.
|
|
||||||
cfs_tpfa_assemble(g, dt, wells, bc_, src,
|
|
||||||
&cq, &trans_[0], gravcapf,
|
|
||||||
wctrl, wcompl,
|
|
||||||
cell_pressure, &porevol_[0],
|
|
||||||
data_);
|
|
||||||
phasemobf_.assign(phasemobf, phasemobf + grid_.numFaces()*3);
|
|
||||||
gravcapf_.assign(gravcapf, gravcapf + grid_.numFaces()*3);
|
|
||||||
state_ = Assembled;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// Encapsulate a sparse linear system in CSR format.
|
|
||||||
struct LinearSystem
|
|
||||||
{
|
|
||||||
int n;
|
|
||||||
int nnz;
|
|
||||||
int* ia;
|
|
||||||
int* ja;
|
|
||||||
double* sa;
|
|
||||||
double* b;
|
|
||||||
double* x;
|
|
||||||
};
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Access the linear system assembled.
|
|
||||||
/// You must call assemble() prior to calling linearSystem().
|
|
||||||
/// @param[out] s The linear system encapsulation to modify.
|
|
||||||
/// After this call, s will point to linear system structures
|
|
||||||
/// that are owned and allocated internally.
|
|
||||||
void linearSystem(LinearSystem& s)
|
|
||||||
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in TPFACompressiblePressureSolver::linearSystem(): "
|
|
||||||
"You must call assemble() prior to calling linearSystem().");
|
|
||||||
}
|
|
||||||
s.n = data_->A->m;
|
|
||||||
s.nnz = data_->A->nnz;
|
|
||||||
s.ia = data_->A->ia;
|
|
||||||
s.ja = data_->A->ja;
|
|
||||||
s.sa = data_->A->sa;
|
|
||||||
s.b = data_->b;
|
|
||||||
s.x = data_->x;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Compute cell pressures and face fluxes.
|
|
||||||
/// You must call assemble() (and solve the linear system accessed
|
|
||||||
/// by calling linearSystem()) prior to calling
|
|
||||||
/// computePressuresAndFluxes().
|
|
||||||
/// @param[out] cell_pressures Cell pressure values.
|
|
||||||
/// @param[out] face_areas Face flux values.
|
|
||||||
void computePressuresAndFluxes(std::vector<double>& cell_pressures,
|
|
||||||
std::vector<double>& face_pressures,
|
|
||||||
std::vector<double>& face_fluxes,
|
|
||||||
std::vector<double>& well_pressures,
|
|
||||||
std::vector<double>& well_fluxes)
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in TPFACompressiblePressureSolver::computePressuresAndFluxes(): "
|
|
||||||
"You must call assemble() (and solve the linear system) "
|
|
||||||
"prior to calling computePressuresAndFluxes().");
|
|
||||||
}
|
|
||||||
int num_cells = grid_.c_grid()->number_of_cells;
|
|
||||||
int num_faces = grid_.c_grid()->number_of_faces;
|
|
||||||
cell_pressures.clear();
|
|
||||||
cell_pressures.resize(num_cells, 0.0);
|
|
||||||
face_pressures.clear();
|
|
||||||
face_pressures.resize(num_faces, 0.0);
|
|
||||||
face_fluxes.clear();
|
|
||||||
face_fluxes.resize(num_faces, 0.0);
|
|
||||||
|
|
||||||
int np = 3; // Number of phases.
|
|
||||||
|
|
||||||
// Wells.
|
|
||||||
well_t* wells = NULL;
|
|
||||||
struct completion_data* wcompl = NULL;
|
|
||||||
double* wpress = 0;
|
|
||||||
double* wflux = 0;
|
|
||||||
if (wells_.number_of_wells != 0) {
|
|
||||||
wells = &wells_;
|
|
||||||
wcompl = &wcompl_;
|
|
||||||
well_pressures.resize(wells_.number_of_wells);
|
|
||||||
well_fluxes.resize(well_cells_storage_.size());
|
|
||||||
wpress = &well_pressures[0];
|
|
||||||
wflux = &well_fluxes[0];
|
|
||||||
}
|
|
||||||
|
|
||||||
cfs_tpfa_press_flux(grid_.c_grid(),
|
|
||||||
bc_, wells,
|
|
||||||
np, &trans_[0], &phasemobf_[0], &gravcapf_[0],
|
|
||||||
wcompl,
|
|
||||||
data_, &cell_pressures[0], &face_fluxes[0],
|
|
||||||
wpress, wflux);
|
|
||||||
cfs_tpfa_fpress(grid_.c_grid(), bc_, np, &htrans_[0],
|
|
||||||
&phasemobf_[0], &gravcapf_[0], data_, &cell_pressures[0],
|
|
||||||
&face_fluxes[0], &face_pressures[0]);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Explicit IMPES time step limit.
|
|
||||||
double explicitTimestepLimit(const double* faceA, // num phases^2 * num faces, fortran ordering!
|
|
||||||
const double* phasemobf,
|
|
||||||
const double* phasemobf_deriv,
|
|
||||||
const double* surf_dens)
|
|
||||||
{
|
|
||||||
compr_quantities cq = { 3, // nphases
|
|
||||||
0, // totcompr
|
|
||||||
0, // voldiscr
|
|
||||||
0, // Ac
|
|
||||||
const_cast<double *>(faceA) ,
|
|
||||||
const_cast<double *>(phasemobf) };
|
|
||||||
return cfs_tpfa_impes_maxtime(grid_.c_grid(), &cq, &trans_[0], &porevol_[0], data_,
|
|
||||||
phasemobf_deriv, surf_dens, gravity_);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Explicit IMPES transport.
|
|
||||||
void explicitTransport(const double dt,
|
|
||||||
double* cell_surfvols)
|
|
||||||
{
|
|
||||||
int np = 3; // Number of phases.
|
|
||||||
|
|
||||||
well_t* wells = NULL;
|
|
||||||
if (wells_.number_of_wells != 0) {
|
|
||||||
wells = &wells_;
|
|
||||||
}
|
|
||||||
cfs_tpfa_expl_mass_transport(grid_.c_grid(), wells, &wcompl_,
|
|
||||||
np, dt, &porevol_[0],
|
|
||||||
data_, cell_surfvols);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Compute cell fluxes from face fluxes.
|
|
||||||
/// You must call assemble() (and solve the linear system accessed
|
|
||||||
/// by calling linearSystem()) prior to calling
|
|
||||||
/// faceFluxToCellFlux().
|
|
||||||
/// @param face_fluxes
|
|
||||||
/// @param face_areas Face flux values (usually output from computePressuresAndFluxes()).
|
|
||||||
/// @param[out] cell_fluxes Cell-wise flux values.
|
|
||||||
/// They are given in cell order, and for each cell there is
|
|
||||||
/// one value for each adjacent face (in the same order as the
|
|
||||||
/// cell-face topology of the grid). Positive values represent
|
|
||||||
/// fluxes out of the cell.
|
|
||||||
void faceFluxToCellFlux(const std::vector<double>& face_fluxes,
|
|
||||||
std::vector<double>& cell_fluxes)
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in TPFACompressiblePressureSolver::faceFluxToCellFlux(): "
|
|
||||||
"You must call assemble() (and solve the linear system) "
|
|
||||||
"prior to calling faceFluxToCellFlux().");
|
|
||||||
}
|
|
||||||
const UnstructuredGrid& g = *(grid_.c_grid());
|
|
||||||
int num_cells = g.number_of_cells;
|
|
||||||
cell_fluxes.resize(g.cell_facepos[num_cells]);
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
|
|
||||||
int face = g.cell_faces[hface];
|
|
||||||
bool pos = (g.face_cells[2*face] == cell);
|
|
||||||
cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Access the number of connections (faces) per cell. Deprecated, will be removed.
|
|
||||||
const std::vector<int>& numCellFaces() const
|
|
||||||
{
|
|
||||||
return ncf_;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
const std::vector<double>& faceTransmissibilities() const
|
|
||||||
{
|
|
||||||
return trans_;
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
// Disabling copy and assigment for now.
|
|
||||||
TPFACompressiblePressureSolver(const TPFACompressiblePressureSolver&);
|
|
||||||
TPFACompressiblePressureSolver& operator=(const TPFACompressiblePressureSolver&);
|
|
||||||
|
|
||||||
enum State { Uninitialized, Initialized, Assembled };
|
|
||||||
State state_;
|
|
||||||
|
|
||||||
// Solver data.
|
|
||||||
cfs_tpfa_data* data_;
|
|
||||||
// Grid.
|
|
||||||
GridAdapter grid_;
|
|
||||||
// Number of faces per cell.
|
|
||||||
std::vector<int> ncf_;
|
|
||||||
// Transmissibility storage.
|
|
||||||
std::vector<double> htrans_;
|
|
||||||
std::vector<double> trans_;
|
|
||||||
// Pore volumes.
|
|
||||||
std::vector<double> porevol_;
|
|
||||||
// Phase mobilities per face.
|
|
||||||
std::vector<double> phasemobf_;
|
|
||||||
// Gravity and capillary contributions (per face).
|
|
||||||
std::vector<double> gravcapf_;
|
|
||||||
// Gravity
|
|
||||||
double gravity_[3];
|
|
||||||
|
|
||||||
// Boundary conditions.
|
|
||||||
FlowBoundaryConditions *bc_;
|
|
||||||
|
|
||||||
// Well data
|
|
||||||
well_t wells_;
|
|
||||||
std::vector<int> well_connpos_storage_;
|
|
||||||
std::vector<int> well_cells_storage_;
|
|
||||||
well_control_t wctrl_;
|
|
||||||
std::vector<well_type> wctrl_type_storage_;
|
|
||||||
std::vector<well_control> wctrl_ctrl_storage_;
|
|
||||||
std::vector<double> wctrl_target_storage_;
|
|
||||||
struct completion_data wcompl_;
|
|
||||||
std::vector<double> well_prodind_storage_;
|
|
||||||
std::vector<double> well_gpot_storage_;
|
|
||||||
std::vector<double> well_A_storage_;
|
|
||||||
std::vector<double> well_phasemob_storage_;
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Initialize wells in solver structure.
|
|
||||||
/// @tparam Wells
|
|
||||||
/// This must conform to the SimpleWells concept.
|
|
||||||
/// @param w
|
|
||||||
/// The well object.
|
|
||||||
template <class Wells>
|
|
||||||
void initWells(const Wells& w)
|
|
||||||
{
|
|
||||||
int num_wells = w.numWells();
|
|
||||||
if (num_wells == 0) {
|
|
||||||
wells_.number_of_wells = 0;
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
wctrl_type_storage_.resize(num_wells);
|
|
||||||
wctrl_ctrl_storage_.resize(num_wells);
|
|
||||||
wctrl_target_storage_.resize(num_wells);
|
|
||||||
for (int i = 0; i < num_wells; ++i) {
|
|
||||||
wctrl_type_storage_[i] = (w.type(i) == Wells::Injector) ? INJECTOR : PRODUCER;
|
|
||||||
wctrl_ctrl_storage_[i] = (w.control(i) == Wells::Rate) ? RATE : BHP;
|
|
||||||
wctrl_target_storage_[i] = w.target(i);
|
|
||||||
int num_perf = w.numPerforations(i);
|
|
||||||
well_connpos_storage_.push_back(well_cells_storage_.size());
|
|
||||||
for (int j = 0; j < num_perf; ++j) {
|
|
||||||
well_cells_storage_.push_back(w.wellCell(i, j));
|
|
||||||
well_prodind_storage_.push_back(w.wellIndex(i, j));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
well_connpos_storage_.push_back(well_cells_storage_.size());
|
|
||||||
int tot_num_perf = well_prodind_storage_.size();
|
|
||||||
well_gpot_storage_.resize(tot_num_perf*3);
|
|
||||||
well_A_storage_.resize(3*3*tot_num_perf);
|
|
||||||
well_phasemob_storage_.resize(3*tot_num_perf);
|
|
||||||
// Setup 'wells_'
|
|
||||||
wells_.number_of_wells = num_wells;
|
|
||||||
wells_.well_connpos = &well_connpos_storage_[0];
|
|
||||||
wells_.well_cells = &well_cells_storage_[0];
|
|
||||||
// Setup 'wctrl_'
|
|
||||||
wctrl_.type = &wctrl_type_storage_[0];
|
|
||||||
wctrl_.ctrl = &wctrl_ctrl_storage_[0];
|
|
||||||
wctrl_.target = &wctrl_target_storage_[0];
|
|
||||||
// Setup 'wcompl_'
|
|
||||||
wcompl_.WI = &well_prodind_storage_[0];
|
|
||||||
wcompl_.gpot = &well_gpot_storage_[0];
|
|
||||||
wcompl_.A = &well_A_storage_[0];
|
|
||||||
wcompl_.phasemob = &well_phasemob_storage_[0];
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
void
|
|
||||||
gather_boundary_conditions(const FlowBCTypes* bctypes ,
|
|
||||||
const double* bcvalues)
|
|
||||||
{
|
|
||||||
if (bc_ == 0) {
|
|
||||||
bc_ = flow_conditions_construct(0);
|
|
||||||
}
|
|
||||||
else {
|
|
||||||
flow_conditions_clear(bc_);
|
|
||||||
}
|
|
||||||
|
|
||||||
int ok = bc_ != 0;
|
|
||||||
|
|
||||||
for (std::size_t i = 0, nf = grid_.numFaces(); ok && (i < nf); ++i) {
|
|
||||||
if (bctypes[ i ] == FBC_PRESSURE) {
|
|
||||||
ok = flow_conditions_append(BC_PRESSURE,
|
|
||||||
static_cast<int>(i),
|
|
||||||
bcvalues[ i ],
|
|
||||||
bc_);
|
|
||||||
}
|
|
||||||
else if (bctypes[ i ] == FBC_FLUX) {
|
|
||||||
ok = flow_conditions_append(BC_FLUX_TOTVOL,
|
|
||||||
static_cast<int>(i),
|
|
||||||
bcvalues[ i ],
|
|
||||||
bc_);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if (! ok) {
|
|
||||||
flow_conditions_destroy(bc_);
|
|
||||||
bc_ = 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
}; // class TPFACompressiblePressureSolver
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#endif // OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
@ -1,294 +0,0 @@
|
|||||||
/*
|
|
||||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
|
||||||
|
|
||||||
OPM is free software: you can redistribute it and/or modify
|
|
||||||
it under the terms of the GNU General Public License as published by
|
|
||||||
the Free Software Foundation, either version 3 of the License, or
|
|
||||||
(at your option) any later version.
|
|
||||||
|
|
||||||
OPM is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
You should have received a copy of the GNU General Public License
|
|
||||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
#define OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
|
|
||||||
|
|
||||||
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
|
|
||||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
|
||||||
#include <opm/core/linalg/sparse_sys.h>
|
|
||||||
#include <opm/core/pressure/flow_bc.h>
|
|
||||||
#include <opm/core/pressure/mimetic/mimetic.h> // for updating gpress
|
|
||||||
#include <opm/core/GridAdapter.hpp>
|
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
|
||||||
#include <stdexcept>
|
|
||||||
#include <cassert>
|
|
||||||
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Encapsulates the ifs_tpfa (= incompressible flow solver
|
|
||||||
/// two-point flux approximation) solver modules.
|
|
||||||
class TPFAPressureSolver
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
/// @brief
|
|
||||||
/// Default constructor, does nothing.
|
|
||||||
TPFAPressureSolver()
|
|
||||||
: state_(Uninitialized), data_(0)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Destructor.
|
|
||||||
~TPFAPressureSolver()
|
|
||||||
{
|
|
||||||
ifs_tpfa_destroy(data_);
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Initialize the solver's structures for a given grid (at some point also well pattern).
|
|
||||||
/// @tparam Grid This must conform to the SimpleGrid concept.
|
|
||||||
/// @param grid The grid object.
|
|
||||||
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
|
|
||||||
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
|
|
||||||
template <class Grid>
|
|
||||||
void init(const Grid& grid, const double* perm, const double* gravity)
|
|
||||||
{
|
|
||||||
// Build C grid structure.
|
|
||||||
grid_.init(grid);
|
|
||||||
|
|
||||||
// Build (empty for now) C well structure.
|
|
||||||
// well_t* w = 0;
|
|
||||||
|
|
||||||
// Initialize data.
|
|
||||||
data_ = ifs_tpfa_construct(grid_.c_grid(), 0);
|
|
||||||
if (!data_) {
|
|
||||||
throw std::runtime_error("Failed to initialize ifs_tpfa solver.");
|
|
||||||
}
|
|
||||||
|
|
||||||
// Compute half-transmissibilities, gravity contributions.
|
|
||||||
int num_cells = grid.numCells();
|
|
||||||
int ngconn = grid_.c_grid()->cell_facepos[num_cells];
|
|
||||||
gpress_.clear();
|
|
||||||
gpress_.resize(ngconn, 0.0);
|
|
||||||
ncf_.resize(num_cells);
|
|
||||||
typename Grid::Vector grav;
|
|
||||||
std::copy(gravity, gravity + Grid::dimension, &grav[0]);
|
|
||||||
int count = 0;
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
int num_local_faces = grid.numCellFaces(cell);
|
|
||||||
ncf_[cell] = num_local_faces;
|
|
||||||
typename Grid::Vector cc = grid.cellCentroid(cell);
|
|
||||||
for (int local_ix = 0; local_ix < num_local_faces; ++local_ix) {
|
|
||||||
int face = grid.cellFace(cell, local_ix);
|
|
||||||
typename Grid::Vector fc = grid.faceCentroid(face);
|
|
||||||
gpress_[count++] = grav*(fc - cc);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
assert(count == ngconn);
|
|
||||||
htrans_.resize(ngconn);
|
|
||||||
tpfa_htrans_compute(grid_.c_grid(), perm, &htrans_[0]);
|
|
||||||
state_ = Initialized;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
enum FlowBCTypes { FBC_UNSET, FBC_PRESSURE, FBC_FLUX };
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Assemble the sparse system.
|
|
||||||
/// You must call init() prior to calling assemble().
|
|
||||||
/// @param sources Source terms, one per cell. Positive numbers
|
|
||||||
/// are sources, negative are sinks.
|
|
||||||
/// @param total_mobilities Scalar total mobilities, one per cell.
|
|
||||||
/// @param omegas Gravity term, one per cell. In a multi-phase
|
|
||||||
/// flow setting this is equal to
|
|
||||||
/// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f]
|
|
||||||
/// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a
|
|
||||||
/// phase density and \f$\lambda_t\f$ is the total mobility.
|
|
||||||
void assemble(const std::vector<double>& sources,
|
|
||||||
const std::vector<double>& total_mobilities,
|
|
||||||
const std::vector<double>& omegas,
|
|
||||||
const std::vector<FlowBCTypes>& bctypes,
|
|
||||||
const std::vector<double>& bcvalues)
|
|
||||||
{
|
|
||||||
if (state_ == Uninitialized) {
|
|
||||||
throw std::runtime_error("Error in TPFAPressureSolver::assemble(): You must call init() prior to calling assemble().");
|
|
||||||
}
|
|
||||||
UnstructuredGrid* g = grid_.c_grid();
|
|
||||||
|
|
||||||
// Source terms from user.
|
|
||||||
double* src = const_cast<double*>(&sources[0]); // Ugly? Yes. Safe? I think so.
|
|
||||||
|
|
||||||
// All well related things are zero.
|
|
||||||
// well_control_t* wctrl = 0;
|
|
||||||
// double* WI = 0;
|
|
||||||
// double* wdp = 0;
|
|
||||||
|
|
||||||
// Compute effective transmissibilities.
|
|
||||||
eff_trans_.resize(grid_.numFaces());
|
|
||||||
tpfa_eff_trans_compute(g, &total_mobilities[0], &htrans_[0], &eff_trans_[0]);
|
|
||||||
|
|
||||||
// Update gravity term.
|
|
||||||
gpress_omegaweighted_.resize(gpress_.size());
|
|
||||||
mim_ip_density_update(g->number_of_cells, g->cell_facepos, &omegas[0],
|
|
||||||
&gpress_[0], &gpress_omegaweighted_[0]);
|
|
||||||
|
|
||||||
|
|
||||||
// Zero the linalg structures.
|
|
||||||
csrmatrix_zero(data_->A);
|
|
||||||
for (std::size_t i = 0; i < data_->A->m; i++) {
|
|
||||||
data_->b[i] = 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
forces_.src = &src[0];
|
|
||||||
forces_.bc = 0;
|
|
||||||
forces_.W = 0;
|
|
||||||
|
|
||||||
// Assemble the embedded linear system.
|
|
||||||
int ok = ifs_tpfa_assemble(g, &forces_, &eff_trans_[0], &gpress_[0], data_);
|
|
||||||
if (!ok) {
|
|
||||||
THROW("Assembly of pressure system failed.");
|
|
||||||
}
|
|
||||||
state_ = Assembled;
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Encapsulate a sparse linear system in CSR format.
|
|
||||||
struct LinearSystem
|
|
||||||
{
|
|
||||||
int n;
|
|
||||||
int nnz;
|
|
||||||
int* ia;
|
|
||||||
int* ja;
|
|
||||||
double* sa;
|
|
||||||
double* b;
|
|
||||||
double* x;
|
|
||||||
};
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Access the linear system assembled.
|
|
||||||
/// You must call assemble() prior to calling linearSystem().
|
|
||||||
/// @param[out] s The linear system encapsulation to modify.
|
|
||||||
/// After this call, s will point to linear system structures
|
|
||||||
/// that are owned and allocated internally.
|
|
||||||
void linearSystem(LinearSystem& s)
|
|
||||||
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in TPFAPressureSolver::linearSystem(): "
|
|
||||||
"You must call assemble() prior to calling linearSystem().");
|
|
||||||
}
|
|
||||||
s.n = data_->A->m;
|
|
||||||
s.nnz = data_->A->nnz;
|
|
||||||
s.ia = data_->A->ia;
|
|
||||||
s.ja = data_->A->ja;
|
|
||||||
s.sa = data_->A->sa;
|
|
||||||
s.b = data_->b;
|
|
||||||
s.x = data_->x;
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Compute cell pressures and face fluxes.
|
|
||||||
/// You must call assemble() (and solve the linear system accessed
|
|
||||||
/// by calling linearSystem()) prior to calling
|
|
||||||
/// computePressuresAndFluxes().
|
|
||||||
/// @param[out] cell_pressures Cell pressure values.
|
|
||||||
/// @param[out] face_areas Face flux values.
|
|
||||||
void computePressuresAndFluxes(std::vector<double>& cell_pressures,
|
|
||||||
std::vector<double>& face_fluxes)
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in TPFAPressureSolver::computePressuresAndFluxes(): "
|
|
||||||
"You must call assemble() (and solve the linear system) "
|
|
||||||
"prior to calling computePressuresAndFluxes().");
|
|
||||||
}
|
|
||||||
int num_cells = grid_.c_grid()->number_of_cells;
|
|
||||||
int num_faces = grid_.c_grid()->number_of_faces;
|
|
||||||
cell_pressures.clear();
|
|
||||||
cell_pressures.resize(num_cells, 0.0);
|
|
||||||
face_fluxes.clear();
|
|
||||||
face_fluxes.resize(num_faces, 0.0);
|
|
||||||
|
|
||||||
ifs_tpfa_solution soln = { 0 };
|
|
||||||
soln.cell_press = &cell_pressures[0];
|
|
||||||
soln.face_flux = &face_fluxes [0];
|
|
||||||
|
|
||||||
ifs_tpfa_press_flux(grid_.c_grid(), &forces_, &eff_trans_[0],
|
|
||||||
data_, &soln);
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Compute cell fluxes from face fluxes.
|
|
||||||
/// You must call assemble() (and solve the linear system accessed
|
|
||||||
/// by calling linearSystem()) prior to calling
|
|
||||||
/// faceFluxToCellFlux().
|
|
||||||
/// @param face_fluxes
|
|
||||||
/// @param face_areas Face flux values (usually output from computePressuresAndFluxes()).
|
|
||||||
/// @param[out] cell_fluxes Cell-wise flux values.
|
|
||||||
/// They are given in cell order, and for each cell there is
|
|
||||||
/// one value for each adjacent face (in the same order as the
|
|
||||||
/// cell-face topology of the grid). Positive values represent
|
|
||||||
/// fluxes out of the cell.
|
|
||||||
void faceFluxToCellFlux(const std::vector<double>& face_fluxes,
|
|
||||||
std::vector<double>& cell_fluxes)
|
|
||||||
{
|
|
||||||
if (state_ != Assembled) {
|
|
||||||
throw std::runtime_error("Error in TPFAPressureSolver::faceFluxToCellFlux(): "
|
|
||||||
"You must call assemble() (and solve the linear system) "
|
|
||||||
"prior to calling faceFluxToCellFlux().");
|
|
||||||
}
|
|
||||||
const UnstructuredGrid& g = *(grid_.c_grid());
|
|
||||||
int num_cells = g.number_of_cells;
|
|
||||||
cell_fluxes.resize(g.cell_facepos[num_cells]);
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
|
|
||||||
int face = g.cell_faces[hface];
|
|
||||||
bool pos = (g.face_cells[2*face] == cell);
|
|
||||||
cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/// @brief
|
|
||||||
/// Access the number of connections (faces) per cell. Deprecated, will be removed.
|
|
||||||
const std::vector<int>& numCellFaces()
|
|
||||||
{
|
|
||||||
return ncf_;
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
// Disabling copy and assigment for now.
|
|
||||||
TPFAPressureSolver(const TPFAPressureSolver&);
|
|
||||||
TPFAPressureSolver& operator=(const TPFAPressureSolver&);
|
|
||||||
|
|
||||||
enum State { Uninitialized, Initialized, Assembled };
|
|
||||||
State state_;
|
|
||||||
|
|
||||||
// Solver data.
|
|
||||||
ifs_tpfa_data* data_;
|
|
||||||
ifs_tpfa_forces forces_;
|
|
||||||
// Grid.
|
|
||||||
GridAdapter grid_;
|
|
||||||
// Number of faces per cell.
|
|
||||||
std::vector<int> ncf_;
|
|
||||||
// Transmissibility storage.
|
|
||||||
std::vector<double> htrans_;
|
|
||||||
std::vector<double> eff_trans_;
|
|
||||||
// Gravity contributions.
|
|
||||||
std::vector<double> gpress_;
|
|
||||||
std::vector<double> gpress_omegaweighted_;
|
|
||||||
// Total mobilities.
|
|
||||||
std::vector<double> totmob_;
|
|
||||||
// Gravity coefficients (\omega = sum_{i = 1}^{num phases}f_i \rho_i[TODO: check this]).
|
|
||||||
std::vector<double> omega_;
|
|
||||||
};
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
#endif // OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED
|
|
||||||
@ -45,7 +45,7 @@
|
|||||||
#include <opm/core/simulator/TwophaseState.hpp>
|
#include <opm/core/simulator/TwophaseState.hpp>
|
||||||
#include <opm/core/simulator/WellState.hpp>
|
#include <opm/core/simulator/WellState.hpp>
|
||||||
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
|
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
|
||||||
#include <opm/core/transport/TransportSolverTwophaseImplicit.hpp>
|
#include <opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp>
|
||||||
#include <boost/filesystem.hpp>
|
#include <boost/filesystem.hpp>
|
||||||
#include <boost/scoped_ptr.hpp>
|
#include <boost/scoped_ptr.hpp>
|
||||||
#include <boost/lexical_cast.hpp>
|
#include <boost/lexical_cast.hpp>
|
||||||
|
|||||||
@ -17,7 +17,7 @@
|
|||||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
#include <opm/core/transport/reorder/DGBasis.hpp>
|
#include <opm/core/tof/DGBasis.hpp>
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
#include <numeric>
|
#include <numeric>
|
||||||
@ -19,8 +19,8 @@
|
|||||||
|
|
||||||
#include <opm/core/grid/CellQuadrature.hpp>
|
#include <opm/core/grid/CellQuadrature.hpp>
|
||||||
#include <opm/core/grid/FaceQuadrature.hpp>
|
#include <opm/core/grid/FaceQuadrature.hpp>
|
||||||
#include <opm/core/transport/reorder/TofDiscGalReorder.hpp>
|
#include <opm/core/tof/TofDiscGalReorder.hpp>
|
||||||
#include <opm/core/transport/reorder/DGBasis.hpp>
|
#include <opm/core/tof/DGBasis.hpp>
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
#include <opm/core/utility/VelocityInterpolation.hpp>
|
#include <opm/core/utility/VelocityInterpolation.hpp>
|
||||||
@ -17,7 +17,7 @@
|
|||||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
#include <opm/core/transport/reorder/TofReorder.hpp>
|
#include <opm/core/tof/TofReorder.hpp>
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
#include <algorithm>
|
#include <algorithm>
|
||||||
@ -36,7 +36,7 @@
|
|||||||
#ifndef OPM_IMPLICITTRANSPORT_HPP_HEADER
|
#ifndef OPM_IMPLICITTRANSPORT_HPP_HEADER
|
||||||
#define OPM_IMPLICITTRANSPORT_HPP_HEADER
|
#define OPM_IMPLICITTRANSPORT_HPP_HEADER
|
||||||
|
|
||||||
#include <opm/core/transport/ImplicitAssembly.hpp>
|
#include <opm/core/transport/implicit/ImplicitAssembly.hpp>
|
||||||
|
|
||||||
namespace Opm {
|
namespace Opm {
|
||||||
namespace ImplicitTransportDetails {
|
namespace ImplicitTransportDetails {
|
||||||
@ -66,5 +66,5 @@ namespace Opm{
|
|||||||
};
|
};
|
||||||
}
|
}
|
||||||
|
|
||||||
#include <opm/core/transport/SimpleFluid2pWrappingProps_impl.hpp>
|
#include <opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp>
|
||||||
#endif // SIMPLEFLUID2PWRAPPINGPROPS_HPP
|
#endif // SIMPLEFLUID2PWRAPPINGPROPS_HPP
|
||||||
@ -28,7 +28,7 @@
|
|||||||
#ifndef SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP
|
#ifndef SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP
|
||||||
#define SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP
|
#define SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP
|
||||||
|
|
||||||
#include <opm/core/transport/SimpleFluid2pWrappingProps.hpp>
|
#include <opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp>
|
||||||
#include <cassert>
|
#include <cassert>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
namespace Opm{
|
namespace Opm{
|
||||||
@ -27,7 +27,7 @@
|
|||||||
*/
|
*/
|
||||||
|
|
||||||
|
|
||||||
#include <opm/core/transport/TransportSolverTwophaseImplicit.hpp>
|
#include <opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp>
|
||||||
#include <opm/core/simulator/TwophaseState.hpp>
|
#include <opm/core/simulator/TwophaseState.hpp>
|
||||||
#include <opm/core/utility/miscUtilities.hpp>
|
#include <opm/core/utility/miscUtilities.hpp>
|
||||||
|
|
||||||
@ -21,24 +21,22 @@
|
|||||||
#ifndef OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED
|
#ifndef OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED
|
||||||
#define OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED
|
#define OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED
|
||||||
|
|
||||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
|
||||||
#include <opm/core/transport/TransportSolverTwophaseInterface.hpp>
|
#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/props/IncompPropertiesInterface.hpp>
|
||||||
#include <opm/core/transport/SimpleFluid2pWrappingProps.hpp>
|
|
||||||
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
|
|
||||||
#include <opm/core/transport/ImplicitTransport.hpp>
|
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/linalg/LinearSolverFactory.hpp>
|
#include <opm/core/linalg/LinearSolverFactory.hpp>
|
||||||
|
|
||||||
#include <opm/core/transport/transport_source.h>
|
|
||||||
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
|
|
||||||
#include <opm/core/transport/NormSupport.hpp>
|
|
||||||
#include <opm/core/transport/ImplicitAssembly.hpp>
|
|
||||||
#include <opm/core/transport/ImplicitTransport.hpp>
|
|
||||||
#include <opm/core/transport/JacobianSystem.hpp>
|
|
||||||
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
|
|
||||||
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
|
|
||||||
|
|
||||||
#include <boost/scoped_ptr.hpp>
|
#include <boost/scoped_ptr.hpp>
|
||||||
#include <vector>
|
#include <vector>
|
||||||
|
|
||||||
@ -36,7 +36,7 @@
|
|||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
#include <string.h>
|
#include <string.h>
|
||||||
|
|
||||||
#include <opm/core/transport/transport_source.h>
|
#include <opm/core/transport/implicit/transport_source.h>
|
||||||
|
|
||||||
|
|
||||||
/* ---------------------------------------------------------------------- */
|
/* ---------------------------------------------------------------------- */
|
||||||
@ -8,7 +8,7 @@
|
|||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
|
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/transport/spu_explicit.h>
|
#include <opm/core/transport/minimal/spu_explicit.h>
|
||||||
|
|
||||||
|
|
||||||
/* Twophase mobility-weighted upwind */
|
/* Twophase mobility-weighted upwind */
|
||||||
@ -11,7 +11,7 @@
|
|||||||
#include <math.h>
|
#include <math.h>
|
||||||
|
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/transport/spu_implicit.h>
|
#include <opm/core/transport/minimal/spu_implicit.h>
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -1,308 +0,0 @@
|
|||||||
/*
|
|
||||||
Copyright (C) 2012 (c) Jostein R. Natvig <jostein natvig at gmail.com>
|
|
||||||
|
|
||||||
Permission is hereby granted, free of charge, to any person obtaining a copy of
|
|
||||||
this software and associated documentation files (the "Software"), to deal in
|
|
||||||
the Software without restriction, including without limitation the rights to
|
|
||||||
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
|
|
||||||
of the Software, and to permit persons to whom the Software is furnished to do
|
|
||||||
so, subject to the following conditions:
|
|
||||||
|
|
||||||
The above copyright notice and this permission notice shall be included in all
|
|
||||||
copies or substantial portions of the Software.
|
|
||||||
|
|
||||||
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
||||||
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
||||||
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
||||||
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
||||||
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
||||||
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
||||||
SOFTWARE.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include <string.h>
|
|
||||||
#include <assert.h>
|
|
||||||
#include <math.h>
|
|
||||||
|
|
||||||
#include <opm/core/transport/reorder/nlsolvers.h>
|
|
||||||
|
|
||||||
static const char no_root_str[]=
|
|
||||||
" In %s:\n"
|
|
||||||
" With G(%5f) =% 5f, G(%5f) =% 5f, G(x) is not bracketed!\n";
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
double
|
|
||||||
find_zero (double (*f)(double, void*), void *data, struct NonlinearSolverCtrl *ctrl)
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
{
|
|
||||||
double zero;
|
|
||||||
switch (ctrl->method) {
|
|
||||||
default:
|
|
||||||
case BISECTION:
|
|
||||||
zero = bisection (f, data, ctrl);
|
|
||||||
break;
|
|
||||||
|
|
||||||
case RIDDERS:
|
|
||||||
zero = ridders (f, data, ctrl);
|
|
||||||
break;
|
|
||||||
|
|
||||||
case REGULAFALSI:
|
|
||||||
zero = regulafalsi(f, data, ctrl);
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
return zero;
|
|
||||||
}
|
|
||||||
|
|
||||||
/* Start with bracket [0,1] with G(0)*G(1)<0. */
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
double
|
|
||||||
ridders (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctrl)
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
{
|
|
||||||
double G0, G1, G2, G3;
|
|
||||||
double s0, s1, s2, s3;
|
|
||||||
double swap, sgn, root;
|
|
||||||
|
|
||||||
ctrl->iterations = 0;
|
|
||||||
|
|
||||||
s2 = ctrl->initialguess;
|
|
||||||
G2 = G(s2, data);
|
|
||||||
if (fabs(G2) < ctrl->nltolerance) { return s2; }
|
|
||||||
|
|
||||||
s0 = ctrl->min_bracket;
|
|
||||||
G0 = G(s0, data);
|
|
||||||
if (fabs(G0) < ctrl->nltolerance) { return s0; }
|
|
||||||
|
|
||||||
s1 = ctrl->max_bracket;
|
|
||||||
G1 = G(s1, data);
|
|
||||||
if (fabs(G1) < ctrl->nltolerance) { return s1; }
|
|
||||||
|
|
||||||
if (G0*G1 > 0)
|
|
||||||
{
|
|
||||||
printf(no_root_str, "ridder", s0, G0, s1, G1);
|
|
||||||
return -1.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (G0>G1)
|
|
||||||
{
|
|
||||||
swap = G0;
|
|
||||||
G0 = G1;
|
|
||||||
G1 = swap;
|
|
||||||
}
|
|
||||||
|
|
||||||
s3 = 0;
|
|
||||||
G3 = 10;
|
|
||||||
|
|
||||||
while ( (fabs(G3) > ctrl->nltolerance) &&
|
|
||||||
(ctrl->iterations++ < ctrl->maxiterations))
|
|
||||||
{
|
|
||||||
/* find zero crossing of line segment [(s0,G0), (s1,G1)] */
|
|
||||||
root = sqrt(G2*G2 - G0*G1);
|
|
||||||
|
|
||||||
sgn = G0>G1 ? 1.0 : -1.0;
|
|
||||||
s3 = s2 + ( s2-s0 )*sgn*G2/root;
|
|
||||||
G3 = G(s3, data);
|
|
||||||
|
|
||||||
|
|
||||||
/* if G2*G3<0 */
|
|
||||||
if (G2*G3 <= 0.0)
|
|
||||||
{
|
|
||||||
if (G2 > G3)
|
|
||||||
{
|
|
||||||
s0 = s3;
|
|
||||||
G0 = G3;
|
|
||||||
s1 = s2;
|
|
||||||
G1 = G2;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
s0 = s2;
|
|
||||||
G0 = G2;
|
|
||||||
s1 = s3;
|
|
||||||
G1 = G3;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
else if (G0*G3 <= 0.0)
|
|
||||||
{
|
|
||||||
s1 = s3;
|
|
||||||
G1 = G3;
|
|
||||||
}
|
|
||||||
else if (G1*G3 <= 0.0)
|
|
||||||
{
|
|
||||||
s0 = s3;
|
|
||||||
G0 = G3;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
printf("In ridder:\nG0=%10.10f, G1=%10.10f, "
|
|
||||||
"G3=%10.10f\n", G0, G1, G3);
|
|
||||||
}
|
|
||||||
s2 = 0.5*(s0+s1);
|
|
||||||
G2 = G(s2, data);
|
|
||||||
}
|
|
||||||
|
|
||||||
ctrl->residual = G3;
|
|
||||||
return s3;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/* Start with bracket [0,1] with G(0)*G(1)<0. Search by finding zero crossing
|
|
||||||
sN of line segment [(sL,GL), (sR,GR)]. Set SL=sN if G(sN<0), sR=sN
|
|
||||||
otherwise.*/
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
double
|
|
||||||
regulafalsi (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctrl)
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
{
|
|
||||||
double Gn, G0, G1;
|
|
||||||
double sn, s0, s1;
|
|
||||||
double swap, gamma_pegasus;
|
|
||||||
|
|
||||||
ctrl->iterations = 0;
|
|
||||||
|
|
||||||
sn = ctrl->initialguess;
|
|
||||||
Gn = G(sn, data);
|
|
||||||
if (fabs(Gn) < ctrl->nltolerance) { return sn; }
|
|
||||||
|
|
||||||
/* Initial guess is interval [s0,s1] */
|
|
||||||
s0 = ctrl->min_bracket;
|
|
||||||
G0 = G(s0, data);
|
|
||||||
if (fabs(G0) < ctrl->nltolerance) { return s0; }
|
|
||||||
|
|
||||||
s1 = ctrl->max_bracket;
|
|
||||||
G1 = G(s1, data);
|
|
||||||
if (fabs(G1) < ctrl->nltolerance) { return s1; }
|
|
||||||
|
|
||||||
if (G0*G1 > 0)
|
|
||||||
{
|
|
||||||
printf(no_root_str, "regulafalsi", s0, G0, s1, G1);
|
|
||||||
return -1.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (G0>G1)
|
|
||||||
{
|
|
||||||
swap = G0;
|
|
||||||
G0 = G1;
|
|
||||||
G1 = swap;
|
|
||||||
}
|
|
||||||
|
|
||||||
while ( (fabs(Gn) > ctrl->nltolerance) &&
|
|
||||||
(ctrl->iterations++ < ctrl->maxiterations))
|
|
||||||
{
|
|
||||||
|
|
||||||
#if 0
|
|
||||||
/* Unmodified Regula-Falsi */
|
|
||||||
/* maintain bracket with G1>G0 */
|
|
||||||
if ( Gn>0 )
|
|
||||||
{
|
|
||||||
G1 = Gn;
|
|
||||||
s1 = sn;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
G0 = Gn;
|
|
||||||
s0 = sn;
|
|
||||||
}
|
|
||||||
#else
|
|
||||||
/* Modified Regula-Falsi*/
|
|
||||||
if ((Gn>0.0)==(G0>0.0))
|
|
||||||
{
|
|
||||||
s0 = s1;
|
|
||||||
G0 = G1;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
/* const double gamma_illinois = 0.5; */
|
|
||||||
gamma_pegasus = G1/(G1+Gn);
|
|
||||||
G0 *= gamma_pegasus;
|
|
||||||
}
|
|
||||||
s1 = sn;
|
|
||||||
G1 = Gn;
|
|
||||||
#endif
|
|
||||||
|
|
||||||
sn = s0 - (s1-s0)*G0/(G1-G0);
|
|
||||||
Gn = G(sn, data);
|
|
||||||
}
|
|
||||||
|
|
||||||
ctrl->residual = Gn;
|
|
||||||
return sn;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/* Start with bracket [0,1] with G(0)*G(1)<0. Search by finding sN=0.5(sL+sR).
|
|
||||||
Set SL=sN if G(sN<0), sR=sN otherwise.*/
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
double
|
|
||||||
bisection (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctrl)
|
|
||||||
/*---------------------------------------------------------------------------*/
|
|
||||||
{
|
|
||||||
double Gn, G0, G1;
|
|
||||||
double sn, s0, s1;
|
|
||||||
double swap;
|
|
||||||
|
|
||||||
ctrl->iterations = 0;
|
|
||||||
|
|
||||||
sn = ctrl->initialguess;
|
|
||||||
Gn = G(sn, data);
|
|
||||||
if (fabs(Gn) < ctrl->nltolerance) { return sn; }
|
|
||||||
|
|
||||||
/* Initial guess is interval [s0,s1] */
|
|
||||||
s0 = ctrl->max_bracket;
|
|
||||||
G0 = G(s0, data);
|
|
||||||
if (fabs(G0) < ctrl->nltolerance) { return s0; }
|
|
||||||
|
|
||||||
s1 = ctrl->max_bracket;
|
|
||||||
G1 = G(s1, data);
|
|
||||||
if (fabs(G1) < ctrl->nltolerance) { return s1; }
|
|
||||||
|
|
||||||
if (G0*G1 > 0.0)
|
|
||||||
{
|
|
||||||
printf(no_root_str, "bisection", s0, G0, s1, G1);
|
|
||||||
return -1.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (G0>G1)
|
|
||||||
{
|
|
||||||
swap = G0;
|
|
||||||
G0 = G1;
|
|
||||||
G1 = swap;
|
|
||||||
}
|
|
||||||
|
|
||||||
while ( (fabs(Gn)>ctrl->nltolerance) &&
|
|
||||||
(ctrl->iterations++ < ctrl->maxiterations) )
|
|
||||||
{
|
|
||||||
if ( Gn>0 )
|
|
||||||
{
|
|
||||||
G1 = Gn;
|
|
||||||
s1 = sn;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
G0 = Gn;
|
|
||||||
s0 = sn;
|
|
||||||
}
|
|
||||||
|
|
||||||
sn = 0.5*(s0+s1);
|
|
||||||
Gn = G(sn, data);
|
|
||||||
}
|
|
||||||
if (ctrl->iterations >= ctrl->maxiterations)
|
|
||||||
{
|
|
||||||
printf("Warning: convergence criterion not met\n");
|
|
||||||
}
|
|
||||||
ctrl->residual = Gn;
|
|
||||||
return sn;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
/* Local Variables: */
|
|
||||||
/* c-basic-offset:4 */
|
|
||||||
/* End: */
|
|
||||||
@ -1,55 +0,0 @@
|
|||||||
/*
|
|
||||||
Copyright (C) 2012 (c) Jostein R. Natvig <jostein natvig at gmail.com>
|
|
||||||
|
|
||||||
Permission is hereby granted, free of charge, to any person obtaining a copy of
|
|
||||||
this software and associated documentation files (the "Software"), to deal in
|
|
||||||
the Software without restriction, including without limitation the rights to
|
|
||||||
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
|
|
||||||
of the Software, and to permit persons to whom the Software is furnished to do
|
|
||||||
so, subject to the following conditions:
|
|
||||||
|
|
||||||
The above copyright notice and this permission notice shall be included in all
|
|
||||||
copies or substantial portions of the Software.
|
|
||||||
|
|
||||||
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
||||||
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
||||||
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
||||||
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
||||||
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
||||||
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
||||||
SOFTWARE.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef NLSOLVERS_H
|
|
||||||
#define NLSOLVERS_H
|
|
||||||
|
|
||||||
#ifdef __cplusplus
|
|
||||||
extern "C" {
|
|
||||||
#endif
|
|
||||||
|
|
||||||
struct NonlinearSolverCtrl
|
|
||||||
{
|
|
||||||
enum Method {RIDDERS, REGULAFALSI, BISECTION} method;
|
|
||||||
double nltolerance;
|
|
||||||
int maxiterations;
|
|
||||||
double min_bracket;
|
|
||||||
double max_bracket;
|
|
||||||
double initialguess;
|
|
||||||
int iterations; /* set by solver */
|
|
||||||
double residual; /* set by solver */
|
|
||||||
};
|
|
||||||
|
|
||||||
double find_zero (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctrl);
|
|
||||||
double bisection (double (*)(double, void*), void*, struct NonlinearSolverCtrl *ctrl);
|
|
||||||
double ridders (double (*)(double, void*), void*, struct NonlinearSolverCtrl *ctrl);
|
|
||||||
double regulafalsi (double (*)(double, void*), void*, struct NonlinearSolverCtrl *ctrl);
|
|
||||||
|
|
||||||
#ifdef __cplusplus
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#endif /* NLSOLVERS_H */
|
|
||||||
|
|
||||||
/* Local Variables: */
|
|
||||||
/* c-basic-offset:4 */
|
|
||||||
/* End: */
|
|
||||||
Loading…
Reference in New Issue
Block a user