130 PrimaryVariables& nextValue,
131 const PrimaryVariables& currentValue,
132 const EqVector& update,
136 nextValue = currentValue;
144 Scalar sumSatDelta = 0.0;
145 Scalar maxSatDelta = 0.0;
146 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
147 maxSatDelta = std::max(std::abs(update[saturation0Idx + phaseIdx]),
149 sumSatDelta += update[saturation0Idx + phaseIdx];
151 maxSatDelta = std::max(std::abs(- sumSatDelta), maxSatDelta);
153 if (maxSatDelta > 0.2) {
154 Scalar alpha = 0.2/maxSatDelta;
155 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
156 nextValue[saturation0Idx + phaseIdx] =
157 currentValue[saturation0Idx + phaseIdx]
158 - alpha*update[saturation0Idx + phaseIdx];
163 clampValue_(nextValue[pressure0Idx],
164 currentValue[pressure0Idx]*0.8,
165 currentValue[pressure0Idx]*1.2);
168 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
169 Scalar& val = nextValue[fugacity0Idx + compIdx];
170 Scalar oldVal = currentValue[fugacity0Idx + compIdx];
174 Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
176 minPhi = std::max(0.001*currentValue[pressure0Idx], minPhi);
180 Scalar maxDelta = 0.7 * minPhi;
181 clampValue_(val, oldVal - maxDelta, oldVal + maxDelta);
184 val = std::max(val, 0.0);
189 if (this->numIterations_ < 3) {
191 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
192 Scalar& val = nextValue[fugacity0Idx + compIdx];
193 Scalar oldVal = currentValue[fugacity0Idx + compIdx];
194 Scalar minPhi = this->problem().model().minActivityCoeff(globalDofIdx, compIdx);
195 if (oldVal < 1.0*minPhi && val > 1.0*minPhi)
197 else if (oldVal > 0.0 && val < 0.0)
202 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
203 Scalar& val = nextValue[saturation0Idx + phaseIdx];
204 Scalar oldVal = currentValue[saturation0Idx + phaseIdx];
205 if (oldVal < 1.0 && val > 1.0)
207 else if (oldVal > 0.0 && val < 0.0)