Skip to content
Snippets Groups Projects
Commit 2944b2f9 authored by podlesny's avatar podlesny
Browse files

added diplacement/state criteria for fpi convergence

parent 6de7438e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/tectonic/spatial-solving/fixedpointiterator.cc
+ 71
− 60
View file @ 2944b2f9
......@@ -5,6 +5,7 @@
#include <dune/common/exceptions.hh>
#include <dune/matrix-vector/axpy.hh>
#include <dune/solvers/norms/sumnorm.hh>
#include <dune/tectonic/utils/reductionfactors.hh>
......@@ -40,6 +41,49 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIte
solverParset_(parset.sub("v.solver")),
errorNorms_(errorNorms) {}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
bool FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::displacementCriterion(const Updaters& updaters, std::vector<Vector>& last_u) const {
bool criterion = true;
std::vector<Vector> u;
updaters.rate_->extractDisplacement(u);
const auto& mats = updaters.rate_->getMatrices();
for (int i=0; i<u.size(); i++) {
const EnergyNorm<Matrix, Vector> ANorm(*mats.elasticity[i]), MNorm(*mats.mass[i]);
const SumNorm<Vector> AMNorm(1.0, ANorm, 1.0, MNorm);
if (u[i].size()==0 || last_u[i].size()==0)
continue;
if (AMNorm.diff(u[i], last_u[i]) >= fixedPointTolerance_) {
criterion = false;
break;
}
}
last_u = u;
return criterion;
}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
bool FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::stateCriterion(const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& newAlpha) const {
bool criterion = true;
for (int i=0; i<alpha.size(); i++) {
if (alpha[i].size()==0 || newAlpha[i].size()==0)
continue;
if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
criterion = false;
break;
}
}
return criterion;
}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
template <class LinearSolver>
FixedPointIterationCounter
......@@ -48,16 +92,14 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
std::vector<Vector>& velocityIterates) {
//std::cout << "FixedPointIterator::run()" << std::endl;
// debugging
/*const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
for (size_t i=0; i<contactCouplings.size(); i++) {
print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
}*/
const auto nBodies = nBodyAssembler_.nGrids();
// debugging
/*const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
for (size_t i=0; i<contactCouplings.size(); i++) {
print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
}*/
FunctionalFactory<std::decay_t<decltype(nBodyAssembler_)>, decltype(globalFriction_), Matrix, Vector> functionalFactory(nBodyAssembler_, globalFriction_, velocityMatrices, velocityRHSs);
functionalFactory.build();
auto functional = functionalFactory.functional();
......@@ -66,91 +108,60 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
size_t fixedPointIteration;
size_t multigridIterations = 0;
std::vector<Vector> last_u;
updaters.rate_->extractOldDisplacement(last_u);
std::vector<ScalarVector> alpha(nBodies);
updaters.state_->extractAlpha(alpha);
Vector totalVelocityIterate, old_v;
nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
old_v = totalVelocityIterate;
std::vector<Vector> v_o;
updaters.rate_->extractOldVelocity(v_o);
Vector total_v_o;
nBodyAssembler_.nodalToTransformed(v_o, total_v_o);
Vector total_v;
nBodyAssembler_.nodalToTransformed(velocityIterates, total_v);
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
//print(alpha, "alpha:");
// contribution from nonlinearity
// update friction
globalFriction_.updateAlpha(alpha);
auto res = solver.solve(solverParset_, totalVelocityIterate);
nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//print(totalVelocityIterate, "totalVelocityIterate:");
//print(velocityIterates, "velocityIterates:");
//DUNE_THROW(Dune::Exception, "Just need to stop here!");
// solve rate problem
auto res = solver.solve(solverParset_, total_v);
multigridIterations += res.iterations;
Vector v_m = old_v;
Vector v_m = total_v_o;
v_m *= 1.0 - lambda_;
Dune::MatrixVector::addProduct(v_m, lambda_, totalVelocityIterate);
Dune::MatrixVector::addProduct(v_m, lambda_, total_v);
// extract relative velocities in mortar basis
std::vector<Vector> v_rel;
split(v_m, v_rel);
//print(v_m, "v_m: ");
//print(v_rel, "v_rel");
//std::cout << "- State problem set" << std::endl;
// solve a state problem
// solve state problem
updaters.state_->solve(v_rel);
//std::cout << "-- Solved" << std::endl;
std::vector<ScalarVector> newAlpha(nBodies);
updaters.state_->extractAlpha(newAlpha);
//print(newAlpha, "new alpha:");
bool breakCriterion = true;
for (int i=0; i<nBodies; i++) {
if (alpha[i].size()==0 || newAlpha[i].size()==0)
continue;
//print(alpha[i], "alpha i:");
//print(newAlpha[i], "new alpha i:");
if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
breakCriterion = false;
//std::cout << "fixedPoint error: " << errorNorms_[i]->diff(alpha[i], newAlpha[i]) << std::endl;
break;
}
}
//std::cout << fixedPointIteration << std::endl;
nBodyAssembler_.postprocess(total_v, velocityIterates);
updaters.rate_->postProcess(velocityIterates);
bool breakCriterion = displacementCriterion(updaters, last_u); //stateCriterion(alpha, newAlpha);
if (lambda_ < 1e-12 or breakCriterion) {
//std::cout << "-FixedPointIteration finished! " << (lambda_ < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
fixedPointIteration++;
break;
}
alpha = newAlpha;
}
//TODO: recently added, might be wrong or superfluous
globalFriction_.updateAlpha(alpha);
//print(alpha, "alpha: ");
//std::cout << "-FixedPointIteration finished with " << fixedPointIteration << " iterations, lambda " << lambda_ << "! " << std::endl;
if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge");
nBodyAssembler_.postprocess(total_v, velocityIterates);
updaters.rate_->postProcess(velocityIterates);
// Cannot use return { fixedPointIteration, multigridIterations };
......
dune/tectonic/spatial-solving/fixedpointiterator.hh
+ 2
− 0
View file @ 2944b2f9
......@@ -45,6 +45,8 @@ class FixedPointIterator {
private:
void split(const Vector& v, std::vector<Vector>& v_rel) const;
bool displacementCriterion(const Updaters& updaters, std::vector<Vector>& last_u) const;
bool stateCriterion(const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& newAlpha) const;
public:
FixedPointIterator(const Dune::ParameterTree& parset,
......
dune/tectonic/time-stepping/rate/rateupdater.cc
+ 6
− 0
View file @ 2944b2f9
......@@ -59,6 +59,12 @@ void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractAccel
acceleration = a;
}
template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
const Matrices<Matrix,2>& RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::getMatrices() const {
return matrices;
}
#include "backward_euler.cc"
#include "newmark.cc"
#include "rateupdater_tmpl.cc"
dune/tectonic/time-stepping/rate/rateupdater.hh
+ 2
− 0
View file @ 2944b2f9
......@@ -36,6 +36,8 @@ class RateUpdater {
void extractOldVelocity(std::vector<Vector>& velocity) const;
void extractAcceleration(std::vector<Vector>& acceleration) const;
const Matrices<Matrix,2>& getMatrices() const;
std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> virtual clone() const = 0;
protected:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment