From 2614079b403991d86bbcd15bd4db9be985abf9f9 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Thu, 10 Jul 2014 23:15:06 +0200
Subject: [PATCH] [Kill] Remove the sliding block problem
---
src/Makefile.am | 6 +-
src/sliding-block-data/boundaryconditions.py | 16 -
src/sliding-block-data/mybody.hh | 54 --
src/sliding-block-data/mygeometry.hh | 27 -
.../myglobalfrictiondata.hh | 45 --
src/sliding-block-data/mygrid.hh | 36 --
src/sliding-block-data/parset.cfg | 69 ---
src/sliding-block.cc | 463 ------------------
8 files changed, 1 insertion(+), 715 deletions(-)
delete mode 100644 src/sliding-block-data/boundaryconditions.py
delete mode 100644 src/sliding-block-data/mybody.hh
delete mode 100644 src/sliding-block-data/mygeometry.hh
delete mode 100644 src/sliding-block-data/myglobalfrictiondata.hh
delete mode 100644 src/sliding-block-data/mygrid.hh
delete mode 100644 src/sliding-block-data/parset.cfg
delete mode 100644 src/sliding-block.cc
diff --git a/src/Makefile.am b/src/Makefile.am
index 6864003b..a939e404 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = sand-wedge sliding-block-2D
+bin_PROGRAMS = sand-wedge
common_sources = \
assemblers.cc \
@@ -12,10 +12,6 @@ sand_wedge_SOURCES = $(common_sources) sand-wedge.cc
sand_wedge_CPPFLAGS = \
$(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
-Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DDIM=2
-sliding_block_2D_SOURCES = $(common_sources) sliding-block.cc
-sliding_block_2D_CPPFLAGS = \
- $(AM_CPPFLAGS) \
- -Ddatadir=\"$(abs_srcdir)/sliding-block-data/\" -DDIM=2
# Some are for clang, others are for gcc
AM_CXXFLAGS = \
diff --git a/src/sliding-block-data/boundaryconditions.py b/src/sliding-block-data/boundaryconditions.py
deleted file mode 100644
index 04a1799a..00000000
--- a/src/sliding-block-data/boundaryconditions.py
+++ /dev/null
@@ -1,16 +0,0 @@
-class neumannCondition:
- def __call__(self, relativeTime):
- return 0
-
-class velocityDirichletCondition:
- def __call__(self, relativeTime):
- if relativeTime <= 0.25:
- factor = 4*relativeTime;
- else:
- factor = 1
- return 2e-4 * factor
-
-Functions = {
- 'neumannCondition' : neumannCondition(),
- 'velocityDirichletCondition' : velocityDirichletCondition()
-}
diff --git a/src/sliding-block-data/mybody.hh b/src/sliding-block-data/mybody.hh
deleted file mode 100644
index b058898b..00000000
--- a/src/sliding-block-data/mybody.hh
+++ /dev/null
@@ -1,54 +0,0 @@
-#ifndef SRC_SLIDING_BLOCK_DATA_MYBODY_HH
-#define SRC_SLIDING_BLOCK_DATA_MYBODY_HH
-
-#include <dune/common/fvector.hh>
-
-#include <dune/fufem/functions/constantfunction.hh>
-
-#include <dune/tectonic/body.hh>
-#include <dune/tectonic/gravity.hh>
-
-#include "mygeometry.hh"
-
-template <int dimension> class MyBody : public Body<dimension> {
- using typename Body<dimension>::ScalarFunction;
- using typename Body<dimension>::VectorField;
- using MyConstantFunction = ConstantFunction<
- Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
-
-public:
- MyBody(Dune::ParameterTree const &parset)
- : poissonRatio_(parset.get<double>("body.poissonRatio")),
- youngModulus_(parset.get<double>("body.youngModulus")),
- shearViscosityField_(parset.get<double>("body.shearViscosity")),
- bulkViscosityField_(parset.get<double>("body.bulkViscosity")),
- densityField_(parset.get<double>("body.density")),
- gravityField_(densityField_, MyGeometry::zenith,
- parset.get<double>("gravity")) {}
-
- double getPoissonRatio() const override { return poissonRatio_; }
-
- double getYoungModulus() const override { return youngModulus_; }
-
- ScalarFunction const &getShearViscosityField() const override {
- return shearViscosityField_;
- }
- ScalarFunction const &getBulkViscosityField() const override {
- return bulkViscosityField_;
- }
-
- ScalarFunction const &getDensityField() const override {
- return densityField_;
- }
- VectorField const &getGravityField() const override { return gravityField_; }
-
-private:
- double const poissonRatio_;
- double const youngModulus_;
- MyConstantFunction const shearViscosityField_;
- MyConstantFunction const bulkViscosityField_;
- MyConstantFunction const densityField_;
- Gravity<dimension> const gravityField_;
-};
-
-#endif
diff --git a/src/sliding-block-data/mygeometry.hh b/src/sliding-block-data/mygeometry.hh
deleted file mode 100644
index 96cd6534..00000000
--- a/src/sliding-block-data/mygeometry.hh
+++ /dev/null
@@ -1,27 +0,0 @@
-#ifndef SRC_SLIDING_BLOCK_DATA_MYGEOMETRY_HH
-#define SRC_SLIDING_BLOCK_DATA_MYGEOMETRY_HH
-
-#include <dune/common/fassign.hh>
-#include <dune/common/fvector.hh>
-
-namespace MyGeometry {
-namespace {
- using LocalVector = Dune::FieldVector<double, 2>;
-
- // kludge because fieldvectors have no initialiser_list constructor, see
- // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
- LocalVector generateVector(double x, double y) {
- LocalVector tmp;
- tmp <<= x, y;
- return tmp;
- }
-}
-
-LocalVector const A = generateVector(0, 0);
-LocalVector const B = generateVector(5, 0);
-LocalVector const C = generateVector(5, 1);
-LocalVector const D = generateVector(0, 1);
-
-LocalVector const zenith = generateVector(0, 1);
-}
-#endif
diff --git a/src/sliding-block-data/myglobalfrictiondata.hh b/src/sliding-block-data/myglobalfrictiondata.hh
deleted file mode 100644
index 8487a125..00000000
--- a/src/sliding-block-data/myglobalfrictiondata.hh
+++ /dev/null
@@ -1,45 +0,0 @@
-#ifndef SRC_SLIDING_BLOCK_DATA_MYGLOBALFRICTIONDATA_HH
-#define SRC_SLIDING_BLOCK_DATA_MYGLOBALFRICTIONDATA_HH
-
-#include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/fufem/functions/constantfunction.hh>
-
-#include <dune/tectonic/globalfrictiondata.hh>
-
-#include "mygeometry.hh"
-
-template <int dimension>
-class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
-private:
- using typename GlobalFrictionData<dimension>::VirtualFunction;
-
-public:
- MyGlobalFrictionData(Dune::ParameterTree const &parset)
- : C_(parset.get<double>("C")),
- L_(parset.get<double>("L")),
- V0_(parset.get<double>("V0")),
- a_(parset.get<double>("a")),
- b_(parset.get<double>("b")),
- mu0_(parset.get<double>("mu0")) {}
-
- double virtual const &C() const override { return C_; }
- double virtual const &L() const override { return L_; }
- double virtual const &V0() const override { return V0_; }
- VirtualFunction virtual const &a() const override { return a_; }
- VirtualFunction virtual const &b() const override { return b_; }
- double virtual const &mu0() const override { return mu0_; }
-
-private:
- using MyConstantFunction = ConstantFunction<
- Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
-
- double const C_;
- double const L_;
- double const V0_;
- MyConstantFunction const a_;
- MyConstantFunction const b_;
- double const mu0_;
-};
-#endif
diff --git a/src/sliding-block-data/mygrid.hh b/src/sliding-block-data/mygrid.hh
deleted file mode 100644
index d38616c1..00000000
--- a/src/sliding-block-data/mygrid.hh
+++ /dev/null
@@ -1,36 +0,0 @@
-#ifndef SRC_SLIDING_BLOCK_DATA_MYGRID_HH
-#define SRC_SLIDING_BLOCK_DATA_MYGRID_HH
-
-#include <dune/grid/common/gridfactory.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include "mygeometry.hh"
-
-template <class Grid> std::shared_ptr<Grid> constructGrid() {
- std::array<unsigned int, 2> elements = { { 5, 1 } };
- return Dune::StructuredGridFactory<Grid>::createSimplexGrid(
- MyGeometry::A, MyGeometry::C, elements);
-}
-
-template <class GridView> class MyFaces {
-private:
- bool isClose(double a, double b) {
- return std::abs(a - b) < 1e-14;
- };
-
-public:
- MyFaces(GridView const &gridView) : lower(gridView), upper(gridView) {
- lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
- return isClose(MyGeometry::A[1], in.geometry().center()[1]);
- });
-
- upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
- return isClose(MyGeometry::C[1], in.geometry().center()[1]);
- });
- }
-
- BoundaryPatch<GridView> lower;
- BoundaryPatch<GridView> upper;
-};
-#endif
diff --git a/src/sliding-block-data/parset.cfg b/src/sliding-block-data/parset.cfg
deleted file mode 100644
index 7f477c49..00000000
--- a/src/sliding-block-data/parset.cfg
+++ /dev/null
@@ -1,69 +0,0 @@
-# -*- mode:conf -*-
-gravity = 9.81
-
-[io]
-printProgress = false
-writeVTK = false
-
-[problem]
-finalTime = 15
-
-[body]
-youngModulus = 5e7
-poissonRatio = 0.3
-density = 5000
-shearViscosity = 0
-bulkViscosity = 0
-
-[boundary.friction]
-C = 0.0
-mu0 = 0.6
-a = 0.010
-b = 0.015
-V0 = 1e-6
-L = 1e-5
-initialState = 4.54e-05 # = exp(-10) = theta, so that alpha = -10
-stateModel = Dieterich
-
-[timeSteps]
-number = 10000
-scheme = newmark
-
-[grid]
-refinements = 4
-
-[u0.solver]
-tolerance = 1e-10
-maximumIterations = 100000
-verbosity = quiet
-
-[a0.solver]
-tolerance = 1e-10
-maximumIterations = 100000
-verbosity = quiet
-
-[v.solver]
-tolerance = 1e-10
-maximumIterations = 100000
-verbosity = quiet
-
-[v.fpi]
-tolerance = 1e-10
-maximumIterations = 10000
-relaxation = 0.5
-requiredReduction = 0.5
-
-[solver.tnnmg.linear]
-maxiumumIterations = 100000
-tolerance = 1e-10
-pre = 3
-cycle = 1 # 1 = V, 2 = W, etc.
-post = 3
-
-[solver.tnnmg.main]
-pre = 1
-multi = 5 # number of multigrid steps
-post = 0
-
-[localsolver]
-steps = 1
diff --git a/src/sliding-block.cc b/src/sliding-block.cc
deleted file mode 100644
index ef52fbe3..00000000
--- a/src/sliding-block.cc
+++ /dev/null
@@ -1,463 +0,0 @@
-#include <Python.h>
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifndef HAVE_PYTHON
-#error Python is required
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#ifndef datadir
-#error datadir unset
-#endif
-
-#ifndef DIM
-#error DIM unset
-#endif
-
-#if !HAVE_ALUGRID
-#error ALUGRID is required
-#endif
-
-#include <cmath>
-#include <exception>
-#include <fstream>
-#include <iostream>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/exceptions.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wignored-qualifiers"
-#include <dune/grid/alugrid.hh>
-#pragma clang diagnostic pop
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wdeprecated-declarations"
-#include <dune/fufem/boundarypatch.hh>
-#pragma clang diagnostic pop
-
-#include <dune/fufem/dunepython.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/sharedpointermap.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/norms/sumnorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/solvers/solver.hh>
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/myblockproblem.hh>
-#include <dune/tectonic/globalnonlinearity.hh>
-
-#include "assemblers.hh"
-#include "enum_parser.cc"
-#include "enum_scheme.cc"
-#include "enum_state_model.cc"
-#include "enum_verbosity.cc"
-#include "enums.hh"
-#include "friction_writer.hh"
-#include "sliding-block-data/mybody.hh"
-#include "sliding-block-data/mygeometry.hh"
-#include "sliding-block-data/myglobalfrictiondata.hh"
-#include "sliding-block-data/mygrid.hh"
-#include "solverfactory.hh"
-#include "state.hh"
-#include "timestepping.hh"
-#include "vtk.hh"
-size_t const dims = DIM;
-
-void initPython() {
- Python::start();
-
- Python::run("import sys");
- Python::run("sys.path.append('" datadir "')");
-}
-
-int main(int argc, char *argv[]) {
- try {
- Dune::ParameterTree parset;
- Dune::ParameterTreeParser::readINITree(datadir "/parset.cfg", parset);
- Dune::ParameterTreeParser::readOptions(argc, argv, parset);
-
- // {{{ Set up grid
- using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
- auto grid = constructGrid<Grid>();
-
- auto const refinements = parset.get<size_t>("grid.refinements");
- grid->globalRefine(refinements);
- size_t const fineVertexCount = grid->size(grid->maxLevel(), dims);
-
- using GridView = Grid::LeafGridView;
- GridView const leafView = grid->leafView();
- // }}}
-
- // Set up myFaces
- MyFaces<GridView> myFaces(leafView);
-
- // Neumann boundary
- BoundaryPatch<GridView> const neumannBoundary(leafView);
-
- // Frictional Boundary
- BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
- Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
- frictionalBoundary.getVertices(frictionalNodes);
-
- // Dirichlet Boundary
- Dune::BitSetVector<dims> noNodes(fineVertexCount);
- Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
- for (size_t i = 0; i < fineVertexCount; ++i) {
- if (myFaces.lower.containsVertex(i))
- dirichletNodes[i][1] = true;
-
- if (myFaces.upper.containsVertex(i))
- dirichletNodes[i] = true;
- }
-
- // Set up functions for time-dependent boundary conditions
- using Function = Dune::VirtualFunction<double, double>;
- using FunctionMap = SharedPointerMap<std::string, Function>;
- FunctionMap functions;
- {
- initPython();
- Python::import("boundaryconditions")
- .get("Functions")
- .toC<typename FunctionMap::Base>(functions);
- }
- auto const &velocityDirichletFunction =
- functions.get("velocityDirichletCondition"),
- &neumannFunction = functions.get("neumannCondition");
-
- using MyAssembler = MyAssembler<GridView, dims>;
- using Matrix = MyAssembler::Matrix;
- using LocalMatrix = Matrix::block_type;
- using Vector = MyAssembler::Vector;
- using LocalVector = Vector::block_type;
- using ScalarMatrix = MyAssembler::ScalarMatrix;
- using ScalarVector = MyAssembler::ScalarVector;
-
- MyAssembler myAssembler(leafView);
-
- MyBody<dims> const body(parset);
-
- Matrix A, C, M;
- myAssembler.assembleElasticity(body.getYoungModulus(),
- body.getPoissonRatio(), A);
- myAssembler.assembleViscosity(body.getShearViscosityField(),
- body.getBulkViscosityField(), C);
- myAssembler.assembleMass(body.getDensityField(), M);
- EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
- // Q: Does it make sense to weigh them in this manner?
- SumNorm<Vector> const AMNorm(1.0, ANorm, 1.0, MNorm);
-
- ScalarMatrix frictionalBoundaryMass;
- myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
- frictionalBoundaryMass);
- EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
- frictionalBoundaryMass);
-
- // Assemble forces
- Vector gravityFunctional;
- myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
-
- // Problem formulation: right-hand side
- auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
- myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
- _relativeTime);
- _ell += gravityFunctional;
- };
- Vector ell(fineVertexCount);
- computeExternalForces(0.0, ell);
-
- // {{{ Initial conditions
- using LinearFactory = SolverFactory<
- dims, BlockNonlinearTNNMGProblem<ConvexProblem<
- ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
- Grid>;
- ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
-
- // TODO: clean up once generic lambdas arrive
- auto const solveLinearProblem = [&](
- Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
- Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> const &_norm,
- Dune::ParameterTree const &_localParset) {
- LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
- refinements, *grid, _dirichletNodes);
-
- typename LinearFactory::ConvexProblem convexProblem(
- 1.0, _matrix, zeroNonlinearity, _rhs, _x);
- typename LinearFactory::BlockProblem problem(parset, convexProblem);
-
- auto multigridStep = factory.getSolver();
- multigridStep->setProblem(_x, problem);
- LoopSolver<Vector> solver(
- multigridStep, _localParset.get<size_t>("maximumIterations"),
- _localParset.get<double>("tolerance"), &_norm,
- _localParset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
-
- solver.preprocess();
- solver.solve();
- };
-
- // Solve the stationary problem
- Vector u_initial(fineVertexCount);
- u_initial = 0.0;
- solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
- parset.sub("u0.solver"));
- Vector ur_initial(fineVertexCount);
- ur_initial = 0.0;
-
- ScalarVector alpha_initial(fineVertexCount);
- alpha_initial =
- std::log(parset.get<double>("boundary.friction.initialState"));
- ScalarVector normalStress(fineVertexCount);
- myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
- body.getYoungModulus(),
- body.getPoissonRatio(), u_initial);
- MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
- auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
- frictionalBoundary, frictionInfo, normalStress);
- myGlobalNonlinearity->updateLogState(alpha_initial);
-
- Vector v_initial(fineVertexCount);
- v_initial = 0.0;
- {
- // Prescribe a homogeneous velocity field in the x-direction
- // This way, the initial condition for v and the Dirichlet
- // condition match up at t=0
- double v_initial_const;
- velocityDirichletFunction.evaluate(0.0, v_initial_const);
- for (auto &x : v_initial)
- x[0] = v_initial_const;
- }
- Vector vr_initial(fineVertexCount);
- vr_initial = 0.0;
-
- Vector a_initial(fineVertexCount);
- a_initial = 0.0;
- {
- // We solve Ma = ell - [Au + Cv + Psi(v)]
- Vector accelerationRHS(fineVertexCount);
- {
- accelerationRHS = 0.0;
- Arithmetic::addProduct(accelerationRHS, A, u_initial);
- Arithmetic::addProduct(accelerationRHS, C, v_initial);
- // NOTE: We assume differentiability of Psi at v0 here!
- myGlobalNonlinearity->addGradient(v_initial, accelerationRHS);
- accelerationRHS *= -1.0;
- accelerationRHS += ell;
- }
- solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
- parset.sub("a0.solver"));
- }
- // }}}
-
- Vector vertexCoordinates(fineVertexCount);
- {
- Dune::MultipleCodimMultipleGeomTypeMapper<
- GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
- for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) {
- auto const geometry = it->geometry();
- assert(geometry.corners() == 1);
- vertexCoordinates[vertexMapper.map(*it)] = geometry.corner(0);
- }
- }
- FrictionWriter<ScalarVector, Vector> frictionWriter(
- vertexCoordinates, frictionalNodes, "friction",
- [](LocalVector const &x) { return x[0]; });
- auto const reportFriction = [&](Vector const &_u, Vector const &_v,
- ScalarVector const &_alpha) {
- ScalarVector c;
- myGlobalNonlinearity->coefficientOfFriction(_v, c);
- frictionWriter.writeKinetics(_u, _v);
- frictionWriter.writeOther(c, _alpha);
- };
- reportFriction(ur_initial, vr_initial, alpha_initial);
-
- MyVTKWriter<typename MyAssembler::VertexBasis,
- typename MyAssembler::CellBasis> const
- vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
-
- if (parset.get<bool>("io.writeVTK")) {
- ScalarVector stress;
- myAssembler.assembleVonMisesStress(
- body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
- vtkWriter.write(0, ur_initial, vr_initial, alpha_initial, stress);
- }
-
- // Set up TNNMG solver
- using NonlinearFactory =
- SolverFactory<dims, MyBlockProblem<ConvexProblem<
- GlobalNonlinearity<Matrix, Vector>, Matrix>>,
- Grid>;
- NonlinearFactory factory(parset.sub("solver.tnnmg"), refinements, *grid,
- dirichletNodes);
- auto multigridStep = factory.getSolver();
-
- std::fstream iterationWriter("iterations", std::fstream::out),
- relaxationWriter("relaxation", std::fstream::out);
-
- auto timeSteppingScheme = initTimeStepper(
- parset.get<Config::scheme>("timeSteps.scheme"),
- velocityDirichletFunction, dirichletNodes, M, A, C, u_initial,
- ur_initial, v_initial, vr_initial, a_initial);
- auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
- parset.get<Config::stateModel>("boundary.friction.stateModel"),
- alpha_initial, frictionalNodes,
- parset.get<double>("boundary.friction.L"));
-
- Vector v = v_initial;
- Vector v_m(fineVertexCount);
- ScalarVector alpha(fineVertexCount);
-
- auto const timeSteps = parset.get<size_t>("timeSteps.number"),
- maximumStateFPI = parset.get<size_t>("v.fpi.maximumIterations"),
- maximumIterations =
- parset.get<size_t>("v.solver.maximumIterations");
- auto const tau = parset.get<double>("problem.finalTime") / timeSteps,
- tolerance = parset.get<double>("v.solver.tolerance"),
- fixedPointTolerance = parset.get<double>("v.fpi.tolerance"),
- relaxation = parset.get<double>("v.fpi.relaxation"),
- requiredReduction =
- parset.get<double>("v.fpi.requiredReduction");
- auto const printProgress = parset.get<bool>("io.printProgress");
- auto const verbosity =
- parset.get<Solver::VerbosityMode>("v.solver.verbosity");
- for (size_t timeStep = 1; timeStep <= timeSteps; ++timeStep) {
- if (printProgress)
- std::cout << std::setw(7) << timeStep << " " << std::flush;
-
- stateUpdater->nextTimeStep();
- timeSteppingScheme->nextTimeStep();
-
- auto const relativeTime = double(timeStep) / double(timeSteps);
- computeExternalForces(relativeTime, ell);
-
- Matrix velocityMatrix;
- Vector velocityRHS(fineVertexCount);
- Vector velocityIterate(fineVertexCount);
-
- stateUpdater->setup(tau);
- timeSteppingScheme->setup(ell, tau, relativeTime, velocityRHS,
- velocityIterate, velocityMatrix);
-
- LoopSolver<Vector> velocityProblemSolver(multigridStep, maximumIterations,
- tolerance, &AMNorm, verbosity,
- false); // absolute error
-
- size_t iterationCounter;
- auto solveVelocityProblem = [&](Vector &_velocityIterate,
- ScalarVector const &_alpha) {
- myGlobalNonlinearity->updateLogState(_alpha);
-
- // NIT: Do we really need to pass u here?
- typename NonlinearFactory::ConvexProblem convexProblem(
- 1.0, velocityMatrix, *myGlobalNonlinearity, velocityRHS,
- _velocityIterate);
- typename NonlinearFactory::BlockProblem velocityProblem(parset,
- convexProblem);
- multigridStep->setProblem(_velocityIterate, velocityProblem);
-
- velocityProblemSolver.preprocess();
- velocityProblemSolver.solve();
- iterationCounter = velocityProblemSolver.getResult().iterations;
- };
-
- Vector u;
- Vector v_saved;
- ScalarVector alpha_saved;
- double lastStateCorrection;
- for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
- timeSteppingScheme->extractOldVelocity(v_m);
- v_m *= 0.5;
- Arithmetic::addProduct(v_m, 0.5, v);
-
- stateUpdater->solve(v_m);
- stateUpdater->extractLogState(alpha);
-
- if (stateFPI == 1)
- relaxationWriter << "N ";
- else {
- double const stateCorrection =
- stateEnergyNorm.diff(alpha, alpha_saved);
- if (stateFPI <= 2 // lastStateCorrection is only set for stateFPI > 2
- or stateCorrection < requiredReduction * lastStateCorrection)
- relaxationWriter << "N ";
- else {
- alpha *= (1.0 - relaxation);
- Arithmetic::addProduct(alpha, relaxation, alpha_saved);
- relaxationWriter << "Y ";
- }
- lastStateCorrection = stateCorrection;
- }
-
- 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);
- if (velocityCorrection < fixedPointTolerance)
- break;
- }
- if (stateFPI == maximumStateFPI)
- DUNE_THROW(Dune::Exception, "FPI failed to converge");
-
- alpha_saved = alpha;
- v_saved = v;
- }
- if (printProgress)
- std::cout << std::endl;
-
- Vector ur;
- Vector vr;
- timeSteppingScheme->postProcessRelativeQuantities();
- timeSteppingScheme->extractRelativeDisplacement(ur);
- timeSteppingScheme->extractRelativeVelocity(vr);
-
- reportFriction(ur, vr, alpha);
- iterationWriter << std::endl;
- relaxationWriter << std::endl;
-
- if (parset.get<bool>("io.writeVTK")) {
- ScalarVector stress;
- myAssembler.assembleVonMisesStress(body.getYoungModulus(),
- body.getPoissonRatio(), u, stress);
- vtkWriter.write(timeStep, ur, vr, alpha, stress);
- }
- }
- iterationWriter.close();
- relaxationWriter.close();
-
- Python::stop();
- }
- catch (Dune::Exception &e) {
- Dune::derr << "Dune reported error: " << e << std::endl;
- }
- catch (std::exception &e) {
- std::cerr << "Standard exception: " << e.what() << std::endl;
- }
-}
--
GitLab