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

Allow for fpi benchmarking

parent 8d51689d
No related branches found
No related tags found
No related merge requests found
...@@ -323,19 +323,26 @@ int main(int argc, char *argv[]) { ...@@ -323,19 +323,26 @@ int main(int argc, char *argv[]) {
VectorType u2 = u1; VectorType u2 = u1;
VectorType u3 = u1; VectorType u3 = u1;
VectorType u4 = u1; VectorType u4 = u1;
VectorType u5 = u1;
SingletonVectorType s4_old(grid->size(grid->maxLevel(), dim)); 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"); s4_old = parset.get<double>("boundary.friction.state.initial");
s5_old = s4_old;
VectorType u1_diff(grid->size(grid->maxLevel(), dim)); VectorType u1_diff(grid->size(grid->maxLevel(), dim));
u1_diff = 0.0; // Has to be zero! u1_diff = 0.0; // Has to be zero!
VectorType u2_diff = u1_diff; VectorType u2_diff = u1_diff;
VectorType u3_diff = u1_diff; VectorType u3_diff = u1_diff;
VectorType u4_diff = u1_diff; VectorType u4_diff = u1_diff;
VectorType u5_diff = u1_diff;
auto s4_new = Dune::make_shared<SingletonVectorType>( auto s4_new = Dune::make_shared<SingletonVectorType>(
grid->size(grid->maxLevel(), dim)); grid->size(grid->maxLevel(), dim));
auto s5_new = Dune::make_shared<SingletonVectorType>(
grid->size(grid->maxLevel(), dim));
*s4_new = s4_old; *s4_new = s4_old;
*s5_new = *s4_new;
SingletonVectorType vonMisesStress; SingletonVectorType vonMisesStress;
...@@ -343,6 +350,7 @@ int main(int argc, char *argv[]) { ...@@ -343,6 +350,7 @@ int main(int argc, char *argv[]) {
VectorType b2; VectorType b2;
VectorType b3; VectorType b3;
VectorType b4; VectorType b4;
VectorType b5;
auto const nodalIntegrals = auto const nodalIntegrals =
assemble_frictional<GridType, GridView, SmallVector, P1Basis>( assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
...@@ -407,12 +415,14 @@ int main(int argc, char *argv[]) { ...@@ -407,12 +415,14 @@ int main(int argc, char *argv[]) {
b2 = b1; b2 = b1;
b3 = b1; b3 = b1;
b4 = b1; b4 = b1;
b5 = b1;
// b -= linear update // b -= linear update
stiffnessMatrix.mmv(u1, b1); stiffnessMatrix.mmv(u1, b1);
stiffnessMatrix.mmv(u2, b2); stiffnessMatrix.mmv(u2, b2);
stiffnessMatrix.mmv(u3, b3); stiffnessMatrix.mmv(u3, b3);
stiffnessMatrix.mmv(u4, b4); stiffnessMatrix.mmv(u4, b4);
stiffnessMatrix.mmv(u5, b5);
if (parset.get<bool>("solver.nonlineargs.use")) { if (parset.get<bool>("solver.nonlineargs.use")) {
auto state = auto state =
...@@ -501,6 +511,55 @@ int main(int argc, char *argv[]) { ...@@ -501,6 +511,55 @@ int main(int argc, char *argv[]) {
u4 += u4_diff; u4 += u4_diff;
s4_old = *s4_new; 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 // Compute von Mises stress and write everything to a file
if (parset.get<bool>("writeVTK")) { if (parset.get<bool>("writeVTK")) {
auto const displacement = auto const displacement =
......
...@@ -9,6 +9,10 @@ printDifference = false ...@@ -9,6 +9,10 @@ printDifference = false
writeVTK = true writeVTK = true
[benchmarks.fpi]
enable = false
iterations = 50
[grid] [grid]
refinements = 7 refinements = 7
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment