diff --git a/opm/core/pressure/HybridPressureSolver.hpp b/opm/core/pressure/HybridPressureSolver.hpp deleted file mode 100644 index 6e7a934e..00000000 --- a/opm/core/pressure/HybridPressureSolver.hpp +++ /dev/null @@ -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 . -*/ - - -#ifndef OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED -#define OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED - -#include -#include -#include -#include -#include - - -/// @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 - 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(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& sources, - const std::vector& total_mobilities, - const std::vector& omegas, - const std::vector& bctypes, - const std::vector& 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::size_type>(grid_.numFaces())); - - FlowBoundaryConditions *bc = gather_boundary_conditions(bctypes, bcvalues); - - // Source terms from user. - double* src = const_cast(&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& cell_pressures, - std::vector& 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& face_fluxes, - std::vector& 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& 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 ncf_; - // B^{-1} storage. - std::vector Binv_; - std::vector Binv_mobilityweighted_; - // Gravity contributions. - std::vector gpress_; - std::vector gpress_omegaweighted_; - - - FlowBoundaryConditions* - gather_boundary_conditions(const std::vector& bctypes , - const std::vector& bcvalues) - { - FlowBoundaryConditions* fbc = flow_conditions_construct(0); - - int ok = fbc != 0; - std::vector::size_type i; - - for (i = 0; ok && (i < bctypes.size()); ++i) { - if (bctypes[ i ] == FBC_PRESSURE) { - ok = flow_conditions_append(BC_PRESSURE, - static_cast(i), - bcvalues[ i ], - fbc); - } - else if (bctypes[ i ] == FBC_FLUX) { - ok = flow_conditions_append(BC_FLUX_TOTVOL, - static_cast(i), - bcvalues[ i ], - fbc); - } - } - - return fbc; - } - - -}; // class HybridPressureSolver - - -#endif // OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED diff --git a/opm/core/pressure/TPFACompressiblePressureSolver.hpp b/opm/core/pressure/TPFACompressiblePressureSolver.hpp deleted file mode 100644 index d38d8260..00000000 --- a/opm/core/pressure/TPFACompressiblePressureSolver.hpp +++ /dev/null @@ -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 . -*/ - -#ifndef OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED -#define OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED - - -#include -#include -#include -#include -#include -#include -#include -#include - - - -/// @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 - 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 - 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(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(totcompr ) , - const_cast(voldiscr ) , - const_cast(cellA ) , - const_cast(faceA ) , - const_cast(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& cell_pressures, - std::vector& face_pressures, - std::vector& face_fluxes, - std::vector& well_pressures, - std::vector& 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(faceA) , - const_cast(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& face_fluxes, - std::vector& 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& numCellFaces() const - { - return ncf_; - } - - - const std::vector& 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 ncf_; - // Transmissibility storage. - std::vector htrans_; - std::vector trans_; - // Pore volumes. - std::vector porevol_; - // Phase mobilities per face. - std::vector phasemobf_; - // Gravity and capillary contributions (per face). - std::vector gravcapf_; - // Gravity - double gravity_[3]; - - // Boundary conditions. - FlowBoundaryConditions *bc_; - - // Well data - well_t wells_; - std::vector well_connpos_storage_; - std::vector well_cells_storage_; - well_control_t wctrl_; - std::vector wctrl_type_storage_; - std::vector wctrl_ctrl_storage_; - std::vector wctrl_target_storage_; - struct completion_data wcompl_; - std::vector well_prodind_storage_; - std::vector well_gpot_storage_; - std::vector well_A_storage_; - std::vector well_phasemob_storage_; - - - /// @brief - /// Initialize wells in solver structure. - /// @tparam Wells - /// This must conform to the SimpleWells concept. - /// @param w - /// The well object. - template - 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(i), - bcvalues[ i ], - bc_); - } - else if (bctypes[ i ] == FBC_FLUX) { - ok = flow_conditions_append(BC_FLUX_TOTVOL, - static_cast(i), - bcvalues[ i ], - bc_); - } - } - - if (! ok) { - flow_conditions_destroy(bc_); - bc_ = 0; - } - } - - -}; // class TPFACompressiblePressureSolver - - - - -#endif // OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED diff --git a/opm/core/pressure/TPFAPressureSolver.hpp b/opm/core/pressure/TPFAPressureSolver.hpp deleted file mode 100644 index 33f7e33a..00000000 --- a/opm/core/pressure/TPFAPressureSolver.hpp +++ /dev/null @@ -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 . -*/ - -#ifndef OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED -#define OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED - - -#include -#include -#include -#include -#include // for updating gpress -#include -#include -#include -#include - - -/// @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 - 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& sources, - const std::vector& total_mobilities, - const std::vector& omegas, - const std::vector& bctypes, - const std::vector& 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(&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& cell_pressures, - std::vector& 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& face_fluxes, - std::vector& 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& 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 ncf_; - // Transmissibility storage. - std::vector htrans_; - std::vector eff_trans_; - // Gravity contributions. - std::vector gpress_; - std::vector gpress_omegaweighted_; - // Total mobilities. - std::vector totmob_; - // Gravity coefficients (\omega = sum_{i = 1}^{num phases}f_i \rho_i[TODO: check this]). - std::vector omega_; -}; - - - -#endif // OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED diff --git a/opm/core/simulator/SimulatorIncompTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp index a3f99625..5b60a8f6 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -45,7 +45,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/opm/core/transport/reorder/DGBasis.cpp b/opm/core/tof/DGBasis.cpp similarity index 99% rename from opm/core/transport/reorder/DGBasis.cpp rename to opm/core/tof/DGBasis.cpp index 8ffb95ec..b2ce5eaa 100644 --- a/opm/core/transport/reorder/DGBasis.cpp +++ b/opm/core/tof/DGBasis.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include diff --git a/opm/core/transport/reorder/DGBasis.hpp b/opm/core/tof/DGBasis.hpp similarity index 100% rename from opm/core/transport/reorder/DGBasis.hpp rename to opm/core/tof/DGBasis.hpp diff --git a/opm/core/transport/reorder/TofDiscGalReorder.cpp b/opm/core/tof/TofDiscGalReorder.cpp similarity index 99% rename from opm/core/transport/reorder/TofDiscGalReorder.cpp rename to opm/core/tof/TofDiscGalReorder.cpp index 10af155b..589393ae 100644 --- a/opm/core/transport/reorder/TofDiscGalReorder.cpp +++ b/opm/core/tof/TofDiscGalReorder.cpp @@ -19,8 +19,8 @@ #include #include -#include -#include +#include +#include #include #include #include diff --git a/opm/core/transport/reorder/TofDiscGalReorder.hpp b/opm/core/tof/TofDiscGalReorder.hpp similarity index 100% rename from opm/core/transport/reorder/TofDiscGalReorder.hpp rename to opm/core/tof/TofDiscGalReorder.hpp diff --git a/opm/core/transport/reorder/TofReorder.cpp b/opm/core/tof/TofReorder.cpp similarity index 99% rename from opm/core/transport/reorder/TofReorder.cpp rename to opm/core/tof/TofReorder.cpp index 8002adc4..48359343 100644 --- a/opm/core/transport/reorder/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include diff --git a/opm/core/transport/reorder/TofReorder.hpp b/opm/core/tof/TofReorder.hpp similarity index 100% rename from opm/core/transport/reorder/TofReorder.hpp rename to opm/core/tof/TofReorder.hpp diff --git a/opm/core/transport/CSRMatrixBlockAssembler.hpp b/opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp similarity index 100% rename from opm/core/transport/CSRMatrixBlockAssembler.hpp rename to opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp diff --git a/opm/core/transport/CSRMatrixUmfpackSolver.hpp b/opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp similarity index 100% rename from opm/core/transport/CSRMatrixUmfpackSolver.hpp rename to opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp diff --git a/opm/core/transport/GravityColumnSolver.hpp b/opm/core/transport/implicit/GravityColumnSolver.hpp similarity index 100% rename from opm/core/transport/GravityColumnSolver.hpp rename to opm/core/transport/implicit/GravityColumnSolver.hpp diff --git a/opm/core/transport/GravityColumnSolver_impl.hpp b/opm/core/transport/implicit/GravityColumnSolver_impl.hpp similarity index 100% rename from opm/core/transport/GravityColumnSolver_impl.hpp rename to opm/core/transport/implicit/GravityColumnSolver_impl.hpp diff --git a/opm/core/transport/ImplicitAssembly.hpp b/opm/core/transport/implicit/ImplicitAssembly.hpp similarity index 100% rename from opm/core/transport/ImplicitAssembly.hpp rename to opm/core/transport/implicit/ImplicitAssembly.hpp diff --git a/opm/core/transport/ImplicitTransport.hpp b/opm/core/transport/implicit/ImplicitTransport.hpp similarity index 99% rename from opm/core/transport/ImplicitTransport.hpp rename to opm/core/transport/implicit/ImplicitTransport.hpp index 7d5cb6e0..3bdd8122 100644 --- a/opm/core/transport/ImplicitTransport.hpp +++ b/opm/core/transport/implicit/ImplicitTransport.hpp @@ -36,7 +36,7 @@ #ifndef OPM_IMPLICITTRANSPORT_HPP_HEADER #define OPM_IMPLICITTRANSPORT_HPP_HEADER -#include +#include namespace Opm { namespace ImplicitTransportDetails { diff --git a/opm/core/transport/JacobianSystem.hpp b/opm/core/transport/implicit/JacobianSystem.hpp similarity index 100% rename from opm/core/transport/JacobianSystem.hpp rename to opm/core/transport/implicit/JacobianSystem.hpp diff --git a/opm/core/transport/NormSupport.hpp b/opm/core/transport/implicit/NormSupport.hpp similarity index 100% rename from opm/core/transport/NormSupport.hpp rename to opm/core/transport/implicit/NormSupport.hpp diff --git a/opm/core/transport/SimpleFluid2pWrappingProps.hpp b/opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp similarity index 96% rename from opm/core/transport/SimpleFluid2pWrappingProps.hpp rename to opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp index 880d8f66..92dbd297 100644 --- a/opm/core/transport/SimpleFluid2pWrappingProps.hpp +++ b/opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp @@ -66,5 +66,5 @@ namespace Opm{ }; } -#include +#include #endif // SIMPLEFLUID2PWRAPPINGPROPS_HPP diff --git a/opm/core/transport/SimpleFluid2pWrappingProps_impl.hpp b/opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp similarity index 97% rename from opm/core/transport/SimpleFluid2pWrappingProps_impl.hpp rename to opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp index c1aba4ad..09dc6731 100644 --- a/opm/core/transport/SimpleFluid2pWrappingProps_impl.hpp +++ b/opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp @@ -28,7 +28,7 @@ #ifndef SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP #define SIMPLEFLUID2PWRAPPINGPROPS_IMPL_HPP -#include +#include #include #include namespace Opm{ diff --git a/opm/core/transport/SinglePointUpwindTwoPhase.hpp b/opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp similarity index 100% rename from opm/core/transport/SinglePointUpwindTwoPhase.hpp rename to opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp diff --git a/opm/core/transport/TransportSolverTwophaseImplicit.cpp b/opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp similarity index 98% rename from opm/core/transport/TransportSolverTwophaseImplicit.cpp rename to opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp index 2b9b6fa3..7ca58932 100644 --- a/opm/core/transport/TransportSolverTwophaseImplicit.cpp +++ b/opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp @@ -27,7 +27,7 @@ */ -#include +#include #include #include diff --git a/opm/core/transport/TransportSolverTwophaseImplicit.hpp b/opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp similarity index 88% rename from opm/core/transport/TransportSolverTwophaseImplicit.hpp rename to opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp index 95017b58..c817f9ee 100644 --- a/opm/core/transport/TransportSolverTwophaseImplicit.hpp +++ b/opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp @@ -21,24 +21,22 @@ #ifndef OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED #define OPM_TRANSPORTSOLVERTWOPHASEIMPLICIT_HEADER_INCLUDED -#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include -#include -#include -#include #include #include -#include -#include -#include -#include -#include -#include -#include -#include - #include #include diff --git a/opm/core/transport/transport_source.c b/opm/core/transport/implicit/transport_source.c similarity index 98% rename from opm/core/transport/transport_source.c rename to opm/core/transport/implicit/transport_source.c index 5e8cdcdf..9fc56e14 100644 --- a/opm/core/transport/transport_source.c +++ b/opm/core/transport/implicit/transport_source.c @@ -36,7 +36,7 @@ #include #include -#include +#include /* ---------------------------------------------------------------------- */ diff --git a/opm/core/transport/transport_source.h b/opm/core/transport/implicit/transport_source.h similarity index 100% rename from opm/core/transport/transport_source.h rename to opm/core/transport/implicit/transport_source.h diff --git a/opm/core/transport/spu_explicit.c b/opm/core/transport/minimal/spu_explicit.c similarity index 97% rename from opm/core/transport/spu_explicit.c rename to opm/core/transport/minimal/spu_explicit.c index 7af40489..4512ace8 100644 --- a/opm/core/transport/spu_explicit.c +++ b/opm/core/transport/minimal/spu_explicit.c @@ -8,7 +8,7 @@ #include #include -#include +#include /* Twophase mobility-weighted upwind */ diff --git a/opm/core/transport/spu_explicit.h b/opm/core/transport/minimal/spu_explicit.h similarity index 100% rename from opm/core/transport/spu_explicit.h rename to opm/core/transport/minimal/spu_explicit.h diff --git a/opm/core/transport/spu_implicit.c b/opm/core/transport/minimal/spu_implicit.c similarity index 99% rename from opm/core/transport/spu_implicit.c rename to opm/core/transport/minimal/spu_implicit.c index d6fb1fd6..a5ad7c5d 100644 --- a/opm/core/transport/spu_implicit.c +++ b/opm/core/transport/minimal/spu_implicit.c @@ -11,7 +11,7 @@ #include #include -#include +#include diff --git a/opm/core/transport/spu_implicit.h b/opm/core/transport/minimal/spu_implicit.h similarity index 100% rename from opm/core/transport/spu_implicit.h rename to opm/core/transport/minimal/spu_implicit.h diff --git a/opm/core/transport/reorder/nlsolvers.c b/opm/core/transport/reorder/nlsolvers.c deleted file mode 100644 index ae372eb2..00000000 --- a/opm/core/transport/reorder/nlsolvers.c +++ /dev/null @@ -1,308 +0,0 @@ -/* -Copyright (C) 2012 (c) Jostein R. Natvig - -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 -#include -#include -#include -#include - -#include - -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: */ diff --git a/opm/core/transport/reorder/nlsolvers.h b/opm/core/transport/reorder/nlsolvers.h deleted file mode 100644 index 2f839d81..00000000 --- a/opm/core/transport/reorder/nlsolvers.h +++ /dev/null @@ -1,55 +0,0 @@ -/* -Copyright (C) 2012 (c) Jostein R. Natvig - -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: */