26 #ifndef EWOMS_NCP_NEWTON_METHOD_HH
27 #define EWOMS_NCP_NEWTON_METHOD_HH
40 template <
class TypeTag>
43 typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
45 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
46 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
47 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
48 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
54 enum { fugacity0Idx = Indices::fugacity0Idx };
55 enum { saturation0Idx = Indices::saturation0Idx };
56 enum { pressure0Idx = Indices::pressure0Idx };
64 choppedIterations_ =
EWOMS_GET_PARAM(TypeTag,
int, NcpNewtonNumChoppedIterations);
65 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
73 ParentType::registerParameters();
76 "The number of Newton iterations for which the "
77 "update gets limited");
91 PrimaryVariables& nextValue,
92 const PrimaryVariables& currentValue,
93 const EqVector& update,
94 const EqVector& currentResidual)
97 nextValue = currentValue;
102 if (this->numIterations_ < choppedIterations_) {
103 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
104 saturationChop_(nextValue[saturation0Idx + phaseIdx],
105 currentValue[saturation0Idx + phaseIdx]);
106 pressureChop_(nextValue[pressure0Idx],
107 currentValue[pressure0Idx]);
110 for (
int compIdx = 0; compIdx < numComponents; ++compIdx) {
111 Scalar &val = nextValue[fugacity0Idx + compIdx];
112 Scalar oldVal = currentValue[fugacity0Idx + compIdx];
117 Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
118 Scalar maxDelta = 0.7 * minPhi;
120 clampValue_(val, oldVal - maxDelta, oldVal + maxDelta);
124 val = std::max(-0.01 * minPhi, val);
125 val = std::min(1.01 * minPhi, val);
131 void clampValue_(Scalar &val, Scalar minVal, Scalar maxVal)
const
132 { val = std::max(minVal, std::min(val, maxVal)); }
134 void pressureChop_(Scalar &val, Scalar oldVal)
const
137 clampValue_(val, oldVal * 0.8, oldVal * 1.2);
140 void saturationChop_(Scalar &val, Scalar oldVal)
const
143 const Scalar maxDelta = 0.20;
144 clampValue_(val, oldVal - maxDelta, oldVal + maxDelta);
147 int choppedIterations_;
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
A Newton solver specific to the NCP model.
Definition: ncpnewtonmethod.hh:41
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:73
NcpNewtonMethod(Simulator &simulator)
Definition: ncpnewtonmethod.hh:62
Definition: baseauxiliarymodule.hh:35
static void registerParameters()
Register all run-time parameters for the Newton method.
Definition: ncpnewtonmethod.hh:71
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:64
Declares the properties required for the NCP compositional multi-phase model.
void updatePrimaryVariables_(int globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector ¤tResidual)
Update a single primary variables object.
Definition: ncpnewtonmethod.hh:90
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:95