From 434513bf85950f33063423cdc576bcb5c00e462a Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 6 Jul 2014 11:59:11 +0200
Subject: [PATCH] [Cleanup] Bug fix: Do not mix old and new BC

During the first fixed point iteration, velocityIterate and the old
velocity are identical except for the new BC: velocityIterate has
already had them imposed.

To simplify the control flow, kill v completely.

Note that since we always make two fixed point iterations and this issue
only occurs during the first step, this change only affects
computational results as (a reduction of) noise
---
 src/sand-wedge.cc | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index a4587a58..5782bb30 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -352,7 +352,6 @@ int main(int argc, char *argv[]) {
         parset.get<double>("boundary.friction.L"),
         parset.get<double>("boundary.friction.V0"));
 
-    Vector v = v_initial;
     Vector v_m(leafVertexCount);
     ScalarVector alpha(leafVertexCount);
 
@@ -413,7 +412,7 @@ int main(int argc, char *argv[]) {
       for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
         timeSteppingScheme->extractOldVelocity(v_m);
         v_m *= 1.0 - lambda;
-        Arithmetic::addProduct(v_m, lambda, v);
+        Arithmetic::addProduct(v_m, lambda, velocityIterate);
 
         stateUpdater->solve(v_m);
         stateUpdater->extractAlpha(alpha);
@@ -421,21 +420,21 @@ int main(int argc, char *argv[]) {
         solveVelocityProblem(velocityIterate, alpha);
         timeSteppingScheme->postProcess(velocityIterate);
         timeSteppingScheme->extractDisplacement(u);
-        timeSteppingScheme->extractVelocity(v);
 
         iterationWriter << iterationCounter << " ";
         if (printProgress)
           std::cout << '.' << std::flush;
 
         if (stateFPI > 1) {
-          double const velocityCorrection = AMNorm.diff(v_saved, v);
+          double const velocityCorrection =
+              AMNorm.diff(v_saved, velocityIterate);
           if (velocityCorrection < fixedPointTolerance)
             break;
         }
         if (stateFPI == maximumStateFPI)
           DUNE_THROW(Dune::Exception, "FPI failed to converge");
 
-        v_saved = v;
+        v_saved = velocityIterate;
       }
       if (printProgress)
         std::cout << std::endl;
-- 
GitLab