Skip to content
Snippets Groups Projects
Commit 7db7854f authored by podlesny's avatar podlesny
Browse files

work in progress: not working remainder

parent 0c08f244
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 290 additions and 424 deletions
dune/tectonic/data-structures/body/body.cc
+ 8
− 0
View file @ 7db7854f
......@@ -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
......
dune/tectonic/data-structures/program_state.hh
+ 15
− 167
View file @ 7db7854f
......@@ -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;
Vector totalX;
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
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");
......
dune/tectonic/spatial-solving/fixedpointiterator.cc
+ 18
− 144
View file @ 7db7854f
......@@ -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);
//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;
FunctionalFactory<std::decay_t<decltype(nBodyAssembler_)>, decltype(globalFriction_), Matrix, Vector> functionalFactory(nBodyAssembler_, globalFriction_, velocityMatrices, velocityRHSs);
functionalFactory.build();
auto functional = functionalFactory.functional();
//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++;
......
dune/tectonic/spatial-solving/fixedpointiterator.hh
+ 2
− 6
View file @ 7db7854f
......@@ -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_;
};
......
src/foam/foam-2D.cfg
+ 4
− 4
View file @ 7db7854f
# -*- 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
......
src/foam/foam.cc
+ 184
− 46
View file @ 7db7854f
......@@ -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);
......
src/foam/foam.cfg
+ 32
− 30
View file @ 7db7854f
# -*- 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
......
src/strikeslip/strikeslip-2D.cfg
+ 2
− 2
View file @ 7db7854f
# -*- 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
......
src/strikeslip/strikeslip.cfg
+ 25
− 25
View file @ 7db7854f
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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment