Loading...
Searching...
No Matches
ImplicitCapillarity_impl.hpp
//===========================================================================
//
// File: ImplicitCapillarity_impl.hpp
//
// Created: Thu May 6 15:36:07 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Jostein R Natvig <jostein.r.natvig@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
Copyright 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS 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.
OpenRS 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 OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
#define OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/Average.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/grid/common/Volumes.hpp>
#include <opm/porsol/common/ReservoirPropertyFixedMobility.hpp>
#include <opm/porsol/euler/MatchSaturatedVolumeFunctor.hpp>
#include <cassert>
#include <cmath>
#include <algorithm>
#include <iostream>
namespace Opm
{
template <class GI, class RP, class BC, template <class, class> class IP>
inline ImplicitCapillarity<GI, RP, BC, IP>::ImplicitCapillarity()
: method_viscous_(true),
method_gravity_(true),
check_sat_(true),
clamp_sat_(false),
residual_tolerance_(1e-8),
linsolver_verbosity_(1),
linsolver_type_(1),
update_relaxation_(1.0)
{
}
template <class GI, class RP, class BC, template <class, class> class IP>
inline ImplicitCapillarity<GI, RP, BC, IP>::ImplicitCapillarity(const GI& g, const RP& r, const BC& b)
: residual_(g, r, b),
method_viscous_(true),
method_gravity_(true),
check_sat_(true),
clamp_sat_(false),
residual_tolerance_(1e-8),
linsolver_verbosity_(1),
linsolver_type_(1),
update_relaxation_(1.0)
{
Dune::FieldVector<double, GI::Dimension> grav(0.0);
psolver_.init(g, r, grav, b);
}
template <class GI, class RP, class BC, template <class, class> class IP>
inline void ImplicitCapillarity<GI, RP, BC, IP>::init(const Opm::ParameterGroup& param)
{
method_viscous_ = param.getDefault("method_viscous", method_viscous_);
method_gravity_ = param.getDefault("method_gravity", method_gravity_);
check_sat_ = param.getDefault("check_sat", check_sat_);
clamp_sat_ = param.getDefault("clamp_sat", clamp_sat_);
residual_tolerance_ = param.getDefault("residual_tolerance", residual_tolerance_);
linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_);
linsolver_type_ = param.getDefault("linsolver_type", linsolver_type_);
update_relaxation_ = param.getDefault("update_relaxation", update_relaxation_);
}
template <class GI, class RP, class BC, template <class, class> class IP>
inline void ImplicitCapillarity<GI, RP, BC, IP>::init(const Opm::ParameterGroup& param,
const GI& g, const RP& r, const BC& b)
{
init(param);
initObj(g, r, b);
}
template <class GI, class RP, class BC, template <class, class> class IP>
inline void ImplicitCapillarity<GI, RP, BC, IP>::initObj(const GI& g, const RP& r, const BC& b)
{
residual_.initObj(g, r, b);
Dune::FieldVector<double, GI::Dimension> grav(0.0);
psolver_.init(g, r, grav, b);
}
namespace ImplicitCapillarityDetails {
void thresholdMobility(double& m, double threshold);
// The matrix variant expects diagonal mobilities.
template <class SomeMatrixType>
void thresholdMobility(SomeMatrixType& m, double threshold)
{
for (int i = 0; i < std::min(m.numRows(), m.numCols()); ++i) {
m(i, i) = std::max(m(i, i), threshold);
}
}
} // anon namespace
template <class GI, class RP, class BC, template <class, class> class IP>
template <class PressureSolution>
void ImplicitCapillarity<GI, RP, BC, IP>::transportSolve(std::vector<double>& saturation,
const double /*time*/,
const typename GI::Vector& gravity,
const PressureSolution& pressure_sol,
const Opm::SparseVector<double>& injection_rates) const
{
// Start a timer.
Opm::time::StopWatch clock;
clock.start();
// Compute capillary mobilities.
typedef typename RP::Mobility Mob;
int num_cells = saturation.size();
std::vector<Mob> cap_mob(num_cells);
for (int c = 0; c < num_cells; ++c) {
Mob& m = cap_mob[c];
residual_.reservoirProperties().phaseMobility(0, c, saturation[c], m.mob);
Mob mob2;
residual_.reservoirProperties().phaseMobility(1, c, saturation[c], mob2.mob);
Mob mob_tot;
mob_tot.setToSum(m, mob2);
Mob mob_totinv;
mob_totinv.setToInverse(mob_tot);
m *= mob_totinv;
m *= mob2;
ImplicitCapillarityDetails::thresholdMobility(m.mob, 1e-10); // @@TODO: User-set limit.
// std::cout << m.mob(0,0) << '\n';
}
ReservoirPropertyFixedMobility<Mob> capillary_mobilities(cap_mob);
// Set up boundary conditions.
BC cap_press_bcs(residual_.boundaryConditions());
for (int i = 0; i < cap_press_bcs.size(); ++i) {
if (cap_press_bcs.flowCond(i).isPeriodic()) {
cap_press_bcs.flowCond(i) = FlowBC(FlowBC::Periodic, 0.0);
}
}
// Compute injection rates from residual.
std::vector<double> injection_rates_residual(num_cells);
residual_.computeResidual(saturation, gravity, pressure_sol, injection_rates,
method_viscous_, method_gravity_, false,
injection_rates_residual);
for (int i = 0; i < num_cells; ++i) {
injection_rates_residual[i] = -injection_rates_residual[i];
}
// Compute capillary pressure.
// Note that the saturation is just a dummy for this call, since the mobilities are fixed.
psolver_.solve(capillary_mobilities, saturation, cap_press_bcs, injection_rates_residual,
residual_tolerance_, linsolver_verbosity_, linsolver_type_);
// Solve for constant to change capillary pressure solution by.
std::vector<double> cap_press(num_cells);
const PressureSolution& pcapsol = psolver_.getSolution();
for (CIt c = residual_.grid().cellbegin(); c != residual_.grid().cellend(); ++c) {
cap_press[c->index()] = pcapsol.pressure(c);
}
MatchSaturatedVolumeFunctor<GI, RP> functor(residual_.grid(),
residual_.reservoirProperties(),
saturation,
cap_press);
double min_cap_press = *std::min_element(cap_press.begin(), cap_press.end());
double max_cap_press = *std::max_element(cap_press.begin(), cap_press.end());
double cap_press_range = max_cap_press - min_cap_press;
double mod_low = 1e100;
double mod_high = -1e100;
Opm::bracketZero(functor, 0.0, cap_press_range, mod_low, mod_high);
const int max_iter = 40;
const double nonlinear_tolerance = 1e-12;
int iterations_used = -1;
typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
double mod_correct = RootFinder::solve(functor, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
std::cout << "Moved capillary pressure solution by " << mod_correct << " after "
<< iterations_used << " iterations." << std::endl;
// saturation = functor.lastSaturations();
const std::vector<double>& sat_new = functor.lastSaturations();
for (int i = 0; i < num_cells; ++i) {
saturation[i] = (1.0 - update_relaxation_)*saturation[i] + update_relaxation_*sat_new[i];
}
// Optionally check and/or clamp results.
if (check_sat_ || clamp_sat_) {
checkAndPossiblyClampSat(saturation);
}
// Stop timer and optionally print seconds taken.
clock.stop();
#ifdef VERBOSE
std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
#endif // VERBOSE
}
template <class GI, class RP, class BC, template <class, class> class IP>
inline void ImplicitCapillarity<GI, RP, BC, IP>::checkAndPossiblyClampSat(std::vector<double>& s) const
{
int num_cells = s.size();
for (int cell = 0; cell < num_cells; ++cell) {
if (s[cell] > 1.0 || s[cell] < 0.0) {
if (clamp_sat_) {
s[cell] = std::max(std::min(s[cell], 1.0), 0.0);
} else if (s[cell] > 1.001 || s[cell] < -0.001) {
OPM_THROW(std::runtime_error,
"Saturation out of range in ImplicitCapillarity: "
"Cell " + std::to_string(cell) +
" sat " + std::to_string(s[cell]));
}
}
}
}
} // end namespace Opm
#endif // OPENRS_IMPLICITCAPILLARITY_IMPL_HEADER
Opm::FlowBC
A class for representing a flow boundary condition.
Definition BoundaryConditions.hpp:122
Opm::ImplicitCapillarity::initObj
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary)
Definition ImplicitCapillarity_impl.hpp:111
Opm::ImplicitCapillarity::ImplicitCapillarity
ImplicitCapillarity()
Definition ImplicitCapillarity_impl.hpp:58
Opm::ImplicitCapillarity::init
void init(const Opm::ParameterGroup &param)
Definition ImplicitCapillarity_impl.hpp:89
Opm::ImplicitCapillarity::transportSolve
void transportSolve(std::vector< double > &saturation, const double time, const typename GridInterface::Vector &gravity, const PressureSolution &pressure_sol, const Opm::SparseVector< double > &injection_rates) const
Solve transport equation.
Opm::IncompFlowSolverHybrid::init
void init(const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
All-in-one initialization routine.
Definition IncompFlowSolverHybrid.hpp:476
Opm::ReservoirPropertyFixedMobility
Definition ReservoirPropertyFixedMobility.hpp:46
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::MatchSaturatedVolumeFunctor
Definition MatchSaturatedVolumeFunctor.hpp:68