Skip to content
Snippets Groups Projects
Commit b84ea96c authored by Elias Pipping's avatar Elias Pipping
Browse files

[Cleanup] Rewrite state update

parent ccfb5aca
No related branches found
No related tags found
No related merge requests found
...@@ -531,22 +531,24 @@ int main(int argc, char *argv[]) { ...@@ -531,22 +531,24 @@ int main(int argc, char *argv[]) {
// Q: is this reasonable? // Q: is this reasonable?
VectorType u; VectorType u;
VectorType u_saved; VectorType u_saved;
SingletonVectorType alpha_saved;
double lastStateCorrection; double lastStateCorrection;
for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) { for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
stateUpdater->solve(v); stateUpdater->solve(v);
{ stateUpdater->extractState(alpha);
SingletonVectorType computed_state;
stateUpdater->extractState(computed_state);
double const stateCorrection =
stateEnergyNorm.diff(computed_state, alpha);
if (state_fpi <= 2 or stateCorrection < if (state_fpi == 1)
requiredReduction * lastStateCorrection) { relaxationWriter << "N ";
alpha = computed_state; else {
double const stateCorrection =
stateEnergyNorm.diff(alpha, alpha_saved);
if (state_fpi <=
2 // lastStateCorrection is only set for state_fpi > 2
or stateCorrection < requiredReduction * lastStateCorrection)
relaxationWriter << "N "; relaxationWriter << "N ";
} else { else {
alpha *= relaxation; alpha *= (1.0 - relaxation);
Arithmetic::addProduct(alpha, 1.0 - relaxation, computed_state); Arithmetic::addProduct(alpha, relaxation, alpha_saved);
relaxationWriter << "Y "; relaxationWriter << "Y ";
} }
lastStateCorrection = stateCorrection; lastStateCorrection = stateCorrection;
...@@ -569,6 +571,7 @@ int main(int argc, char *argv[]) { ...@@ -569,6 +571,7 @@ int main(int argc, char *argv[]) {
if (state_fpi == state_fpi_max) if (state_fpi == state_fpi_max)
DUNE_THROW(Dune::Exception, "FPI failed to converge"); DUNE_THROW(Dune::Exception, "FPI failed to converge");
alpha_saved = alpha;
u_saved = u; u_saved = u;
} }
if (printProgress) if (printProgress)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment