From 7db7854f85ee95a7f1406aa14cf8cf21421c1dc5 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Tue, 9 Feb 2021 00:00:42 +0100
Subject: [PATCH] work in progress: not working remainder
---
dune/tectonic/data-structures/body/body.cc | 8 +
.../tectonic/data-structures/program_state.hh | 182 ++------------
.../spatial-solving/fixedpointiterator.cc | 162 ++----------
.../spatial-solving/fixedpointiterator.hh | 8 +-
src/foam/foam-2D.cfg | 8 +-
src/foam/foam.cc | 230 ++++++++++++++----
src/foam/foam.cfg | 62 ++---
src/strikeslip/strikeslip-2D.cfg | 4 +-
src/strikeslip/strikeslip.cfg | 50 ++--
9 files changed, 290 insertions(+), 424 deletions(-)
diff --git a/dune/tectonic/data-structures/body/body.cc b/dune/tectonic/data-structures/body/body.cc
index cd63c1d5..83fb4fa8 100644
--- a/dune/tectonic/data-structures/body/body.cc
+++ b/dune/tectonic/data-structures/body/body.cc
@@ -5,6 +5,8 @@
#include <dune/common/hybridutilities.hh>
#include "body.hh"
+#include <dune/tectonic/utils/debugutils.hh>
+
// -------------------------------
// --- LeafBody Implementation ---
// -------------------------------
@@ -50,11 +52,17 @@ void LeafBody<GridTEMPLATE, VectorType>::assemble() {
// assemble matrices
assembler_->assembleElasticity(bodyData_->getYoungModulus(), bodyData_->getPoissonRatio(), *matrices_.elasticity);
+ //print(*matrices_.elasticity, "elasticity");
+
assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
+ //print(*matrices_.damping, "viscosity");
+
assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
+ //print(*matrices_.mass, "mass");
// assemble forces
assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
+ //print(gravityFunctional_, "gravity");
}
// setter and getter
diff --git a/dune/tectonic/data-structures/program_state.hh b/dune/tectonic/data-structures/program_state.hh
index 9e1f63f5..d6a921c2 100644
--- a/dune/tectonic/data-structures/program_state.hh
+++ b/dune/tectonic/data-structures/program_state.hh
@@ -19,8 +19,10 @@
#include "network/contactnetwork.hh"
#include "matrices.hh"
#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+#include "../spatial-solving/makelinearsolver.hh"
+#include "../spatial-solving/nonlinearsolver.hh"
#include "../spatial-solving/solverfactory.hh"
-#include "../spatial-solving/solverfactory.cc"
+#include "../spatial-solving/functionalfactory.hh"
#include "../spatial-solving/tnnmg/functional.hh"
#include "../spatial-solving/tnnmg/zerononlinearity.hh"
#include "../utils/debugutils.hh"
@@ -99,60 +101,15 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
// Set up initial conditions
template <class ContactNetwork>
void setupInitialConditions(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
- using Matrix = typename ContactNetwork::Matrix;
- const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
std::cout << "-- setupInitialConditions --" << std::endl;
std::cout << "----------------------------" << std::endl;
- // make linear solver for linear correction in TNNMGStep
- using Norm = EnergyNorm<Matrix, Vector>;
- using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
- using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
-
- /* const auto& preconditionerParset = parset.sub("solver.tnnmg.preconditioner");
-
- Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
- Preconditioner preconditioner(preconditionerParset, contactNetwork, activeLevels);
- preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
-
- std::cout << "Building preconditioner..." << std::endl;
- preconditioner.build();
-
- auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
- cgStep->setPreconditioner(preconditioner);
-
- Norm norm(*cgStep);
-
- auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset.get<int>("solver.tnnmg.main.multi"),
- parset.get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"),
- norm, Solver::QUIET);
-
-*/
- // set multigrid solver
- auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
-
- using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
- using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
-
- TransferOperators transfer(contactNetwork.nLevels()-1);
- for (size_t i=0; i<transfer.size(); i++) {
- transfer[i] = std::make_shared<TransferOperator>();
- transfer[i]->setup(contactNetwork, i, i+1);
- }
-
- // Remove any recompute filed so that initially the full transferoperator is assembled
- for (size_t i=0; i<transfer.size(); i++)
- std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
-
- auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
- linearMultigridStep->setMGType(1, 3, 3);
- linearMultigridStep->setSmoother(smoother);
- linearMultigridStep->setTransferOperators(transfer);
-
- Norm norm(*linearMultigridStep);
+ using Matrix = typename ContactNetwork::Matrix;
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
- auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, parset.get<int>("solver.tnnmg.main.multi"), parset.get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+ auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
+ auto nonlinearity = ZeroNonlinearity();
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
@@ -160,126 +117,17 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
Dune::ParameterTree const &_localParset) {
- std::vector<const Matrix*> matrices_ptr(_matrices.size());
- for (size_t i=0; i<matrices_ptr.size(); i++) {
- matrices_ptr[i] = _matrices[i].get();
- }
-
- /*std::vector<Matrix> matrices(velocityMatrices.size());
- std::vector<Vector> rhs(velocityRHSs.size());
- for (size_t i=0; i<globalFriction_.size(); i++) {
- matrices[i] = velocityMatrices[i];
- rhs[i] = velocityRHSs[i];
-
- globalFriction_[i]->addHessian(v_rel[i], matrices[i]);
- globalFriction_[i]->addGradient(v_rel[i], rhs[i]);
-
- matrices_ptr[i] = &matrices[i];
- }*/
-
- // assemble full global contact problem
- Matrix bilinearForm;
-
- nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
-
- Vector totalRhs, oldTotalRhs;
- nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
- oldTotalRhs = totalRhs;
-
- Vector totalX, oldTotalX;
- nBodyAssembler.nodalToTransformed(_x, totalX);
- oldTotalX = totalX;
-
- // get lower and upper obstacles
- const auto& totalObstacles = nBodyAssembler.totalObstacles_;
- Vector lower(totalObstacles.size());
- Vector upper(totalObstacles.size());
-
- for (size_t j=0; j<totalObstacles.size(); ++j) {
- const auto& totalObstaclesj = totalObstacles[j];
- auto& lowerj = lower[j];
- auto& upperj = upper[j];
- for (size_t d=0; d<dims; ++d) {
- lowerj[d] = totalObstaclesj[d][0];
- upperj[d] = totalObstaclesj[d][1];
- }
- }
-
- // print problem
- /* print(bilinearForm, "bilinearForm");
- print(totalRhs, "totalRhs");
- print(_dirichletNodes, "ignore");
- print(totalObstacles, "totalObstacles");
- print(lower, "lower");
- print(upper, "upper");*/
-
- // set up functional
- using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
- Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper); //TODO
-
- // set up TNMMG solver
- using Factory = SolverFactory<Functional, BitVector>;
- //Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
- Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
- factory.build(linearSolver);
-
- /* std::vector<BitVector> bodyDirichletNodes;
- nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
- for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
- print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
- }*/
-
- /* print(bilinearForm, "matrix: ");
- print(totalX, "totalX: ");
- print(totalRhs, "totalRhs: ");*/
-
- auto tnnmgStep = factory.step();
- factory.setProblem(totalX);
-
- const EnergyNorm<Matrix, Vector> norm(bilinearForm);
-
- LoopSolver<Vector> solver(
- *tnnmgStep.get(), _localParset.get<size_t>("maximumIterations"),
- _localParset.get<double>("tolerance"), norm,
- _localParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
-
- solver.preprocess();
- solver.solve();
-
- std::cout << "ProgramState: Energy of TNNMG solution: " << J(tnnmgStep->getSol()) << std::endl;
-
- nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
-
- Vector res = totalRhs;
- bilinearForm.mmv(tnnmgStep->getSol(), res);
- std::cout << "TNNMG Res - energy norm: " << norm.operator ()(res) << std::endl;
-
- // TODO: remove after debugging
- /* using DeformedGridType = typename LevelContactNetwork<GridType, dims>::DeformedGridType;
- using OldLinearFactory = SolverFactoryOld<DeformedGridType, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
- OldLinearFactory oldFactory(parset.sub("solver.tnnmg"), nBodyAssembler, _dirichletNodes);
-
- auto oldStep = oldFactory.getStep();
- oldStep->setProblem(bilinearForm, oldTotalX, oldTotalRhs);
-
- LoopSolver<Vector> oldSolver(
- oldStep.get(), _localParset.get<size_t>("maximumIterations"),
- _localParset.get<double>("tolerance"), &norm,
- _localParset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
+ Vector totalX;
+ nBodyAssembler.nodalToTransformed(_x, totalX);
- oldSolver.preprocess();
- oldSolver.solve();
+ FunctionalFactory<std::decay_t<decltype(nBodyAssembler)>, decltype(nonlinearity), Matrix, Vector> functionalFactory(nBodyAssembler, nonlinearity, _matrices, _rhs);
+ functionalFactory.build();
+ auto functional = functionalFactory.functional();
- Vector oldRes = totalRhs;
- bilinearForm.mmv(oldStep->getSol(), oldRes);
- std::cout << "Old Res - energy norm: " << norm.operator ()(oldRes) << std::endl;*/
+ NonlinearSolver<std::remove_reference_t<decltype(*functional)>, BitVector> solver(parset.sub("solver.tnnmg"), linearSolver, functional, _dirichletNodes);
+ solver.solve(_localParset, totalX);
- // print(tnnmgStep->getSol(), "TNNMG Solution: ");
- /* print(oldStep->getSol(), "Old Solution: ");
- auto diff = tnnmgStep->getSol();
- diff -= oldStep->getSol();
- std::cout << "Energy norm: " << norm.operator ()(diff) << std::endl;*/
+ nBodyAssembler.postprocess(totalX, _x);
};
timeStep = parset.get<size_t>("initialTime.timeStep");
diff --git a/dune/tectonic/spatial-solving/fixedpointiterator.cc b/dune/tectonic/spatial-solving/fixedpointiterator.cc
index a36d76f2..bef1fea2 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator.cc
+++ b/dune/tectonic/spatial-solving/fixedpointiterator.cc
@@ -6,41 +6,13 @@
#include <dune/matrix-vector/axpy.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-
-#include <dune/contact/assemblers/nbodyassembler.hh>
-#include <dune/contact/common/dualbasisadapter.hh>
-
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
-
-#include <dune/functions/gridfunctions/gridfunction.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-#include <dune/geometry/type.hh>
-#include <dune/geometry/referenceelements.hh>
-
-#include <dune/fufem/functions/basisgridfunction.hh>
-
-#include "../data-structures/enums.hh"
-#include "../data-structures/enumparser.hh"
-
+#include <dune/tectonic/utils/reductionfactors.hh>
+#include "functionalfactory.hh"
+#include "nonlinearsolver.hh"
-#include "../utils/tobool.hh"
#include "../utils/debugutils.hh"
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-
-#include "../spatial-solving/contact/nbodycontacttransfer.hh"
-
-#include "tnnmg/functional.hh"
-#include "tnnmg/zerononlinearity.hh"
-
-#include <dune/tectonic/utils/reductionfactors.hh>
-
#include "fixedpointiterator.hh"
void FixedPointIterationCounter::operator+=(
@@ -65,9 +37,7 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIte
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
lambda_(parset.get<double>("v.fpi.lambda")),
- velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
- velocityTolerance_(parset.get<double>("v.solver.tolerance")),
- verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
+ solverParset_(parset.sub("v.solver")),
errorNorms_(errorNorms) {}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
@@ -80,110 +50,28 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
//std::cout << "FixedPointIterator::run()" << std::endl;
+ 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:");
}*/
- const auto nBodies = nBodyAssembler_.nGrids();
-
- std::vector<const Matrix*> matrices_ptr(nBodies);
- for (int i=0; i<nBodies; i++) {
- matrices_ptr[i] = &velocityMatrices[i];
- }
-
- // assemble full global contact problem
- Matrix bilinearForm;
- nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
+ FunctionalFactory<std::decay_t<decltype(nBodyAssembler_)>, decltype(globalFriction_), Matrix, Vector> functionalFactory(nBodyAssembler_, globalFriction_, velocityMatrices, velocityRHSs);
+ functionalFactory.build();
+ auto functional = functionalFactory.functional();
- //print(bilinearForm, "bilinearForm:");
-
- Vector totalRhs;
- nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
-
- //print(totalRhs, "totalRhs:");
-
- // get lower and upper obstacles
- const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
- Vector lower(totalObstacles.size());
- Vector upper(totalObstacles.size());
-
- for (size_t j=0; j<totalObstacles.size(); ++j) {
- const auto& totalObstaclesj = totalObstacles[j];
- auto& lowerj = lower[j];
- auto& upperj = upper[j];
- for (size_t d=0; d<dims; ++d) {
- lowerj[d] = totalObstaclesj[d][0];
- upperj[d] = totalObstaclesj[d][1];
- }
- }
-
- //print(totalObstacles, "totalObstacles:");
-
- //print(lower, "lower obstacles:");
- //print(upper, "upper obstacles:");
-
- // compute velocity obstacles
- /*Vector vLower, vUpper;
- std::vector<Vector> u0, v0;
- updaters.rate_->extractOldVelocity(v0);
- updaters.rate_->extractOldDisplacement(u0);
-
- Vector totalu0, totalv0;
- nBodyAssembler_.concatenateVectors(u0, totalu0);
- nBodyAssembler_.concatenateVectors(v0, totalv0);
-
- updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
- updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper); */
-
- //print(vLower, "vLower obstacles:");
- //print(vUpper, "vUpper obstacles:");
-
- //std::cout << "- Problem assembled: success" << std::endl;
-
- //print(ignoreNodes_, "ignoreNodes:");
-
- // set up functional and TNMMG solver
- //using ZeroSolverFactory = SolverFactory<Functional, IgnoreVector>;
- //Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), vLower, vUpper);
- //ZeroSolverFactory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);
- Functional J(bilinearForm, totalRhs, globalFriction_, lower, upper);
- Factory solverFactory(parset_.sub("solver.tnnmg"), J, ignoreNodes_);
- solverFactory.build(linearSolver);
-
- auto step = solverFactory.step();
-
- //std::cout << "- Functional and TNNMG step setup: success" << std::endl;
-
- EnergyNorm<Matrix, Vector> energyNorm(bilinearForm);
- LoopSolver<Vector> velocityProblemSolver(*step.get(), velocityMaxIterations_,
- velocityTolerance_, energyNorm,
- verbosity_);
+ NonlinearSolver<std::remove_reference_t<decltype(*functional)>, IgnoreVector> solver(parset_.sub("solver.tnnmg"), linearSolver, functional, ignoreNodes_);
size_t fixedPointIteration;
size_t multigridIterations = 0;
std::vector<ScalarVector> alpha(nBodies);
updaters.state_->extractAlpha(alpha);
- Vector totalVelocityIterate;
+ Vector totalVelocityIterate, old_v;
nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
-
- // project in onto admissible set
- const size_t blocksize = Vector::block_type::dimension;
- for (size_t i=0; i<totalVelocityIterate.size(); i++) {
- for (size_t j=0; j<blocksize; j++) {
- if (totalVelocityIterate[i][j] < lower[i][j]) {
- totalVelocityIterate[i][j] = lower[i][j];
- }
-
- if (totalVelocityIterate[i][j] > upper[i][j]) {
- totalVelocityIterate[i][j] = upper[i][j];
- }
- }
- }
-
- Vector old_v = totalVelocityIterate;
+ old_v = totalVelocityIterate;
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) {
@@ -193,25 +81,9 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
// contribution from nonlinearity
globalFriction_.updateAlpha(alpha);
- //print(velocityIterates, "velocityIterates:");
- //print(totalVelocityIterate, "totalVelocityIterate:");
- //std::cout << "- FixedPointIteration iterate" << std::endl;
-
- // solve a velocity problem
- solverFactory.setProblem(totalVelocityIterate);
+ auto res = solver.solve(solverParset_, totalVelocityIterate);
- //std::cout << "- Velocity problem set" << std::endl;
-
- velocityProblemSolver.preprocess();
- //std::cout << "-- Preprocessed" << std::endl;
- velocityProblemSolver.solve();
- //std::cout << "-- Solved" << std::endl;
-
- const auto& tnnmgSol = step->getSol();
-
- //std::cout << "FixPointIterator: Energy of TNNMG solution: " << J(tnnmgSol) << std::endl;
-
- nBodyAssembler_.postprocess(tnnmgSol, velocityIterates);
+ nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//print(totalVelocityIterate, "totalVelocityIterate:");
@@ -219,11 +91,11 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
//DUNE_THROW(Dune::Exception, "Just need to stop here!");
- multigridIterations += velocityProblemSolver.getResult().iterations;
+ multigridIterations += res.iterations;
Vector v_m = old_v;
v_m *= 1.0 - lambda_;
- Dune::MatrixVector::addProduct(v_m, lambda_, tnnmgSol);
+ Dune::MatrixVector::addProduct(v_m, lambda_, totalVelocityIterate);
// extract relative velocities in mortar basis
std::vector<Vector> v_rel;
@@ -259,6 +131,8 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
}
}
+ //std::cout << fixedPointIteration << std::endl;
+
if (lambda_ < 1e-12 or breakCriterion) {
//std::cout << "-FixedPointIteration finished! " << (lambda_ < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
fixedPointIteration++;
diff --git a/dune/tectonic/spatial-solving/fixedpointiterator.hh b/dune/tectonic/spatial-solving/fixedpointiterator.hh
index 183f2f8d..3742ed98 100644
--- a/dune/tectonic/spatial-solving/fixedpointiterator.hh
+++ b/dune/tectonic/spatial-solving/fixedpointiterator.hh
@@ -13,8 +13,6 @@
#include <dune/contact/assemblers/nbodyassembler.hh>
-#include "tnnmg/functional.hh"
-#include "tnnmg/zerononlinearity.hh"
#include "solverfactory.hh"
struct FixedPointIterationCounter {
@@ -33,7 +31,6 @@ class FixedPointIterator {
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
- using Functional = typename Factory::Functional; //Functional<Matrix&, Vector&, &, Vector&, Vector&, double>; //;
using Nonlinearity = typename Factory::Functional::Nonlinearity;
const static int dims = Vector::block_type::dimension;
@@ -76,9 +73,8 @@ class FixedPointIterator {
size_t fixedPointMaxIterations_;
double fixedPointTolerance_;
double lambda_;
- size_t velocityMaxIterations_;
- double velocityTolerance_;
- Solver::VerbosityMode verbosity_;
+ const Dune::ParameterTree& solverParset_;
+
const ErrorNorms& errorNorms_;
};
diff --git a/src/foam/foam-2D.cfg b/src/foam/foam-2D.cfg
index d911706c..6f7c2ecc 100644
--- a/src/foam/foam-2D.cfg
+++ b/src/foam/foam-2D.cfg
@@ -1,12 +1,12 @@
# -*- mode:conf -*-
[body0]
-smallestDiameter = 0.05 # 2e-3 [m]
+smallestDiameter = 0.5 # 2e-3 [m]
[body1]
-smallestDiameter = 0.01 # 2e-3 [m]
+smallestDiameter = 0.5 # 2e-3 [m]
[timeSteps]
-refinementTolerance = 1e-5 # 1e-5
+refinementTolerance = 1e-3 # 1e-5
[u0.solver]
tolerance = 1e-8
@@ -18,7 +18,7 @@ tolerance = 1e-8
tolerance = 1e-8
[v.fpi]
-tolerance = 1e-5
+tolerance = 1e-3
[solver.tnnmg.preconditioner.basesolver]
tolerance = 1e-10
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index 59660df5..324bcd88 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -69,8 +69,9 @@
//#include <dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh>
#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
#include <dune/tectonic/spatial-solving/solverfactory.hh>
+#include <dune/tectonic/spatial-solving/makelinearsolver.hh>
-#include <dune/tectonic/time-stepping/adaptivetimestepper.hh>
+#include <dune/tectonic/time-stepping/uniformtimestepper.hh>
#include <dune/tectonic/time-stepping/rate.hh>
#include <dune/tectonic/time-stepping/state.hh>
#include <dune/tectonic/time-stepping/updaters.hh>
@@ -79,6 +80,13 @@
#include <dune/tectonic/utils/diameter.hh>
#include <dune/tectonic/utils/geocoordinate.hh>
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include "dune/tectonic/spatial-solving/tnnmg/zerononlinearity.hh"
+
+
+
// for getcwd
#include <unistd.h>
@@ -188,6 +196,172 @@ int main(int argc, char *argv[]) {
programState.setupInitialConditions(parset, contactNetwork);
}
+ /*
+ // setup initial conditions in program state
+ programState.alpha[0] = 0;
+ programState.alpha[1] = parset.get<double>("boundary.friction.initialAlpha"); */
+
+ /*auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
+ frictionalBoundary, frictionData);
+ myGlobalNonlinearity->updateLogState(alpha_initial);
+*/
+
+ /*
+ // Solving a linear problem with a multigrid solver
+ auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+
+ auto const solveLinearProblem = [&](
+ const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+ const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
+ Dune::ParameterTree const &_localParset) {
+
+ std::vector<const Matrix*> matrices_ptr(_matrices.size());
+ for (size_t i=0; i<matrices_ptr.size(); i++) {
+ matrices_ptr[i] = _matrices[i].get();
+ }
+
+ // assemble full global contact problem
+ Matrix bilinearForm;
+ nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+ Vector totalRhs, oldTotalRhs;
+ nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
+ oldTotalRhs = totalRhs;
+
+ Vector totalX, oldTotalX;
+ nBodyAssembler.nodalToTransformed(_x, totalX);
+ oldTotalX = totalX;
+
+ // get lower and upper obstacles
+ const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+ Vector lower(totalObstacles.size());
+ Vector upper(totalObstacles.size());
+
+ for (size_t j=0; j<totalObstacles.size(); ++j) {
+ const auto& totalObstaclesj = totalObstacles[j];
+ auto& lowerj = lower[j];
+ auto& upperj = upper[j];
+ for (size_t d=0; d<dims; ++d) {
+ lowerj[d] = totalObstaclesj[d][0];
+ upperj[d] = totalObstaclesj[d][1];
+ }
+ }
+
+ // set up functional
+ using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+ Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper);
+
+ // set up TNMMG solver
+ using Factory = SolverFactory<Functional, BitVector>;
+ //Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
+ Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
+ factory.build(linearSolver); */
+
+ /* std::vector<BitVector> bodyDirichletNodes;
+ nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
+ for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
+ print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
+ }*/
+
+ /*
+ auto tnnmgStep = factory.step();
+ factory.setProblem(totalX);
+
+ const EnergyNorm<Matrix, Vector> norm(bilinearForm);
+
+ LoopSolver<Vector> solver(
+ *tnnmgStep.get(), _localParset.get<size_t>("maximumIterations"),
+ _localParset.get<double>("tolerance"), norm,
+ _localParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+ solver.preprocess();
+ solver.solve();
+
+ nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
+ }; */
+
+ /*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> _norm,
+ Dune::ParameterTree const &_localParset) {
+
+
+
+ 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
+ programState.u[0] = 0.0;
+ programState.u[1] = 0.0;
+
+ solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, programState.u,
+ parset.sub("u0.solver"));
+
+ programState.v[0] = 0.0;
+ programState.v[1] = 0.0;
+
+
+ programState.a[0] = 0.0;
+ programState.a[1] = 0.0;
+ {
+ // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+ // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+ std::vector<Vector> accelerationRHS = ell0;
+
+ for (size_t i=0; i<bodyCount; i++) {
+ const auto& body = contactNetwork.body(i);
+ Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
+ }
+
+ BitVector noNodes(dirichletNodes.size(), false);
+ solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, programState.a,
+ parset.sub("a0.solver"));
+
+ }
+
+ for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
+ const auto& coupling = contactNetwork.coupling(i);
+ const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
+
+ const auto nonmortarIdx = coupling->gridIdx_[0];
+ const auto& body = contactNetwork.body(nonmortarIdx);
+
+ ScalarVector frictionBoundaryStress(weightedNormalStress[nonmortarIdx].size());
+
+ body->assembler()->assembleWeightedNormalStress(
+ contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
+ body->data()->getPoissonRatio(), u[nonmortarIdx]);
+
+ weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+ }
+
+
+*/
+
+
+
auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
contactNetwork.body(i)->setDeformation(programState.u[i]);
@@ -198,6 +372,8 @@ int main(int argc, char *argv[]) {
// ------------------------
// assemble global friction
// ------------------------
+ print(programState.weightedNormalStress, "weightedNormalStress");
+
contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
auto& globalFriction = contactNetwork.globalFriction();
@@ -279,43 +455,6 @@ int main(int argc, char *argv[]) {
const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
- auto const mustRefine = [&](Updaters &coarseUpdater,
- Updaters &fineUpdater) {
-
- //return false;
- //std::cout << "Step 1" << std::endl;
-
- std::vector<ScalarVector> coarseAlpha;
- coarseAlpha.resize(bodyCount);
- coarseUpdater.state_->extractAlpha(coarseAlpha);
-
- //print(coarseAlpha, "coarseAlpha:");
-
- std::vector<ScalarVector> fineAlpha;
- fineAlpha.resize(bodyCount);
- fineUpdater.state_->extractAlpha(fineAlpha);
-
- //print(fineAlpha, "fineAlpha:");
-
- //std::cout << "Step 3" << std::endl;
-
- ScalarVector::field_type energyNorm = 0;
- for (size_t i=0; i<stateEnergyNorms.size(); i++) {
- //std::cout << "for " << i << std::endl;
-
- //std::cout << not stateEnergyNorms[i] << std::endl;
-
- if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
- continue;
-
- energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
- }
- //std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance << std::endl;
- //std::cout << "must refine: " << (energyNorm > refinementTolerance) << std::endl;
- return energyNorm > refinementTolerance;
- };
-
-
std::signal(SIGXCPU, handleSignal);
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
@@ -348,21 +487,20 @@ int main(int argc, char *argv[]) {
stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
externalForces, stateEnergyNorms);
- AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
- adaptiveTimeStepper(stepBase, contactNetwork, current,
- programState.relativeTime, programState.relativeTau,
- mustRefine);
+ UniformTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+ timeStepper(stepBase, contactNetwork, current,
+ programState.relativeTime, programState.relativeTau);
size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
- while (!adaptiveTimeStepper.reachedEnd()) {
+ while (!timeStepper.reachedEnd()) {
programState.timeStep++;
//preconditioner.build();
- iterationCount = adaptiveTimeStepper.advance();
+ iterationCount = timeStepper.advance();
- programState.relativeTime = adaptiveTimeStepper.relativeTime_;
- programState.relativeTau = adaptiveTimeStepper.relativeTau_;
+ programState.relativeTime = timeStepper.relativeTime_;
+ programState.relativeTau = timeStepper.relativeTau_;
current.rate_->extractDisplacement(programState.u);
current.rate_->extractVelocity(programState.v);
current.rate_->extractAcceleration(programState.a);
diff --git a/src/foam/foam.cfg b/src/foam/foam.cfg
index b5f4216c..cf375466 100644
--- a/src/foam/foam.cfg
+++ b/src/foam/foam.cfg
@@ -1,77 +1,79 @@
# -*- mode:conf -*-
+outPath = test # output written to ./output/outPath
+
gravity = 9.81 # [m/s^2]
[body0]
-length = 0.4 # [m]
-height = 0.04 # [m]
-depth = 0.04 # [m]
-bulkModulus = 1.5e5 # [Pa] #2190
-poissonRatio = 0.11 # [1] #0.11
+length = 5.0 # [m]
+height = 1.0 # [m]
+bulkModulus = 4.12e7 # [Pa] #2190
+poissonRatio = 0.3 # [1] #0.11
[body0.elastic]
-density = 1300 # [kg/m^3] #750
+density = 5e3 # [kg/m^3] #750
shearViscosity = 0 # [Pas]
bulkViscosity = 0 # [Pas]
[body0.viscoelastic]
-density = 1300 # [kg/m^3]
+density = 5e3 # [kg/m^3] #750
shearViscosity = 0 # [Pas]
bulkViscosity = 0 # [Pas]
[body1]
-length = 0.04 # [m]
-height = 0.04 # [m]
-depth = 0.04 # [m]
-bulkModulus = 1.5e5 # [Pa]
-poissonRatio = 0.11 # [1]
+length = 5.00 # [m]
+height = 1.00 # [m]
+bulkModulus = 4.12e7 # [Pa]
+poissonRatio = 0.3 # [1]
[body1.elastic]
-density = 1300 # [kg/m^3]
+density = 5e3 # [kg/m^3]
shearViscosity = 0 # [Pas]
bulkViscosity = 0 # [Pas]
[body1.viscoelastic]
-density = 1300 # [kg/m^3]
+density = 5e3 # [kg/m^3]
shearViscosity = 0 # [Pas]
bulkViscosity = 0 # [Pas]
+
[boundary.friction]
-C = 6 # [Pa]
-mu0 = 0.48 # [ ]
-V0 = 1e-3 # [m/s]
-L = 1e-6 # [m]
-initialAlpha = 0 # [ ]
+C = 0.0 # [Pa]
+mu0 = 0.6 # [ ]
+V0 = 1e-6 # [m/s]
+L = 1e-5 # [m]
+initialAlpha = -10.0 # [ ]
stateModel = AgeingLaw
frictionModel = Truncated #Regularised
[boundary.friction.weakening]
-a = 0.054 # [ ]
-b = 0.074 # [ ]
+a = 0.010 # [ ]
+b = 0.015 # [ ]
[boundary.friction.strengthening]
-a = 0.054 # [ ]
-b = 0.074 # [ ]
+a = 0.010 # [ ]
+b = 0.015 # [ ]
+
[boundary.neumann]
-sigmaN = 10000.0 # [Pa]
+sigmaN = 0.0 # [Pa]
[boundary.dirichlet]
-finalVelocity = 3e-3 # [m/s]
+finalVelocity = 2e-4 # [m/s]
[io]
-data.write = false
+data.write = true
printProgress = true
restarts.first = 0
restarts.spacing= 20
-restarts.write = false #true
+restarts.write = true #true
vtk.write = true
[problem]
-finalTime = 1000 # [s] #1000
+finalTime = 15 # [s] #1000
bodyCount = 2
[initialTime]
timeStep = 0
relativeTime = 0.0
-relativeTau = 5e-6 # 1e-6
+relativeTau = 1e-4 # 1e-6
[timeSteps]
scheme = newmark
-timeSteps = 1
+timeSteps = 1e4
[u0.solver]
maximumIterations = 100
diff --git a/src/strikeslip/strikeslip-2D.cfg b/src/strikeslip/strikeslip-2D.cfg
index cabe495a..9fd7ce5c 100644
--- a/src/strikeslip/strikeslip-2D.cfg
+++ b/src/strikeslip/strikeslip-2D.cfg
@@ -1,9 +1,9 @@
# -*- mode:conf -*-
[body0]
-smallestDiameter = 0.006 # 2e-3 [m]
+smallestDiameter = 0.02 # 2e-3 [m]
[body1]
-smallestDiameter = 0.006 # 2e-3 [m]
+smallestDiameter = 0.02 # 2e-3 [m]
[timeSteps]
refinementTolerance = 1e-3# 1e-5
diff --git a/src/strikeslip/strikeslip.cfg b/src/strikeslip/strikeslip.cfg
index dc73ee1e..dc92ad1a 100644
--- a/src/strikeslip/strikeslip.cfg
+++ b/src/strikeslip/strikeslip.cfg
@@ -1,21 +1,21 @@
outPath = test # output written to ./output/strikeslip/outPath
# -*- mode:conf -*-
-gravity = 9.81 # [m/s^2]
+gravity = 0.0 # [m/s^2]
[body0]
length = 0.5 # [m]
height = 0.5 # [m]
depth = 0.12 # [m]
-bulkModulus = 1.5e5 # [Pa] #2190
-poissonRatio = 0.11 # [1] #0.11
+bulkModulus = 0.5e5 # [Pa] #2190
+poissonRatio = 0.3 # [1] #0.11
[body0.elastic]
distFromDiag = 0.2
-density = 1300 # [kg/m^3] #750
-shearViscosity = 0 # [Pas]
-bulkViscosity = 0 # [Pas]
+density = 900 # [kg/m^3] #750
+shearViscosity = 1e3 # [Pas]
+bulkViscosity = 1e3 # [Pas]
[body0.viscoelastic]
-density = 1300 # [kg/m^3]
+density = 1000 # [kg/m^3]
shearViscosity = 1e4 # [Pas]
bulkViscosity = 1e4 # [Pas]
@@ -23,32 +23,32 @@ bulkViscosity = 1e4 # [Pas]
length = 0.5 # [m]
height = 0.5 # [m]
depth = 0.12 # [m]
-bulkModulus = 1.5e5 # [Pa]
-poissonRatio = 0.11 # [1]
+bulkModulus = 1e5 # [Pa]
+poissonRatio = 0.0 # [1]
[body1.elastic]
distFromDiag = 0.2
-density = 1300 # [kg/m^3]
-shearViscosity = 0 # [Pas]
-bulkViscosity = 0 # [Pas]
+density = 900 # [kg/m^3]
+shearViscosity = 0.0 # [Pas]
+bulkViscosity = 0.0 # [Pas]
[body1.viscoelastic]
-density = 1300 # [kg/m^3]
-shearViscosity = 1e4 # [Pas]
-bulkViscosity = 1e4 # [Pas]
+density = 1000 # [kg/m^3]
+shearViscosity = 0.0 # [Pas]
+bulkViscosity = 0.0 # [Pas]
[boundary.friction]
-C = 6 # [Pa]
-mu0 = 0.48 # [ ]
-V0 = 1e-3 # [m/s]
-L = 1e-6 # [m]
+C = 10 # [Pa]
+mu0 = 0.7 # [ ]
+V0 = 5e-5 # [m/s]
+L = 2.25e-5 # [m]
initialAlpha = 0 # [ ]
stateModel = AgeingLaw
frictionModel = Truncated #Regularised
[boundary.friction.weakening]
-a = 0.054 # [ ]
-b = 0.074 # [ ]
+a = 0.002 # [ ]
+b = 0.017 # [ ]
[boundary.friction.strengthening]
-a = 0.054 # [ ]
-b = 0.074 # [ ]
+a = 0.020 # [ ]
+b = 0.005 # [ ]
[boundary.neumann]
sigmaN = 200.0 # [Pa]
@@ -60,7 +60,7 @@ finalVelocity = 1e-4 # [m/s]
data.write = true
printProgress = true
restarts.first = 0
-restarts.spacing= 1
+restarts.spacing= 50
restarts.write = true #true
vtk.write = true
@@ -75,7 +75,7 @@ relativeTau = 2e-2 # 1e-6
[timeSteps]
scheme = newmark
-timeSteps = 15
+timeSteps = 10
[u0.solver]
maximumIterations = 100
--
GitLab