From d01966c780bf151f43cdd35272a61317b2488077 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 14 Jan 2011 16:29:28 +0000
Subject: [PATCH] Rewrote the convergence monitoring using only the interface
values.
After all, we are thinking 'Steklov-Poincare' here, which means
we are solving an equation on the interface. The old code
taking into account the subdomain solutions wasn't really satisfying:
in particular, it was not symmetric wrt the subdomains.
[[Imported from SVN: r6762]]
---
dirneucoupling.cc | 53 ++++++++---------------------------------------
1 file changed, 9 insertions(+), 44 deletions(-)
diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index 0137963c..740b029a 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -344,7 +344,7 @@ int main (int argc, char *argv[]) try
FieldVector<double,3> lambdaTorque(0);
//
- double normOfOldCorrection = 0;
+ double normOfOldCorrection = 1;
int dnStepsActuallyTaken = 0;
for (int i=0; i<maxDDIterations; i++) {
@@ -352,9 +352,8 @@ int main (int argc, char *argv[]) try
std::cout << " Domain-Decomposition- Step Number: " << i << std::endl;
std::cout << "----------------------------------------------------" << std::endl;
- // Backup of the current solution for the error computation later on
- VectorType oldSolution3d = x3d;
- RodSolutionType oldSolutionRod = rodX;
+ // Backup of the current iterate for the error computation later on
+ RigidBodyMotion<3> oldLambda = lambda;
if (ddType=="FixedPointIteration") {
@@ -516,54 +515,20 @@ int main (int argc, char *argv[]) try
// Compute error in the energy norm
// ////////////////////////////////////////////
- // the 3d body
- double oldNorm = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, stiffnessMatrix3d);
- oldSolution3d -= x3d;
- double normOfCorrection = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, stiffnessMatrix3d);
+ double lengthOfCorrection = RigidBodyMotion<3>::distance(oldLambda, lambda);
- double max3dRelCorrection = 0;
- for (size_t j=0; j<x3d.size(); j++)
- for (int k=0; k<dim; k++)
- max3dRelCorrection = std::max(max3dRelCorrection,
- std::fabs(oldSolution3d[j][k])/ std::fabs(x3d[j][k]));
-
- // the rod
- RodDifferenceType rodDiff = computeGeodesicDifference(oldSolutionRod, rodX);
- double maxRodRelCorrection = 0;
- for (size_t j=0; j<rodX.size(); j++)
- for (int k=0; k<dim; k++)
- maxRodRelCorrection = std::max(maxRodRelCorrection,
- std::fabs(rodDiff[j][k])/ std::fabs(rodX[j].r[k]));
-
- // Absolute corrections
- double maxRodCorrection = computeGeodesicDifference(oldSolutionRod, rodX).infinity_norm();
- double max3dCorrection = oldSolution3d.infinity_norm();
-
-
- std::cout << "rod correction: " << maxRodCorrection
- << " rod rel correction: " << maxRodRelCorrection
- << " 3d correction: " << max3dCorrection
- << " 3d rel correction: " << max3dRelCorrection << std::endl;
-
-
- oldNorm = std::sqrt(oldNorm);
-
- normOfCorrection = std::sqrt(normOfCorrection);
-
- double relativeError = normOfCorrection / oldNorm;
-
- double convRate = normOfCorrection / normOfOldCorrection;
+ double convRate = lengthOfCorrection / normOfOldCorrection;
- normOfOldCorrection = normOfCorrection;
+ normOfOldCorrection = lengthOfCorrection;
// Output
- std::cout << "DD iteration: " << i << " -- ||u^{n+1} - u^n|| / ||u^n||: " << relativeError << ", "
+ std::cout << "DD iteration: " << i << ", "
+ << "interface correction: " << lengthOfCorrection << ", "
<< "convrate " << convRate << "\n";
dnStepsActuallyTaken = i;
- //if (relativeError < ddTolerance)
- if (std::max(max3dRelCorrection,maxRodRelCorrection) < ddTolerance)
+ if (lengthOfCorrection < ddTolerance)
break;
}
--
GitLab