Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • podlesny/dune-tectonic
  • agnumpde/dune-tectonic
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 884 additions and 399 deletions
src/sand-wedge-data/twopiece.hh → src/one-body-problem-data/segmented-function.hh
View file @ 8fe950bc
#ifndef SRC_SAND_WEDGE_DATA_TWOPIECE_HH #ifndef SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
#define SRC_SAND_WEDGE_DATA_TWOPIECE_HH #define SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
#include <dune/common/function.hh> #include <dune/common/function.hh>
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
...@@ -7,28 +7,24 @@ ...@@ -7,28 +7,24 @@
#include "mygeometry.hh" #include "mygeometry.hh"
class TwoPieceFunction class SegmentedFunction
: public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>, : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
Dune::FieldVector<double, 1>> { Dune::FieldVector<double, 1>> {
private: private:
bool liesBelow(Dune::FieldVector<double, MY_DIM> const &x, bool liesBelow(Dune::FieldVector<double, MY_DIM> const &x,
Dune::FieldVector<double, MY_DIM> const &y, Dune::FieldVector<double, MY_DIM> const &y,
Dune::FieldVector<double, MY_DIM> const &z) const { Dune::FieldVector<double, MY_DIM> const &z) const {
return (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1] - x[1]; return x[1] + (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1];
}; };
bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const { bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
return liesBelow(_K, _M, z); return liesBelow(MyGeometry::K, MyGeometry::M, z);
}; };
Dune::FieldVector<double, MY_DIM> const &_K;
Dune::FieldVector<double, MY_DIM> const &_M;
double const _v1; double const _v1;
double const _v2; double const _v2;
public: public:
TwoPieceFunction(double v1, double v2) SegmentedFunction(double v1, double v2) : _v1(v1), _v2(v2) {}
: _K(MyGeometry::K), _M(MyGeometry::M), _v1(v1), _v2(v2) {}
void evaluate(Dune::FieldVector<double, MY_DIM> const &x, void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
Dune::FieldVector<double, 1> &y) const { Dune::FieldVector<double, 1> &y) const {
......
src/sand-wedge-data/weakpatch.hh → src/one-body-problem-data/weakpatch.hh
View file @ 8fe950bc
#ifndef SRC_SAND_WEDGE_DATA_WEAKPATCH_HH #ifndef SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
#define SRC_SAND_WEDGE_DATA_WEAKPATCH_HH #define SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
template <class LocalVector> ConvexPolyhedron<LocalVector> getWeakPatch() { template <class LocalVector>
ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset) {
ConvexPolyhedron<LocalVector> weakPatch; ConvexPolyhedron<LocalVector> weakPatch;
#if MY_DIM == 3 #if MY_DIM == 3
weakPatch.vertices.resize(4); weakPatch.vertices.resize(4);
...@@ -11,6 +12,16 @@ template <class LocalVector> ConvexPolyhedron<LocalVector> getWeakPatch() { ...@@ -11,6 +12,16 @@ template <class LocalVector> ConvexPolyhedron<LocalVector> getWeakPatch() {
weakPatch.vertices[k][2] = -MyGeometry::depth / 2.0; weakPatch.vertices[k][2] = -MyGeometry::depth / 2.0;
weakPatch.vertices[k + 2][2] = MyGeometry::depth / 2.0; weakPatch.vertices[k + 2][2] = MyGeometry::depth / 2.0;
} }
switch (parset.get<Config::PatchType>("patchType")) {
case Config::Rectangular:
break;
case Config::Trapezoidal:
weakPatch.vertices[1][0] += 0.05 * MyGeometry::lengthScale;
weakPatch.vertices[3][0] -= 0.05 * MyGeometry::lengthScale;
break;
default:
assert(false);
}
#else #else
weakPatch.vertices.resize(2); weakPatch.vertices.resize(2);
weakPatch.vertices[0] = MyGeometry::X; weakPatch.vertices[0] = MyGeometry::X;
......
src/sand-wedge.cc → src/one-body-problem.cc
View file @ 8fe950bc
#include <Python.h>
#ifdef HAVE_CONFIG_H #ifdef HAVE_CONFIG_H
#include "config.h" #include "config.h"
#endif #endif
#ifndef HAVE_PYTHON
#error Python is required
#endif
#ifdef HAVE_IPOPT #ifdef HAVE_IPOPT
#undef HAVE_IPOPT #undef HAVE_IPOPT
#endif #endif
#ifndef datadir #include <atomic>
#error datadir unset
#endif
#include <cmath> #include <cmath>
#include <csignal>
#include <exception> #include <exception>
#include <fstream> #include <fstream>
#include <iostream> #include <iostream>
...@@ -27,6 +19,7 @@ ...@@ -27,6 +19,7 @@
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/function.hh> #include <dune/common/function.hh>
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh> #include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh> #include <dune/common/parametertreeparser.hh>
...@@ -34,61 +27,59 @@ ...@@ -34,61 +27,59 @@
#include <dune/istl/bcrsmatrix.hh> #include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh> #include <dune/istl/bvector.hh>
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wdeprecated-declarations"
#include <dune/fufem/boundarypatch.hh> #include <dune/fufem/boundarypatch.hh>
#pragma clang diagnostic pop #include <dune/fufem/formatstring.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/sharedpointermap.hh>
#include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/sumnorm.hh>
#include <dune/solvers/norms/twonorm.hh>
#include <dune/solvers/solvers/loopsolver.hh> #include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/solvers/solver.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/tnnmg/problem-classes/convexproblem.hh>
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/myblockproblem.hh> #include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/globalfriction.hh> #include <dune/tectonic/globalfriction.hh>
#include <dune/fufem/hdf5/file.hh>
#include "adaptivetimestepper.hh" // FIXME
#include "assemblers.hh" #include "assemblers.hh"
#include "tobool.hh" #include "diameter.hh"
#include "coupledtimestepper.hh"
#include "distances.hh"
#include "enumparser.hh" #include "enumparser.hh"
#include "enums.hh" #include "enums.hh"
#include "fixedpointiterator.hh"
#include "friction_writer.hh"
#include "gridselector.hh" #include "gridselector.hh"
#include "sand-wedge-data/mybody.hh" #include "hdf5-writer.hh"
#include "sand-wedge-data/mygeometry.hh" #include "hdf5/restart-io.hh"
#include "sand-wedge-data/myglobalfrictiondata.hh" #include "matrices.hh"
#include "sand-wedge-data/mygrid.hh" #include "program_state.hh"
#include "sand-wedge-data/special_writer.hh" #include "one-body-problem-data/bc.hh"
#include "sand-wedge-data/weakpatch.hh" #include "one-body-problem-data/mybody.hh"
#include "solverfactory.hh" #include "one-body-problem-data/mygeometry.hh"
#include "state.hh" #include "one-body-problem-data/myglobalfrictiondata.hh"
#include "timestepping.hh" #include "one-body-problem-data/mygrid.hh"
#include "one-body-problem-data/weakpatch.hh"
#include "spatial-solving/solverfactory.hh"
#include "time-stepping/adaptivetimestepper.hh"
#include "time-stepping/rate.hh"
#include "time-stepping/state.hh"
#include "time-stepping/updaters.hh"
#include "vtk.hh" #include "vtk.hh"
size_t const dims = MY_DIM; size_t const dims = MY_DIM;
void initPython() { Dune::ParameterTree getParameters(int argc, char *argv[]) {
Python::start(); Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
Python::run("import sys"); Dune::ParameterTreeParser::readINITree(
Python::run("sys.path.append('" datadir "')"); Dune::Fufem::formatString("one-body-problem-%dD.cfg", dims), parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
return parset;
} }
static std::atomic<bool> terminationRequested(false);
void handleSignal(int signum) { terminationRequested = true; }
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
try { try {
Dune::ParameterTree parset; Dune::MPIHelper::instance(argc, argv);
Dune::ParameterTreeParser::readINITree(datadir "/parset.cfg", parset); auto const parset = getParameters(argc, argv);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
MyGeometry::render(); MyGeometry::render();
MyGeometry::write(); MyGeometry::write();
...@@ -96,14 +87,13 @@ int main(int argc, char *argv[]) { ...@@ -96,14 +87,13 @@ int main(int argc, char *argv[]) {
using GridView = Grid::LeafGridView; using GridView = Grid::LeafGridView;
using MyAssembler = MyAssembler<GridView, dims>; using MyAssembler = MyAssembler<GridView, dims>;
using Matrix = MyAssembler::Matrix; using Matrix = MyAssembler::Matrix;
using LocalMatrix = Matrix::block_type;
using Vector = MyAssembler::Vector; using Vector = MyAssembler::Vector;
using LocalVector = Vector::block_type; using LocalVector = Vector::block_type;
using ScalarMatrix = MyAssembler::ScalarMatrix; using ScalarMatrix = MyAssembler::ScalarMatrix;
using ScalarVector = MyAssembler::ScalarVector; using ScalarVector = MyAssembler::ScalarVector;
using LocalScalarVector = ScalarVector::block_type;
auto weakPatch = getWeakPatch<LocalVector>(); auto const weakPatch =
getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"));
// {{{ Set up grid // {{{ Set up grid
GridConstructor<Grid> gridConstructor; GridConstructor<Grid> gridConstructor;
...@@ -112,13 +102,10 @@ int main(int argc, char *argv[]) { ...@@ -112,13 +102,10 @@ int main(int argc, char *argv[]) {
refine(*grid, weakPatch, refine(*grid, weakPatch,
parset.get<double>("boundary.friction.smallestDiameter")); parset.get<double>("boundary.friction.smallestDiameter"));
auto const refinements = grid->maxLevel();
double minDiameter = std::numeric_limits<double>::infinity(); double minDiameter = std::numeric_limits<double>::infinity();
double maxDiameter = 0.0; double maxDiameter = 0.0;
for (auto it = grid->template leafbegin<0>(); for (auto &&e : elements(grid->leafGridView())) {
it != grid->template leafend<0>(); ++it) { auto const geometry = e.geometry();
auto const geometry = it->geometry();
auto const diam = diameter(geometry); auto const diam = diameter(geometry);
minDiameter = std::min(minDiameter, diam); minDiameter = std::min(minDiameter, diam);
maxDiameter = std::max(maxDiameter, diam); maxDiameter = std::max(maxDiameter, diam);
...@@ -126,25 +113,16 @@ int main(int argc, char *argv[]) { ...@@ -126,25 +113,16 @@ int main(int argc, char *argv[]) {
std::cout << "min diameter: " << minDiameter << std::endl; std::cout << "min diameter: " << minDiameter << std::endl;
std::cout << "max diameter: " << maxDiameter << std::endl; std::cout << "max diameter: " << maxDiameter << std::endl;
GridView const leafView = grid->leafGridView(); auto const leafView = grid->leafGridView();
size_t const leafVertexCount = leafView.size(dims); auto const leafVertexCount = leafView.size(dims);
std::cout << "Number of DOFs: " << leafVertexCount << std::endl; std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
auto myFaces = gridConstructor.constructFaces(leafView); auto myFaces = gridConstructor.constructFaces(leafView);
// Neumann boundary
BoundaryPatch<GridView> const neumannBoundary(leafView); BoundaryPatch<GridView> const neumannBoundary(leafView);
// Frictional Boundary
BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower; BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
Dune::BitSetVector<1> frictionalNodes(leafVertexCount);
frictionalBoundary.getVertices(frictionalNodes);
// Surface
BoundaryPatch<GridView> const &surface = myFaces.upper; BoundaryPatch<GridView> const &surface = myFaces.upper;
Dune::BitSetVector<1> surfaceNodes(leafVertexCount);
surface.getVertices(surfaceNodes);
// Dirichlet Boundary // Dirichlet Boundary
Dune::BitSetVector<dims> noNodes(leafVertexCount); Dune::BitSetVector<dims> noNodes(leafVertexCount);
...@@ -164,257 +142,188 @@ int main(int argc, char *argv[]) { ...@@ -164,257 +142,188 @@ int main(int argc, char *argv[]) {
// Set up functions for time-dependent boundary conditions // Set up functions for time-dependent boundary conditions
using Function = Dune::VirtualFunction<double, double>; using Function = Dune::VirtualFunction<double, double>;
using FunctionMap = SharedPointerMap<std::string, Function>; Function const &velocityDirichletFunction = VelocityDirichletCondition();
FunctionMap functions; Function const &neumannFunction = NeumannCondition();
{
initPython();
Python::import("boundaryconditions")
.get("Functions")
.toC<typename FunctionMap::Base>(functions);
}
auto const &velocityDirichletFunction =
functions.get("velocityDirichletCondition"),
&neumannFunction = functions.get("neumannCondition");
std::fstream timeStepWriter("timeSteps", std::fstream::out);
timeStepWriter << std::fixed << std::setprecision(10);
auto const reportTimeStep = [&](double _relativeTime, double _relativeTau) {
timeStepWriter << _relativeTime << " " << _relativeTau << std::endl;
};
MyAssembler myAssembler(leafView); MyAssembler const myAssembler(leafView);
MyBody<dims> const body(parset); MyBody<dims> const body(parset);
Matrix A, C, M; Matrices<Matrix> matrices;
myAssembler.assembleElasticity(body.getYoungModulus(), myAssembler.assembleElasticity(body.getYoungModulus(),
body.getPoissonRatio(), A); body.getPoissonRatio(), matrices.elasticity);
myAssembler.assembleViscosity(body.getShearViscosityField(), myAssembler.assembleViscosity(body.getShearViscosityField(),
body.getBulkViscosityField(), C); body.getBulkViscosityField(),
myAssembler.assembleMass(body.getDensityField(), M); matrices.damping);
EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M); myAssembler.assembleMass(body.getDensityField(), matrices.mass);
ScalarMatrix frictionalBoundaryMass; ScalarMatrix relativeFrictionalBoundaryMass;
myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary, myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
frictionalBoundaryMass); relativeFrictionalBoundaryMass);
relativeFrictionalBoundaryMass /= frictionalBoundary.area();
EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm( EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
frictionalBoundaryMass); relativeFrictionalBoundaryMass);
// Assemble forces // Assemble forces
Vector gravityFunctional; Vector gravityFunctional;
myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional); myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
// Problem formulation: right-hand side // Problem formulation: right-hand side
std::function<void(double, Vector &)> computeExternalForces = [&]( std::function<void(double, Vector &)> computeExternalForces =
double _relativeTime, Vector &_ell) { [&](double _relativeTime, Vector &_ell) {
myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction, myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
_relativeTime); _relativeTime);
_ell += gravityFunctional; _ell += gravityFunctional;
}; };
Vector ell0(leafVertexCount);
computeExternalForces(0.0, ell0); using MyProgramState = ProgramState<Vector, ScalarVector>;
MyProgramState programState(leafVertexCount);
// {{{ Initial conditions auto const firstRestart = parset.get<size_t>("io.restarts.first");
using LinearFactory = SolverFactory< auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
dims, BlockNonlinearTNNMGProblem<ConvexProblem< auto const writeRestarts = parset.get<bool>("io.restarts.write");
ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>, auto const writeData = parset.get<bool>("io.data.write");
Grid>; bool const handleRestarts = writeRestarts or firstRestart > 0;
ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
auto dataFile =
auto const solveLinearProblem = [&]( writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix, auto restartFile = handleRestarts
Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> const &_norm, ? std::make_unique<HDF5::File>(
Dune::ParameterTree const &_localParset) { "restarts.h5",
LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME writeRestarts ? HDF5::Access::READWRITE
refinements, *grid, _dirichletNodes); : HDF5::Access::READONLY)
: nullptr;
typename LinearFactory::ConvexProblem convexProblem( auto restartIO = handleRestarts
1.0, _matrix, zeroNonlinearity, _rhs, _x); ? std::make_unique<RestartIO<MyProgramState>>(
typename LinearFactory::BlockProblem problem(parset, convexProblem); *restartFile, leafVertexCount)
: nullptr;
auto multigridStep = factory.getSolver();
multigridStep->setProblem(_x, problem); if (firstRestart > 0) // automatically adjusts the time and timestep
LoopSolver<Vector> solver( restartIO->read(firstRestart, programState);
multigridStep, _localParset.get<size_t>("maximumIterations"), else
_localParset.get<double>("tolerance"), &_norm, programState.setupInitialConditions(parset, computeExternalForces,
_localParset.get<Solver::VerbosityMode>("verbosity"), matrices, myAssembler, dirichletNodes,
false); // absolute error noNodes, frictionalBoundary, body);
solver.preprocess();
solver.solve();
};
// Solve the stationary problem
Vector u_initial(leafVertexCount);
u_initial = 0.0;
solveLinearProblem(dirichletNodes, A, ell0, u_initial, ANorm,
parset.sub("u0.solver"));
ScalarVector alpha_initial(leafVertexCount);
alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
ScalarVector normalStress(leafVertexCount);
myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
body.getYoungModulus(),
body.getPoissonRatio(), u_initial);
MyGlobalFrictionData<LocalVector> frictionInfo( MyGlobalFrictionData<LocalVector> frictionInfo(
parset.sub("boundary.friction"), weakPatch); parset.sub("boundary.friction"), weakPatch);
auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity( auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
frictionalBoundary, frictionInfo, normalStress); frictionalBoundary, frictionInfo, programState.weightedNormalStress);
myGlobalFriction->updateAlpha(alpha_initial); myGlobalFriction->updateAlpha(programState.alpha);
Vector v_initial(leafVertexCount);
v_initial = 0.0;
{
double v_initial_const;
velocityDirichletFunction.evaluate(0.0, v_initial_const);
assert(std::abs(v_initial_const) < 1e-14);
}
Vector a_initial(leafVertexCount);
a_initial = 0.0;
{
// We solve Ma = ell0 - [Au + Cv + DPsi(v)]
Vector accelerationRHS(leafVertexCount);
{
accelerationRHS = 0.0;
Arithmetic::addProduct(accelerationRHS, A, u_initial);
Arithmetic::addProduct(accelerationRHS, C, v_initial);
// NOTE: We assume differentiability of Psi at 0 here!
myGlobalFriction->addGradient(v_initial, accelerationRHS);
accelerationRHS *= -1.0;
accelerationRHS += ell0;
}
solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
parset.sub("a0.solver"));
}
// }}}
Vector vertexCoordinates(leafVertexCount); Vector vertexCoordinates(leafVertexCount);
{ {
Dune::MultipleCodimMultipleGeomTypeMapper< Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView); GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) { for (auto &&v : vertices(leafView))
auto const geometry = it->geometry(); vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
assert(geometry.corners() == 1);
vertexCoordinates[vertexMapper.index(*it)] = geometry.corner(0);
}
} }
FrictionWriter<ScalarVector, Vector> frictionWriter( using MyVertexBasis = typename MyAssembler::VertexBasis;
vertexCoordinates, frictionalNodes, "friction", auto dataWriter =
MyGeometry::horizontalProjection); writeData ? std::make_unique<
BoundaryWriter<ScalarVector, Vector> verticalSurfaceWriter( HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
vertexCoordinates, surfaceNodes, "verticalSurface", *dataFile, vertexCoordinates, myAssembler.vertexBasis,
MyGeometry::verticalProjection); surface, frictionalBoundary, weakPatch)
BoundaryWriter<ScalarVector, Vector> horizontalSurfaceWriter( : nullptr;
vertexCoordinates, surfaceNodes, "horizontalSurface", MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
MyGeometry::horizontalProjection); myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
auto const report = [&](Vector const &_u, Vector const &_v, IterationRegister iterationCount;
ScalarVector const &_alpha) { auto const report = [&](bool initial = false) {
horizontalSurfaceWriter.writeKinetics(_u, _v); if (writeData) {
verticalSurfaceWriter.writeKinetics(_u, _v); dataWriter->reportSolution(programState, *myGlobalFriction);
if (!initial)
ScalarVector c; dataWriter->reportIterations(programState, iterationCount);
myGlobalFriction->coefficientOfFriction(_v, c); dataFile->flush();
frictionWriter.writeKinetics(_u, _v); }
frictionWriter.writeOther(c, _alpha);
};
report(u_initial, v_initial, alpha_initial);
MyVTKWriter<typename MyAssembler::VertexBasis, if (writeRestarts and !initial and
typename MyAssembler::CellBasis> const programState.timeStep % restartSpacing == 0) {
vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs"); restartIO->write(programState);
restartFile->flush();
}
if (parset.get<bool>("io.writeVTK")) { if (parset.get<bool>("io.printProgress"))
ScalarVector stress; std::cout << "timeStep = " << std::setw(6) << programState.timeStep
myAssembler.assembleVonMisesStress( << ", time = " << std::setw(12) << programState.relativeTime
body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress); << ", tau = " << std::setw(12) << programState.relativeTau
vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress); << std::endl;
}
SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities", if (parset.get<bool>("io.vtk.write")) {
leafView); ScalarVector stress;
SpecialWriter<GridView, dims> specialDisplacementWriter( myAssembler.assembleVonMisesStress(body.getYoungModulus(),
"specialDisplacements", leafView); body.getPoissonRatio(),
programState.u, stress);
vtkWriter.write(programState.timeStep, programState.u, programState.v,
programState.alpha, stress);
}
};
report(true);
// Set up TNNMG solver // Set up TNNMG solver
using NonlinearFactory = SolverFactory< using NonlinearFactory = SolverFactory<
dims, dims,
MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>, MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
Grid>; Grid>;
NonlinearFactory factory(parset.sub("solver.tnnmg"), refinements, *grid, NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes);
dirichletNodes);
using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
using UpdaterPair = std::pair< StateUpdater<ScalarVector, Vector>>;
std::shared_ptr<StateUpdater<ScalarVector, Vector>>, MyUpdater current(
std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dims>>>; initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
UpdaterPair current( velocityDirichletFunction, dirichletNodes, matrices,
programState.u, programState.v, programState.a),
initStateUpdater<ScalarVector, Vector>( initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"), parset.get<Config::stateModel>("boundary.friction.stateModel"),
alpha_initial, frictionalNodes, programState.alpha, *frictionalBoundary.getVertices(),
parset.get<double>("boundary.friction.L"), parset.get<double>("boundary.friction.L"),
parset.get<double>("boundary.friction.V0")), parset.get<double>("boundary.friction.V0")));
initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, dirichletNodes, M, A, C,
u_initial, v_initial, a_initial));
auto const refinementTolerance = auto const refinementTolerance =
parset.get<double>("timeSteps.refinementTolerance"); parset.get<double>("timeSteps.refinementTolerance");
auto const mustRefine = [&](UpdaterPair &coarseUpdater, auto const mustRefine = [&](MyUpdater &coarseUpdater,
UpdaterPair &fineUpdater) { MyUpdater &fineUpdater) {
ScalarVector coarseAlpha; ScalarVector coarseAlpha;
coarseUpdater.first->extractAlpha(coarseAlpha); coarseUpdater.state_->extractAlpha(coarseAlpha);
ScalarVector fineAlpha; ScalarVector fineAlpha;
fineUpdater.first->extractAlpha(fineAlpha); fineUpdater.state_->extractAlpha(fineAlpha);
return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance; return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
}; };
size_t timeStep = 1; std::signal(SIGXCPU, handleSignal);
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
AdaptiveTimeStepper<NonlinearFactory, UpdaterPair, AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
EnergyNorm<ScalarMatrix, ScalarVector>> EnergyNorm<ScalarMatrix, ScalarVector>>
adaptiveTimeStepper(factory, parset, myGlobalFriction, current, adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
computeExternalForces, stateEnergyNorm, mustRefine); programState.relativeTime, programState.relativeTau,
computeExternalForces, stateEnergyNorm, mustRefine);
while (!adaptiveTimeStepper.reachedEnd()) { while (!adaptiveTimeStepper.reachedEnd()) {
adaptiveTimeStepper.advance(); programState.timeStep++;
reportTimeStep(adaptiveTimeStepper.getRelativeTime(),
adaptiveTimeStepper.getRelativeTau());
Vector u, v;
ScalarVector alpha;
current.second->extractDisplacement(u);
current.second->extractVelocity(v);
current.first->extractAlpha(alpha);
report(u, v, alpha);
{
BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
myAssembler.vertexBasis, v);
BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
displacement(myAssembler.vertexBasis, u);
specialVelocityWriter.write(velocity);
specialDisplacementWriter.write(displacement);
}
if (parset.get<bool>("io.writeVTK")) { iterationCount = adaptiveTimeStepper.advance();
ScalarVector stress;
myAssembler.assembleVonMisesStress(body.getYoungModulus(), programState.relativeTime = adaptiveTimeStepper.relativeTime_;
body.getPoissonRatio(), u, stress); programState.relativeTau = adaptiveTimeStepper.relativeTau_;
vtkWriter.write(timeStep, u, v, alpha, stress); current.rate_->extractDisplacement(programState.u);
current.rate_->extractVelocity(programState.v);
current.rate_->extractAcceleration(programState.a);
current.state_->extractAlpha(programState.alpha);
report();
if (terminationRequested) {
std::cerr << "Terminating prematurely" << std::endl;
break;
} }
timeStep++;
} }
timeStepWriter.close(); } catch (Dune::Exception &e) {
Python::stop();
}
catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl; Dune::derr << "Dune reported error: " << e << std::endl;
} } catch (std::exception &e) {
catch (std::exception &e) {
std::cerr << "Standard exception: " << e.what() << std::endl; std::cerr << "Standard exception: " << e.what() << std::endl;
} }
} }
src/sand-wedge-data/parset.cfg → src/one-body-problem.cfg
View file @ 8fe950bc
...@@ -2,11 +2,15 @@ ...@@ -2,11 +2,15 @@
gravity = 9.81 # [m/s^2] gravity = 9.81 # [m/s^2]
[io] [io]
data.write = true
printProgress = false printProgress = false
writeVTK = false restarts.first = 0
restarts.spacing= 20
restarts.write = true
vtk.write = false
[problem] [problem]
finalTime = 1800 # [s] finalTime = 1000 # [s]
[body] [body]
bulkModulus = 0.5e5 # [Pa] bulkModulus = 0.5e5 # [Pa]
...@@ -21,49 +25,41 @@ shearViscosity = 1e4 # [Pas] ...@@ -21,49 +25,41 @@ shearViscosity = 1e4 # [Pas]
bulkViscosity = 1e4 # [Pas] bulkViscosity = 1e4 # [Pas]
[boundary.friction] [boundary.friction]
smallestDiameter= 2e-3 # [m]
C = 10 # [Pa] C = 10 # [Pa]
mu0 = 0.7 # [ ] mu0 = 0.7 # [ ]
V0 = 5e-5 # [m/s] V0 = 5e-5 # [m/s]
L = 2.5e-5# [m] L = 2.25e-5 # [m]
initialAlpha = 0 # [ ] initialAlpha = 0 # [ ]
stateModel = AgeingLaw stateModel = AgeingLaw
frictionModel = Truncated frictionModel = Regularised
[boundary.friction.weakening] [boundary.friction.weakening]
a = 0.002 # [ ] a = 0.002 # [ ]
b = 0.014 # [ ] b = 0.017 # [ ]
[boundary.friction.strengthening] [boundary.friction.strengthening]
a = 0.025 # [ ] a = 0.020 # [ ]
b = 0.005 # [ ] b = 0.005 # [ ]
[timeSteps] [timeSteps]
refinementTolerance = 1e-5
number = 100000
scheme = newmark scheme = newmark
[u0.solver] [u0.solver]
tolerance = 1e-8
maximumIterations = 100000 maximumIterations = 100000
verbosity = quiet verbosity = quiet
[a0.solver] [a0.solver]
tolerance = 1e-8
maximumIterations = 100000 maximumIterations = 100000
verbosity = quiet verbosity = quiet
[v.solver] [v.solver]
tolerance = 1e-8
maximumIterations = 100000 maximumIterations = 100000
verbosity = quiet verbosity = quiet
[v.fpi] [v.fpi]
tolerance = 1e-5
maximumIterations = 10000 maximumIterations = 10000
lambda = 0.5 lambda = 0.5
[solver.tnnmg.linear] [solver.tnnmg.linear]
maxiumumIterations = 100000 maxiumumIterations = 100000
tolerance = 1e-10
pre = 3 pre = 3
cycle = 1 # 1 = V, 2 = W, etc. cycle = 1 # 1 = V, 2 = W, etc.
post = 3 post = 3
......
src/polyhedrondistance.hh deleted 100644 → 0
View file @ 4a97dcf2
#ifndef SRC_POLYHEDRONDISTANCE_HH
#define SRC_POLYHEDRONDISTANCE_HH
// Based on the closest point projection from dune-contact
#include <dune/common/dynmatrix.hh>
#include <dune/solvers/common/arithmetic.hh>
template <class Coordinate> struct ConvexPolyhedron {
std::vector<Coordinate> vertices;
template <class DynamicVector>
Coordinate barycentricToGlobal(DynamicVector const &b) const {
assert(b.size() == vertices.size());
Coordinate r(0);
for (size_t i = 0; i < b.size(); i++)
Arithmetic::addProduct(r, b[i], vertices[i]);
return r;
}
template <class DynamicVector>
void sanityCheck(DynamicVector const &b, double tol) const {
double sum = 0;
for (size_t i = 0; i < b.size(); i++) {
if (b[i] < -tol)
DUNE_THROW(Dune::Exception, "No barycentric coords " << b[1] << " nr. "
<< i);
sum += b[i];
}
if (sum > 1 + tol)
DUNE_THROW(Dune::Exception, "No barycentric coords, sum: " << sum);
}
};
template <class Coordinate, class Matrix>
void populateMatrix(ConvexPolyhedron<Coordinate> const &segment,
Matrix &matrix) {
size_t const nCorners = segment.vertices.size();
for (size_t i = 0; i < nCorners; i++)
for (size_t j = 0; j <= i; j++)
matrix[i][j] = segment.vertices[i] * segment.vertices[j];
// make symmetric
for (size_t i = 0; i < nCorners; i++)
for (size_t j = i + 1; j < nCorners; j++)
matrix[i][j] = matrix[j][i];
}
/**
* The functional to be minimized is
*
* 0.5*norm(\sum_i lambda_i corner_i - target)^2
*
* where lambda_i are barycentric coordinates, i.e.
*
* 0 \le lambda_i \le 1 and \sum_i lambda_i=1
*
* The resulting quadratic minimization problem is given by
*
* 0.5 lambda^T*A*lambda - l^T*lambda
*
* with A_ij = (corner_i,corner_j) and l_i = (corner_i, target)
*/
template <class Coordinate>
void approximate(Coordinate const &target, Dune::DynamicMatrix<double> const &A,
Dune::DynamicVector<double> const &l,
ConvexPolyhedron<Coordinate> const &segment,
Dune::DynamicVector<double> &sol) {
size_t nCorners = segment.vertices.size();
// compute the residual
Dune::DynamicVector<double> rhs = l;
// compute rhs -= A*sol_
for (size_t k = 0; k < nCorners; k++)
for (size_t j = 0; j < nCorners; j++)
rhs[k] -= A[k][j] * sol[j];
// use the edge vectors as search directions
double alpha = 0.0;
for (size_t j1 = 0; j1 < nCorners - 1; ++j1) {
for (size_t j2 = j1 + 1; j2 < nCorners; ++j2) {
// compute matrix entry and rhs for edge direction
double a = A[j1][j1] - A[j1][j2] - A[j2][j1] + A[j2][j2];
double b = rhs[j1] - rhs[j2];
// compute minimizer for correction
if (a > 0) {
alpha = b / a;
double const largestAlpha = sol[j2];
double const smallestAlpha = -sol[j1];
// project alpha such that we stay positive
if (alpha > largestAlpha)
alpha = largestAlpha;
else if (alpha < smallestAlpha)
alpha = smallestAlpha;
} else {
// if the quadratic term vanishes, we have a linear function, which
// attains its extrema on the boundary
double sum = sol[j1] + sol[j2];
double lValue = 0.5 * A[j1][j1] * sum * sum - rhs[j1] * sum;
double uValue = 0.5 * A[j2][j2] * sum * sum - rhs[j2] * sum;
alpha = lValue < uValue ? sol[j2] : -sol[j1];
}
// apply correction
sol[j1] += alpha;
sol[j2] -= alpha;
// update the local residual for corrections
for (size_t p = 0; p < nCorners; p++)
rhs[p] -= (A[p][j1] - A[p][j2]) * alpha;
}
}
}
// Point-to-Polyhedron
template <class Coordinate>
double distance(const Coordinate &target,
const ConvexPolyhedron<Coordinate> &segment,
double correctionTolerance) {
size_t nCorners = segment.vertices.size();
double const barycentricTolerance = 1e-14;
Dune::DynamicMatrix<double> A(nCorners, nCorners);
populateMatrix(segment, A);
Dune::DynamicVector<double> l(nCorners);
for (size_t i = 0; i < nCorners; i++)
l[i] = target * segment.vertices[i];
// choose a feasible initial solution
Dune::DynamicVector<double> sol(nCorners);
for (size_t i = 0; i < nCorners; i++)
sol[i] = 1.0 / nCorners;
Coordinate oldCoordinates = segment.barycentricToGlobal(sol);
Coordinate newCoordinates;
size_t maxIterations = 1000;
size_t iterations;
for (iterations = 0; iterations < maxIterations; ++iterations) {
approximate(target, A, l, segment, sol);
newCoordinates = segment.barycentricToGlobal(sol);
if ((oldCoordinates - newCoordinates).two_norm() < correctionTolerance)
break;
oldCoordinates = newCoordinates;
}
assert(iterations < maxIterations);
segment.sanityCheck(sol, barycentricTolerance);
return (target - newCoordinates).two_norm();
}
// Polyhedron-to-Polyhedron
template <class Coordinate>
double distance(const ConvexPolyhedron<Coordinate> &s1,
const ConvexPolyhedron<Coordinate> &s2,
double valueCorrectionTolerance) {
size_t nCorners1 = s1.vertices.size();
size_t nCorners2 = s2.vertices.size();
double const barycentricTolerance = 1e-14;
// choose feasible initial solutions
Dune::DynamicVector<double> sol1(nCorners1);
for (size_t i = 0; i < s1.vertices.size(); i++)
sol1[i] = 1.0 / s1.vertices.size();
auto c1 = s1.barycentricToGlobal(sol1);
Dune::DynamicVector<double> sol2(nCorners2);
for (size_t i = 0; i < nCorners2; i++)
sol2[i] = 1.0 / nCorners2;
auto c2 = s2.barycentricToGlobal(sol2);
double oldDistance = (c2 - c1).two_norm();
Dune::DynamicMatrix<double> A1(nCorners1, nCorners1);
populateMatrix(s1, A1);
Dune::DynamicMatrix<double> A2(nCorners2, nCorners2);
populateMatrix(s2, A2);
size_t maxIterations = 1000;
size_t iterations;
for (iterations = 0; iterations < maxIterations; ++iterations) {
Dune::DynamicVector<double> l2(nCorners2);
for (size_t i = 0; i < nCorners2; i++)
l2[i] = c1 * s2.vertices[i];
approximate(c1, A2, l2, s2, sol2);
c2 = s2.barycentricToGlobal(sol2);
Dune::DynamicVector<double> l1(nCorners1);
for (size_t i = 0; i < nCorners1; i++)
l1[i] = c2 * s1.vertices[i];
approximate(c2, A1, l1, s1, sol1);
c1 = s1.barycentricToGlobal(sol1);
double const newDistance = (c2 - c1).two_norm();
assert(newDistance - oldDistance <= 1e-12); // debugging
if (oldDistance - newDistance <= valueCorrectionTolerance) {
s1.sanityCheck(sol1, barycentricTolerance);
s2.sanityCheck(sol2, barycentricTolerance);
return newDistance;
}
oldDistance = newDistance;
}
assert(false);
}
#endif
src/program_state.hh 0 → 100644
View file @ 8fe950bc
#ifndef SRC_PROGRAM_STATE_HH
#define SRC_PROGRAM_STATE_HH
#include <dune/common/parametertree.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
#include <dune/tectonic/body.hh>
#include "assemblers.hh"
#include "matrices.hh"
#include "spatial-solving/solverfactory.hh"
template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
public:
using Vector = VectorTEMPLATE;
using ScalarVector = ScalarVectorTEMPLATE;
ProgramState(int leafVertexCount)
: u(leafVertexCount),
v(leafVertexCount),
a(leafVertexCount),
alpha(leafVertexCount),
weightedNormalStress(leafVertexCount) {}
// Set up initial conditions
template <class Matrix, class GridView>
void setupInitialConditions(
Dune::ParameterTree const &parset,
std::function<void(double, Vector &)> externalForces,
Matrices<Matrix> const matrices,
MyAssembler<GridView, Vector::block_type::dimension> const &myAssembler,
Dune::BitSetVector<Vector::block_type::dimension> const &dirichletNodes,
Dune::BitSetVector<Vector::block_type::dimension> const &noNodes,
BoundaryPatch<GridView> const &frictionalBoundary,
Body<Vector::block_type::dimension> const &body) {
using LocalVector = typename Vector::block_type;
using LocalMatrix = typename Matrix::block_type;
auto constexpr dims = LocalVector::dimension;
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
Vector const &_rhs, Vector &_x,
Dune::ParameterTree const &_localParset) {
using LinearFactory = SolverFactory<
dims, BlockNonlinearTNNMGProblem<ConvexProblem<
ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
typename GridView::Grid>;
ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
myAssembler.gridView.grid(), _dirichletNodes);
typename LinearFactory::ConvexProblem convexProblem(
1.0, _matrix, zeroNonlinearity, _rhs, _x);
typename LinearFactory::BlockProblem problem(parset, convexProblem);
auto multigridStep = factory.getStep();
multigridStep->setProblem(_x, problem);
EnergyNorm<Matrix, Vector> const norm(_matrix);
LoopSolver<Vector> solver(
multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
_localParset.get<double>("tolerance"), &norm,
_localParset.get<Solver::VerbosityMode>("verbosity"),
false); // absolute error
solver.preprocess();
solver.solve();
};
timeStep = 0;
relativeTime = 0.0;
relativeTau = 1e-6;
Vector ell0(u.size());
externalForces(relativeTime, ell0);
// Initial velocity
v = 0.0;
// Initial displacement: Start from a situation of minimal stress,
// which is automatically attained in the case [v = 0 = a].
// Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
solveLinearProblem(dirichletNodes, matrices.elasticity, ell0, u,
parset.sub("u0.solver"));
// Initial acceleration: Computed in agreement with Ma = ell0 - Au
// (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
Vector accelerationRHS = ell0;
Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity, u);
solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
parset.sub("a0.solver"));
// Initial state
alpha = parset.get<double>("boundary.friction.initialAlpha");
// Initial normal stress
myAssembler.assembleWeightedNormalStress(
frictionalBoundary, weightedNormalStress, body.getYoungModulus(),
body.getPoissonRatio(), u);
}
public:
Vector u;
Vector v;
Vector a;
ScalarVector alpha;
ScalarVector weightedNormalStress;
double relativeTime;
double relativeTau;
size_t timeStep;
};
#endif
src/sand-wedge-data/boundaryconditions.py deleted 100644 → 0
View file @ 4a97dcf2
import math
class neumannCondition:
def __call__(self, relativeTime):
return 0
class velocityDirichletCondition:
def __call__(self, relativeTime):
finalVelocity = -5e-5
if relativeTime <= 0.1:
return finalVelocity * ( 1.0 - math.cos(relativeTime * math.pi / 0.1) ) / 2.0
return finalVelocity
Functions = {
'neumannCondition' : neumannCondition(),
'velocityDirichletCondition' : velocityDirichletCondition()
}
src/sand-wedge-data/special_writer.hh deleted 100644 → 0
View file @ 4a97dcf2
#ifndef SRC_SPECIAL_WRITER_HH
#define SRC_SPECIAL_WRITER_HH
#include <fstream>
#include <utility>
#include <dune/common/fvector.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include "mygeometry.hh"
template <class GridView, int dimension> class SpecialWriter {
using LocalVector = Dune::FieldVector<double, dimension>;
using Element = typename GridView::Grid::template Codim<0>::Entity;
using ElementPointer =
typename GridView::Grid::template Codim<0>::EntityPointer;
void writeHorizontal(LocalVector const &v) {
writer_ << MyGeometry::horizontalProjection(v) << " ";
}
void writeVertical(LocalVector const &v) {
writer_ << MyGeometry::verticalProjection(v) << " ";
}
std::pair<ElementPointer, LocalVector> globalToLocal(LocalVector const &x)
const {
auto const element = hsearch_.findEntity(x);
return std::make_pair(element, element->geometry().local(x));
}
std::fstream writer_;
Dune::HierarchicSearch<typename GridView::Grid, GridView> const hsearch_;
std::pair<ElementPointer, LocalVector> const G;
std::pair<ElementPointer, LocalVector> const H;
std::pair<ElementPointer, LocalVector> const J;
std::pair<ElementPointer, LocalVector> const I;
std::pair<ElementPointer, LocalVector> const U;
std::pair<ElementPointer, LocalVector> const Z;
public:
SpecialWriter(std::string filename, GridView const &gridView)
: writer_(filename, std::fstream::out),
hsearch_(gridView.grid(), gridView),
G(globalToLocal(MyGeometry::G)),
H(globalToLocal(MyGeometry::H)),
J(globalToLocal(MyGeometry::J)),
I(globalToLocal(MyGeometry::I)),
U(globalToLocal(MyGeometry::U)),
Z(globalToLocal(MyGeometry::Z)) {
writer_ << "Gh Hh Jh Ih Uv Uh Zv Zh" << std::endl;
}
void write(VirtualGridFunction<typename GridView::Grid, LocalVector> const &
specialField) {
LocalVector value;
specialField.evaluateLocal(*G.first, G.second, value);
writeHorizontal(value);
specialField.evaluateLocal(*H.first, H.second, value);
writeHorizontal(value);
specialField.evaluateLocal(*J.first, J.second, value);
writeHorizontal(value);
specialField.evaluateLocal(*I.first, I.second, value);
writeHorizontal(value);
specialField.evaluateLocal(*U.first, U.second, value);
writeVertical(value);
writeHorizontal(value);
specialField.evaluateLocal(*Z.first, Z.second, value);
writeVertical(value);
writeHorizontal(value);
writer_ << std::endl;
}
};
#endif
src/fixedpointiterator.cc → src/spatial-solving/fixedpointiterator.cc
View file @ 8fe950bc
...@@ -8,18 +8,22 @@ ...@@ -8,18 +8,22 @@
#include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh> #include <dune/solvers/solvers/loopsolver.hh>
#include "enums.hh" #include "../enums.hh"
#include "enumparser.hh" #include "../enumparser.hh"
#include "fixedpointiterator.hh" #include "fixedpointiterator.hh"
template <class Factory, class StateUpdater, class VelocityUpdater, void FixedPointIterationCounter::operator+=(
class ErrorNorm> FixedPointIterationCounter const &other) {
FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>:: iterations += other.iterations;
FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset, multigridIterations += other.multigridIterations;
std::shared_ptr<Nonlinearity> globalFriction, }
ErrorNorm const &errorNorm)
: factory_(factory), template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
: step_(factory.getStep()),
parset_(parset), parset_(parset),
globalFriction_(globalFriction), globalFriction_(globalFriction),
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")), fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
...@@ -30,23 +34,20 @@ FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>:: ...@@ -30,23 +34,20 @@ FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::
verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")), verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
errorNorm_(errorNorm) {} errorNorm_(errorNorm) {}
template <class Factory, class StateUpdater, class VelocityUpdater, template <class Factory, class Updaters, class ErrorNorm>
class ErrorNorm> FixedPointIterationCounter
int FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run( FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
std::shared_ptr<StateUpdater> stateUpdater, Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
std::shared_ptr<VelocityUpdater> velocityUpdater,
Matrix const &velocityMatrix, Vector const &velocityRHS,
Vector &velocityIterate) { Vector &velocityIterate) {
auto multigridStep = factory_.getSolver();
EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix); EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
LoopSolver<Vector> velocityProblemSolver( LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
multigridStep, velocityMaxIterations_, velocityTolerance_, &energyNorm, velocityTolerance_, &energyNorm,
verbosity_, false); // absolute error verbosity_, false); // absolute error
size_t fixedPointIteration; size_t fixedPointIteration;
size_t multigridIterations = 0;
ScalarVector alpha; ScalarVector alpha;
stateUpdater->extractAlpha(alpha); updaters.state_->extractAlpha(alpha);
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_; for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
++fixedPointIteration) { ++fixedPointIteration) {
// solve a velocity problem // solve a velocity problem
...@@ -54,21 +55,24 @@ int FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run( ...@@ -54,21 +55,24 @@ int FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run(
ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_, ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
velocityRHS, velocityIterate); velocityRHS, velocityIterate);
BlockProblem velocityProblem(parset_, convexProblem); BlockProblem velocityProblem(parset_, convexProblem);
multigridStep->setProblem(velocityIterate, velocityProblem); step_->setProblem(velocityIterate, velocityProblem);
velocityProblemSolver.preprocess(); velocityProblemSolver.preprocess();
velocityProblemSolver.solve(); velocityProblemSolver.solve();
multigridIterations += velocityProblemSolver.getResult().iterations;
Vector v_m; Vector v_m;
velocityUpdater->extractOldVelocity(v_m); updaters.rate_->extractOldVelocity(v_m);
v_m *= 1.0 - lambda_; v_m *= 1.0 - lambda_;
Arithmetic::addProduct(v_m, lambda_, velocityIterate); Arithmetic::addProduct(v_m, lambda_, velocityIterate);
// solve a state problem // solve a state problem
stateUpdater->solve(v_m); updaters.state_->solve(v_m);
ScalarVector newAlpha; ScalarVector newAlpha;
stateUpdater->extractAlpha(newAlpha); updaters.state_->extractAlpha(newAlpha);
if (errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) { if (lambda_ < 1e-12 or
errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
fixedPointIteration++; fixedPointIteration++;
break; break;
} }
...@@ -77,9 +81,20 @@ int FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run( ...@@ -77,9 +81,20 @@ int FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run(
if (fixedPointIteration == fixedPointMaxIterations_) if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge"); DUNE_THROW(Dune::Exception, "FPI failed to converge");
velocityUpdater->postProcess(velocityIterate); updaters.rate_->postProcess(velocityIterate);
// Cannot use return { fixedPointIteration, multigridIterations };
// with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
FixedPointIterationCounter ret;
ret.iterations = fixedPointIteration;
ret.multigridIterations = multigridIterations;
return ret;
}
return fixedPointIteration; std::ostream &operator<<(std::ostream &stream,
FixedPointIterationCounter const &fpic) {
return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
<< ")";
} }
#include "fixedpointiterator_tmpl.cc" #include "fixedpointiterator_tmpl.cc"
src/fixedpointiterator.hh → src/spatial-solving/fixedpointiterator.hh
View file @ 8fe950bc
#ifndef SRC_FIXEDPOINTITERATOR_HH #ifndef SRC_SPATIAL_SOLVING_FIXEDPOINTITERATOR_HH
#define SRC_FIXEDPOINTITERATOR_HH #define SRC_SPATIAL_SOLVING_FIXEDPOINTITERATOR_HH
#include <memory> #include <memory>
...@@ -8,10 +8,19 @@ ...@@ -8,10 +8,19 @@
#include <dune/solvers/norms/norm.hh> #include <dune/solvers/norms/norm.hh>
#include <dune/solvers/solvers/solver.hh> #include <dune/solvers/solvers/solver.hh>
template <class Factory, class StateUpdater, class VelocityUpdater, struct FixedPointIterationCounter {
class ErrorNorm> size_t iterations = 0;
size_t multigridIterations = 0;
void operator+=(FixedPointIterationCounter const &other);
};
std::ostream &operator<<(std::ostream &stream,
FixedPointIterationCounter const &fpic);
template <class Factory, class Updaters, class ErrorNorm>
class FixedPointIterator { class FixedPointIterator {
using ScalarVector = typename StateUpdater::ScalarVector; using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
using Vector = typename Factory::Vector; using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix; using Matrix = typename Factory::Matrix;
using ConvexProblem = typename Factory::ConvexProblem; using ConvexProblem = typename Factory::ConvexProblem;
...@@ -23,13 +32,13 @@ class FixedPointIterator { ...@@ -23,13 +32,13 @@ class FixedPointIterator {
std::shared_ptr<Nonlinearity> globalFriction, std::shared_ptr<Nonlinearity> globalFriction,
ErrorNorm const &errorNorm_); ErrorNorm const &errorNorm_);
int run(std::shared_ptr<StateUpdater> stateUpdater, FixedPointIterationCounter run(Updaters updaters,
std::shared_ptr<VelocityUpdater> velocityUpdater, Matrix const &velocityMatrix,
Matrix const &velocityMatrix, Vector const &velocityRHS, Vector const &velocityRHS,
Vector &velocityIterate); Vector &velocityIterate);
private: private:
Factory &factory_; std::shared_ptr<typename Factory::Step> step_;
Dune::ParameterTree const &parset_; Dune::ParameterTree const &parset_;
std::shared_ptr<Nonlinearity> globalFriction_; std::shared_ptr<Nonlinearity> globalFriction_;
......
src/spatial-solving/fixedpointiterator_tmpl.cc 0 → 100644
View file @ 8fe950bc
#ifndef MY_DIM
#error MY_DIM unset
#endif
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
#include <dune/common/function.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/tnnmg/problem-classes/convexproblem.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/myblockproblem.hh>
#include "../time-stepping/rate/rateupdater.hh"
#include "../time-stepping/state/stateupdater.hh"
#include "../time-stepping/updaters.hh"
#include "solverfactory.hh"
using Function = Dune::VirtualFunction<double, double>;
using Factory = SolverFactory<
MY_DIM,
MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
Grid>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
template class FixedPointIterator<Factory, MyUpdaters, ErrorNorm>;
src/solverfactory.cc → src/spatial-solving/solverfactory.cc
View file @ 8fe950bc
...@@ -13,15 +13,16 @@ ...@@ -13,15 +13,16 @@
template <size_t dim, class BlockProblem, class Grid> template <size_t dim, class BlockProblem, class Grid>
SolverFactory<dim, BlockProblem, Grid>::SolverFactory( SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
Dune::ParameterTree const &parset, size_t refinements, Grid const &grid, Dune::ParameterTree const &parset, Grid const &grid,
Dune::BitSetVector<dim> const &ignoreNodes) Dune::BitSetVector<dim> const &ignoreNodes)
: baseEnergyNorm(linearBaseSolverStep), : baseEnergyNorm(linearBaseSolverStep),
linearBaseSolver(&linearBaseSolverStep, linearBaseSolver(&linearBaseSolverStep,
parset.get<size_t>("linear.maxiumumIterations"), parset.get<size_t>("linear.maxiumumIterations"),
parset.get<double>("linear.tolerance"), &baseEnergyNorm, parset.get<double>("linear.tolerance"), &baseEnergyNorm,
Solver::QUIET), Solver::QUIET),
transferOperators(refinements), transferOperators(grid.maxLevel()),
multigridStep(new Solver(linearIterationStep, nonlinearSmoother)) { multigridStep(
std::make_shared<Step>(linearIterationStep, nonlinearSmoother)) {
// linear iteration step // linear iteration step
linearIterationStep.setMGType(parset.get<int>("linear.cycle"), linearIterationStep.setMGType(parset.get<int>("linear.cycle"),
parset.get<int>("linear.pre"), parset.get<int>("linear.pre"),
...@@ -47,12 +48,11 @@ template <size_t dim, class BlockProblem, class Grid> ...@@ -47,12 +48,11 @@ template <size_t dim, class BlockProblem, class Grid>
SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() { SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
for (auto &&x : transferOperators) for (auto &&x : transferOperators)
delete x; delete x;
delete multigridStep;
} }
template <size_t dim, class BlockProblem, class Grid> template <size_t dim, class BlockProblem, class Grid>
auto SolverFactory<dim, BlockProblem, Grid>::getSolver() -> Solver *{ auto SolverFactory<dim, BlockProblem, Grid>::getStep()
-> std::shared_ptr<Step> {
return multigridStep; return multigridStep;
} }
......
src/solverfactory.hh → src/spatial-solving/solverfactory.hh
View file @ 8fe950bc
#ifndef SRC_SOLVERFACTORY_HH #ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
#define SRC_SOLVERFACTORY_HH #define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
#include <dune/common/bitsetvector.hh> #include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh> #include <dune/common/parametertree.hh>
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wsign-compare"
#include <dune/solvers/iterationsteps/multigridstep.hh> #include <dune/solvers/iterationsteps/multigridstep.hh>
#pragma clang diagnostic pop
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh> #include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh> #include <dune/solvers/solvers/loopsolver.hh>
...@@ -26,16 +22,17 @@ class SolverFactory { ...@@ -26,16 +22,17 @@ class SolverFactory {
private: private:
using NonlinearSmoother = GenericNonlinearGS<BlockProblem>; using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
using Solver =
TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
public: public:
SolverFactory(Dune::ParameterTree const &parset, size_t refinements, using Step =
Grid const &grid, Dune::BitSetVector<dim> const &ignoreNodes); TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
Dune::BitSetVector<dim> const &ignoreNodes);
~SolverFactory(); ~SolverFactory();
Solver *getSolver(); std::shared_ptr<Step> getStep();
private: private:
TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep; TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
...@@ -46,6 +43,6 @@ class SolverFactory { ...@@ -46,6 +43,6 @@ class SolverFactory {
MultigridStep<Matrix, Vector> linearIterationStep; MultigridStep<Matrix, Vector> linearIterationStep;
std::vector<CompressedMultigridTransfer<Vector> *> transferOperators; std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
NonlinearSmoother nonlinearSmoother; NonlinearSmoother nonlinearSmoother;
Solver *multigridStep; std::shared_ptr<Step> multigridStep;
}; };
#endif #endif
src/solverfactory_tmpl.cc → src/spatial-solving/solverfactory_tmpl.cc
View file @ 8fe950bc
...@@ -2,8 +2,8 @@ ...@@ -2,8 +2,8 @@
#error MY_DIM unset #error MY_DIM unset
#endif #endif
#include "explicitgrid.hh" #include "../explicitgrid.hh"
#include "explicitvectors.hh" #include "../explicitvectors.hh"
#include <dune/tnnmg/nonlinearities/zerononlinearity.hh> #include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh> #include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
...@@ -17,6 +17,6 @@ template class SolverFactory< ...@@ -17,6 +17,6 @@ template class SolverFactory<
MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>, MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
Grid>; Grid>;
template class SolverFactory< template class SolverFactory<
MY_DIM, BlockNonlinearTNNMGProblem< MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
ConvexProblem<ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>, ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
Grid>; Grid>;
src/state_tmpl.cc deleted 100644 → 0
View file @ 4a97dcf2
#include "explicitvectors.hh"
template std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater<
ScalarVector, Vector>(Config::stateModel model,
ScalarVector const &alpha_initial,
Dune::BitSetVector<1> const &frictionalNodes,
double L, double V0);
src/time-stepping/adaptivetimestepper.cc 0 → 100644
View file @ 8fe950bc
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "adaptivetimestepper.hh"
void IterationRegister::registerCount(FixedPointIterationCounter count) {
totalCount += count;
}
void IterationRegister::registerFinalCount(FixedPointIterationCounter count) {
finalCount = count;
}
void IterationRegister::reset() {
totalCount = FixedPointIterationCounter();
finalCount = FixedPointIterationCounter();
}
template <class Factory, class Updaters, class ErrorNorm>
AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, Updaters &current,
double relativeTime, double relativeTau,
std::function<void(double, Vector &)> externalForces,
ErrorNorm const &errorNorm,
std::function<bool(Updaters &, Updaters &)> mustRefine)
: relativeTime_(relativeTime),
relativeTau_(relativeTau),
finalTime_(parset.get<double>("problem.finalTime")),
factory_(factory),
parset_(parset),
globalFriction_(globalFriction),
current_(current),
R1_(),
externalForces_(externalForces),
mustRefine_(mustRefine),
errorNorm_(errorNorm) {}
template <class Factory, class Updaters, class ErrorNorm>
bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
return relativeTime_ >= 1.0;
}
template <class Factory, class Updaters, class ErrorNorm>
IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
/*
| C | We check here if making the step R1 of size tau is a
| R1 | R2 | good idea. To check if we can coarsen, we compare
|F1|F2| | the result of (R1+R2) with C, i.e. two steps of size
tau with one of size 2*tau. To check if we need to
refine, we compare the result of (F1+F2) with R1, i.e. two steps
of size tau/2 with one of size tau. The method makes multiple
coarsening/refining attempts, with coarsening coming first. */
if (R1_.updaters == Updaters())
R1_ = step(current_, relativeTime_, relativeTau_);
bool didCoarsen = false;
iterationRegister_.reset();
UpdatersWithCount R2;
UpdatersWithCount C;
while (relativeTime_ + relativeTau_ <= 1.0) {
R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
C = step(current_, relativeTime_, 2 * relativeTau_);
if (mustRefine_(R2.updaters, C.updaters))
break;
didCoarsen = true;
relativeTau_ *= 2;
R1_ = C;
}
UpdatersWithCount F1;
UpdatersWithCount F2;
if (!didCoarsen) {
while (true) {
F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
relativeTau_ / 2.0);
if (!mustRefine_(F2.updaters, R1_.updaters))
break;
relativeTau_ /= 2.0;
R1_ = F1;
R2 = F2;
}
}
iterationRegister_.registerFinalCount(R1_.count);
relativeTime_ += relativeTau_;
current_ = R1_.updaters;
R1_ = R2;
return iterationRegister_;
}
template <class Factory, class Updaters, class ErrorNorm>
typename AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::UpdatersWithCount
AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
Updaters const &oldUpdaters, double rTime, double rTau) {
UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
newUpdatersAndCount.count =
MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
newUpdatersAndCount.updaters, errorNorm_,
externalForces_)
.step(rTime, rTau);
iterationRegister_.registerCount(newUpdatersAndCount.count);
return newUpdatersAndCount;
}
#include "adaptivetimestepper_tmpl.cc"
src/time-stepping/adaptivetimestepper.hh 0 → 100644
View file @ 8fe950bc
#ifndef SRC_TIME_STEPPING_ADAPTIVETIMESTEPPER_HH
#define SRC_TIME_STEPPING_ADAPTIVETIMESTEPPER_HH
#include <fstream>
#include "coupledtimestepper.hh"
struct IterationRegister {
void registerCount(FixedPointIterationCounter count);
void registerFinalCount(FixedPointIterationCounter count);
void reset();
FixedPointIterationCounter totalCount;
FixedPointIterationCounter finalCount;
};
template <class Factory, class Updaters, class ErrorNorm>
class AdaptiveTimeStepper {
struct UpdatersWithCount {
Updaters updaters;
FixedPointIterationCounter count;
};
using Vector = typename Factory::Vector;
using ConvexProblem = typename Factory::ConvexProblem;
using Nonlinearity = typename ConvexProblem::NonlinearityType;
using MyCoupledTimeStepper = CoupledTimeStepper<Factory, Updaters, ErrorNorm>;
public:
AdaptiveTimeStepper(Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction,
Updaters &current, double relativeTime,
double relativeTau,
std::function<void(double, Vector &)> externalForces,
ErrorNorm const &errorNorm,
std::function<bool(Updaters &, Updaters &)> mustRefine);
bool reachedEnd();
IterationRegister advance();
double relativeTime_;
double relativeTau_;
private:
UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
double rTau);
double finalTime_;
Factory &factory_;
Dune::ParameterTree const &parset_;
std::shared_ptr<Nonlinearity> globalFriction_;
Updaters &current_;
UpdatersWithCount R1_;
std::function<void(double, Vector &)> externalForces_;
std::function<bool(Updaters &, Updaters &)> mustRefine_;
ErrorNorm const &errorNorm_;
IterationRegister iterationRegister_;
};
#endif
src/fixedpointiterator_tmpl.cc → src/time-stepping/adaptivetimestepper_tmpl.cc
View file @ 8fe950bc
...@@ -2,6 +2,9 @@ ...@@ -2,6 +2,9 @@
#error MY_DIM unset #error MY_DIM unset
#endif #endif
#include "../explicitgrid.hh"
#include "../explicitvectors.hh"
#include <dune/common/function.hh> #include <dune/common/function.hh>
#include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/norms/energynorm.hh>
...@@ -10,12 +13,10 @@ ...@@ -10,12 +13,10 @@
#include <dune/tectonic/globalfriction.hh> #include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/myblockproblem.hh> #include <dune/tectonic/myblockproblem.hh>
#include "explicitgrid.hh" #include "../spatial-solving/solverfactory.hh"
#include "explicitvectors.hh" #include "rate/rateupdater.hh"
#include "solverfactory.hh"
#include "state/stateupdater.hh" #include "state/stateupdater.hh"
#include "timestepping.hh" #include "updaters.hh"
using Function = Dune::VirtualFunction<double, double>; using Function = Dune::VirtualFunction<double, double>;
using Factory = SolverFactory< using Factory = SolverFactory<
...@@ -23,9 +24,9 @@ using Factory = SolverFactory< ...@@ -23,9 +24,9 @@ using Factory = SolverFactory<
MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>, MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
Grid>; Grid>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>; using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyVelocityUpdater = TimeSteppingScheme<Vector, Matrix, Function, MY_DIM>; using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>; using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
template class FixedPointIterator<Factory, MyStateUpdater, MyVelocityUpdater, template class AdaptiveTimeStepper<Factory, MyUpdaters, ErrorNorm>;
ErrorNorm>;
src/time-stepping/coupledtimestepper.cc 0 → 100644
View file @ 8fe950bc
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include "coupledtimestepper.hh"
template <class Factory, class Updaters, class ErrorNorm>
CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
double finalTime, Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, Updaters updaters,
ErrorNorm const &errorNorm,
std::function<void(double, Vector &)> externalForces)
: finalTime_(finalTime),
factory_(factory),
parset_(parset),
globalFriction_(globalFriction),
updaters_(updaters),
externalForces_(externalForces),
errorNorm_(errorNorm) {}
template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
double relativeTau) {
updaters_.state_->nextTimeStep();
updaters_.rate_->nextTimeStep();
auto const newRelativeTime = relativeTime + relativeTau;
Vector ell;
externalForces_(newRelativeTime, ell);
Matrix velocityMatrix;
Vector velocityRHS;
Vector velocityIterate;
auto const tau = relativeTau * finalTime_;
updaters_.state_->setup(tau);
updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS,
velocityIterate, velocityMatrix);
FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
factory_, parset_, globalFriction_, errorNorm_);
auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
velocityRHS, velocityIterate);
return iterations;
}
#include "coupledtimestepper_tmpl.cc"
src/coupledtimestepper.hh → src/time-stepping/coupledtimestepper.hh
View file @ 8fe950bc
#ifndef SRC_COUPLEDTIMESTEPPER_HH #ifndef SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
#define SRC_COUPLEDTIMESTEPPER_HH #define SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
#include <functional> #include <functional>
#include <memory> #include <memory>
#include <dune/common/parametertree.hh> #include <dune/common/parametertree.hh>
template <class Factory, class StateUpdater, class VelocityUpdater, #include "../spatial-solving/fixedpointiterator.hh"
class ErrorNorm>
template <class Factory, class Updaters, class ErrorNorm>
class CoupledTimeStepper { class CoupledTimeStepper {
using Vector = typename Factory::Vector; using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix; using Matrix = typename Factory::Matrix;
...@@ -18,20 +19,17 @@ class CoupledTimeStepper { ...@@ -18,20 +19,17 @@ class CoupledTimeStepper {
CoupledTimeStepper(double finalTime, Factory &factory, CoupledTimeStepper(double finalTime, Factory &factory,
Dune::ParameterTree const &parset, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, std::shared_ptr<Nonlinearity> globalFriction,
std::shared_ptr<StateUpdater> stateUpdater, Updaters updaters, ErrorNorm const &errorNorm,
std::shared_ptr<VelocityUpdater> velocityUpdater,
ErrorNorm const &errorNorm,
std::function<void(double, Vector &)> externalForces); std::function<void(double, Vector &)> externalForces);
int step(double relativeTime, double relativeTau); FixedPointIterationCounter step(double relativeTime, double relativeTau);
private: private:
double finalTime_; double finalTime_;
Factory &factory_; Factory &factory_;
Dune::ParameterTree const &parset_; Dune::ParameterTree const &parset_;
std::shared_ptr<Nonlinearity> globalFriction_; std::shared_ptr<Nonlinearity> globalFriction_;
std::shared_ptr<StateUpdater> stateUpdater_; Updaters updaters_;
std::shared_ptr<VelocityUpdater> velocityUpdater_;
std::function<void(double, Vector &)> externalForces_; std::function<void(double, Vector &)> externalForces_;
ErrorNorm const &errorNorm_; ErrorNorm const &errorNorm_;
}; };
......