From a23bfec105007c70ba10903d5f656f3e34ee18ab Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 1 Mar 2012 13:25:59 +0100
Subject: [PATCH] Allow for fpi benchmarking
---
src/one-body-sample.cc | 59 ++++++++++++++++++++++++++++++++++++++
src/one-body-sample.parset | 4 +++
2 files changed, 63 insertions(+)
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 637549e9..e4974849 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -323,19 +323,26 @@ int main(int argc, char *argv[]) {
VectorType u2 = u1;
VectorType u3 = u1;
VectorType u4 = u1;
+ VectorType u5 = u1;
SingletonVectorType s4_old(grid->size(grid->maxLevel(), dim));
+ SingletonVectorType s5_old(grid->size(grid->maxLevel(), dim));
s4_old = parset.get<double>("boundary.friction.state.initial");
+ s5_old = s4_old;
VectorType u1_diff(grid->size(grid->maxLevel(), dim));
u1_diff = 0.0; // Has to be zero!
VectorType u2_diff = u1_diff;
VectorType u3_diff = u1_diff;
VectorType u4_diff = u1_diff;
+ VectorType u5_diff = u1_diff;
auto s4_new = Dune::make_shared<SingletonVectorType>(
grid->size(grid->maxLevel(), dim));
+ auto s5_new = Dune::make_shared<SingletonVectorType>(
+ grid->size(grid->maxLevel(), dim));
*s4_new = s4_old;
+ *s5_new = *s4_new;
SingletonVectorType vonMisesStress;
@@ -343,6 +350,7 @@ int main(int argc, char *argv[]) {
VectorType b2;
VectorType b3;
VectorType b4;
+ VectorType b5;
auto const nodalIntegrals =
assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
@@ -407,12 +415,14 @@ int main(int argc, char *argv[]) {
b2 = b1;
b3 = b1;
b4 = b1;
+ b5 = b1;
// b -= linear update
stiffnessMatrix.mmv(u1, b1);
stiffnessMatrix.mmv(u2, b2);
stiffnessMatrix.mmv(u3, b3);
stiffnessMatrix.mmv(u4, b4);
+ stiffnessMatrix.mmv(u5, b5);
if (parset.get<bool>("solver.nonlineargs.use")) {
auto state =
@@ -501,6 +511,55 @@ int main(int argc, char *argv[]) {
u4 += u4_diff;
s4_old = *s4_new;
+ if (parset.get<bool>("benchmarks.fpi.enable")) {
+ for (int state_fpi = 0;
+ state_fpi < parset.get<int>("benchmarks.fpi.iterations");
+ ++state_fpi) {
+ auto myGlobalNonlinearity =
+ assemble_nonlinearity<VectorType, OperatorType>(
+ grid->size(grid->maxLevel(), dim), parset, nodalIntegrals,
+ s5_new, h);
+ MyConvexProblemType const myConvexProblem(stiffnessMatrix,
+ *myGlobalNonlinearity, b5);
+
+ MyBlockProblemType myBlockProblem(parset, myConvexProblem);
+ multigridStep.setProblem(u5_diff, myBlockProblem);
+
+ LoopSolver<VectorType> overallSolver(
+ &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
+ solver_tolerance, &energyNorm, verbosity);
+ overallSolver.solve();
+
+ if (!parset.get<bool>("boundary.friction.state.evolve"))
+ continue;
+
+ for (size_t i = 0; i < frictionalNodes.size(); ++i) {
+ if (frictionalNodes[i][0]) {
+ double const L = parset.get<double>("boundary.friction.ruina.L");
+ double const unorm = u5_diff[i].two_norm();
+
+ double ret1 =
+ compute_state_update_bisection(h, unorm, L, s5_old[i][0]);
+ assert(abs(1.0 / h * ret1 - (s5_old[i] - unorm / L) / h -
+ exp(-ret1)) < 1e-12);
+ double ret2 =
+ compute_state_update_lambert(h, unorm, L, s5_old[i][0]);
+ assert(abs(1.0 / h * ret2 - (s5_old[i] - unorm / L) / h -
+ exp(-ret2)) < 1e-12);
+ assert(abs(ret2 - ret1) < 1e-14);
+
+ (*s5_new)[i][0] = ret1;
+ }
+ }
+ }
+ }
+
+ u5 += u5_diff;
+ s5_old = *s5_new;
+
+ if (parset.get<bool>("benchmarks.fpi.enable"))
+ std::cout << energyNorm.diff(u5, u4) / energyNorm(u5) << std::endl;
+
// Compute von Mises stress and write everything to a file
if (parset.get<bool>("writeVTK")) {
auto const displacement =
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index 317c1792..d4baea4c 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -9,6 +9,10 @@ printDifference = false
writeVTK = true
+[benchmarks.fpi]
+enable = false
+iterations = 50
+
[grid]
refinements = 7
--
GitLab