NonLinearSolver: use Scalar type

This commit is contained in:
Arne Morten Kvarving 2024-02-22 14:44:48 +01:00
parent cdf227bcbd
commit 39554ab7dd
2 changed files with 40 additions and 29 deletions

View File

@ -31,8 +31,9 @@
namespace Opm {
namespace detail {
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const int it, const int numPhases, const double relaxRelTol,
template<class Scalar>
void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
const int it, const int numPhases, const Scalar relaxRelTol,
bool& oscillate, bool& stagnate)
{
// The detection of oscillation in two primary variable results in the report of the detection
@ -52,8 +53,8 @@ void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const auto& F1 = residualHistory[it - 1];
const auto& F2 = residualHistory[it - 2];
for (int p = 0; p < numPhases; ++p) {
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
const Scalar d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const Scalar d2 = std::abs((F0[p] - F1[p]) / F0[p]);
oscillatePhase += (d1 < relaxRelTol) && (relaxRelTol < d2);
@ -65,9 +66,9 @@ void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
oscillate = (oscillatePhase > 1);
}
template <class BVector>
template <class BVector, class Scalar>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
const double omega,
const Scalar omega,
NonlinearRelaxType relaxType)
{
// The dxOld is updated with dx.
@ -104,16 +105,25 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
return;
}
template<int Size> using BV = Dune::BlockVector<Dune::FieldVector<double,Size>>;
#define INSTANCE(Size) \
template void stabilizeNonlinearUpdate<BV<Size>>(BV<Size>&, BV<Size>&, \
const double, NonlinearRelaxType);
INSTANCE(1)
INSTANCE(2)
INSTANCE(3)
INSTANCE(4)
INSTANCE(5)
INSTANCE(6)
template<class Scalar, int Size>
using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Size>>;
#define INSTANCE(T,Size) \
template void stabilizeNonlinearUpdate<BV<T,Size>,T>(BV<T,Size>&, BV<T,Size>&, \
const T, NonlinearRelaxType);
#define INSTANCE_TYPE(T) \
template void detectOscillations(const std::vector<std::vector<T>>&, \
const int, const int, const T, \
bool&, bool&); \
INSTANCE(T,1) \
INSTANCE(T,2) \
INSTANCE(T,3) \
INSTANCE(T,4) \
INSTANCE(T,5) \
INSTANCE(T,6)
INSTANCE_TYPE(double)
} // namespace detail
} // namespace Opm

View File

@ -89,15 +89,16 @@ enum class NonlinearRelaxType {
namespace detail {
/// Detect oscillation or stagnation in a given residual history.
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const int it, const int numPhases, const double relaxRelTol,
template<class Scalar>
void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
const int it, const int numPhases, const Scalar relaxRelTol,
bool& oscillate, bool& stagnate);
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
/// Implemention for Dune block vectors.
template <class BVector>
template <class BVector, class Scalar>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
const double omega, NonlinearRelaxType relaxType);
const Scalar omega, NonlinearRelaxType relaxType);
}
@ -113,9 +114,9 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
struct SolverParameters
{
NonlinearRelaxType relaxType_;
double relaxMax_;
double relaxIncrement_;
double relaxRelTol_;
Scalar relaxMax_;
Scalar relaxIncrement_;
Scalar relaxRelTol_;
int maxIter_; // max nonlinear iterations
int minIter_; // min nonlinear iterations
@ -277,7 +278,7 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
int wellIterationsLastStep() const
{ return wellIterationsLast_; }
std::vector<std::vector<double> >
std::vector<std::vector<Scalar> >
computeFluidInPlace(const std::vector<int>& fipnum) const
{ return model_->computeFluidInPlace(fipnum); }
@ -290,7 +291,7 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
{ return *model_; }
/// Detect oscillation or stagnation in a given residual history.
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
const int it, bool& oscillate, bool& stagnate) const
{
detail::detectOscillations(residualHistory, it, model_->numPhases(),
@ -301,17 +302,17 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
/// Implemention for Dune block vectors.
template <class BVector>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const Scalar omega) const
{
detail::stabilizeNonlinearUpdate(dx, dxOld, omega, this->relaxType());
}
/// The greatest relaxation factor (i.e. smallest factor) allowed.
double relaxMax() const
Scalar relaxMax() const
{ return param_.relaxMax_; }
/// The step-change size for the relaxation factor.
double relaxIncrement() const
Scalar relaxIncrement() const
{ return param_.relaxIncrement_; }
/// The relaxation type (Dampen or SOR).
@ -319,7 +320,7 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
{ return param_.relaxType_; }
/// The relaxation relative tolerance.
double relaxRelTol() const
Scalar relaxRelTol() const
{ return param_.relaxRelTol_; }
/// The maximum number of nonlinear iterations allowed.