OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4497 from akva2/nonlinearsolverebos_compile_unit

Browse Source 
NonlinearsolverEbos: add compile unit
This commit is contained in:
Bård Skaflestad
2023-03-01 16:14:22 +01:00
committed by GitHub
parent 8c038fd36c cbaaca4f14
commit 00da723b59
3 changed files with 151 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@@ -47,6 +47,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/flow/FlowMainEbos.cpp
opm/simulators/flow/KeywordValidation.cpp
opm/simulators/flow/Main.cpp
opm/simulators/flow/NonlinearSolverEbos.cpp
opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.cpp
opm/simulators/flow/ValidationFunctions.cpp
opm/simulators/linalg/bda/WellContributions.cpp

120
opm/simulators/flow/NonlinearSolverEbos.cpp Normal file
View File

@@ -0,0 +1,120 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
Copyright 2015 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/flow/NonlinearSolverEbos.hpp>
#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <opm/common/ErrorMacros.hpp>
#include <cmath>
#include <stdexcept>
namespace Opm {
namespace detail {
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const int it, const int numPhases, const double relaxRelTol,
bool& oscillate, bool& stagnate)
{
// The detection of oscillation in two primary variable results in the report of the detection
// of oscillation for the solver.
// Only the saturations are used for oscillation detection for the black oil model.
// Stagnate is not used for any treatment here.
if (it < 2) {
oscillate = false;
stagnate = false;
return;
}
stagnate = true;
int oscillatePhase = 0;
const auto& F0 = residualHistory[it];
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]);
oscillatePhase += (d1 < relaxRelTol) && (relaxRelTol < d2);
// Process is 'stagnate' unless at least one phase
// exhibits significant residual change.
stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
}
oscillate = (oscillatePhase > 1);
}
template <class BVector>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
const double omega,
NonlinearRelaxType relaxType)
{
// The dxOld is updated with dx.
// If omega is equal to 1., no relaxtion will be appiled.
BVector tempDxOld = dxOld;
dxOld = dx;
switch (relaxType) {
case NonlinearRelaxType::Dampen: {
if (omega == 1.) {
return;
}
for (auto& d : dx) {
d *= omega;
}
return;
}
case NonlinearRelaxType::SOR: {
if (omega == 1.) {
return;
}
auto i = dx.size();
for (i = 0; i < dx.size(); ++i) {
dx[i] *= omega;
tempDxOld[i] *= (1.-omega);
dx[i] += tempDxOld[i];
}
return;
}
default:
OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
}
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)
} // namespace detail
} // namespace Opm

103
opm/simulators/flow/NonlinearSolverEbos.hpp
View File

@@ -82,6 +82,27 @@ namespace Opm {
class WellState;
// Available relaxation scheme types.
enum class NonlinearRelaxType {
Dampen,
SOR
};
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,
bool& oscillate, bool& stagnate);
/// 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, NonlinearRelaxType relaxType);
}
/// A nonlinear solver class suitable for general fully-implicit models,
/// as well as pressure, transport and sequential models.
template <class TypeTag, class PhysicalModel>
@@ -90,16 +111,10 @@ class WellState;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
// Available relaxation scheme types.
enum RelaxType {
Dampen,
SOR
};
// Solver parameters controlling nonlinear process.
struct SolverParameters
{
RelaxType relaxType_;
NonlinearRelaxType relaxType_;
double relaxMax_;
double relaxIncrement_;
double relaxRelTol_;
@@ -118,9 +133,9 @@ class WellState;
const auto& relaxationTypeString = EWOMS_GET_PARAM(TypeTag, std::string, NewtonRelaxationType);
if (relaxationTypeString == "dampen") {
relaxType_ = Dampen;
relaxType_ = NonlinearRelaxType::Dampen;
} else if (relaxationTypeString == "sor") {
relaxType_ = SOR;
relaxType_ = NonlinearRelaxType::SOR;
} else {
OPM_THROW(std::runtime_error,
"Unknown Relaxtion Type " + relaxationTypeString);
@@ -138,7 +153,7 @@ class WellState;
void reset()
{
// default values for the solver parameters
relaxType_ = Dampen;
relaxType_ = NonlinearRelaxType::Dampen;
relaxMax_ = 0.5;
relaxIncrement_ = 0.1;
relaxRelTol_ = 0.2;
@@ -279,34 +294,8 @@ class WellState;
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const int it, bool& oscillate, bool& stagnate) const
{
// The detection of oscillation in two primary variable results in the report of the detection
// of oscillation for the solver.
// Only the saturations are used for oscillation detection for the black oil model.
// Stagnate is not used for any treatment here.
if ( it < 2 ) {
oscillate = false;
stagnate = false;
return;
}
stagnate = true;
int oscillatePhase = 0;
const std::vector<double>& F0 = residualHistory[it];
const std::vector<double>& F1 = residualHistory[it - 1];
const std::vector<double>& F2 = residualHistory[it - 2];
for (int p= 0; p < model_->numPhases(); ++p){
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2);
// Process is 'stagnate' unless at least one phase
// exhibits significant residual change.
stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
}
oscillate = (oscillatePhase > 1);
detail::detectOscillations(residualHistory, it, model_->numPhases(),
this->relaxRelTol(), oscillate, stagnate);
}
@@ -315,40 +304,7 @@ class WellState;
template <class BVector>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
{
// The dxOld is updated with dx.
// If omega is equal to 1., no relaxtion will be appiled.
BVector tempDxOld = dxOld;
dxOld = dx;
switch (relaxType()) {
case Dampen: {
if (omega == 1.) {
return;
}
auto i = dx.size();
for (i = 0; i < dx.size(); ++i) {
dx[i] *= omega;
}
return;
}
case SOR: {
if (omega == 1.) {
return;
}
auto i = dx.size();
for (i = 0; i < dx.size(); ++i) {
dx[i] *= omega;
tempDxOld[i] *= (1.-omega);
dx[i] += tempDxOld[i];
}
return;
}
default:
OPM_THROW(std::runtime_error, "Can only handle Dampen and SOR relaxation type.");
}
return;
detail::stabilizeNonlinearUpdate(dx, dxOld, omega, this->relaxType());
}
/// The greatest relaxation factor (i.e. smallest factor) allowed.
@@ -360,7 +316,7 @@ class WellState;
{ return param_.relaxIncrement_; }
/// The relaxation type (Dampen or SOR).
enum RelaxType relaxType() const
NonlinearRelaxType relaxType() const
{ return param_.relaxType_; }
/// The relaxation relative tolerance.
@@ -392,6 +348,7 @@ class WellState;
int linearIterationsLast_;
int wellIterationsLast_;
};
} // namespace Opm
#endif // OPM_NON_LINEAR_SOLVER_EBOS_HPP