Forked from
agnumpde / dune-tectonic
Up to date with the upstream repository.
-
Elias Pipping authored
- The configuration parameter io.writeVTK is now io.vtk.write - The configuration parameters restarts.first and restarts.spacing were moved into the io namespace - Restarts are stored in a separate file - VTK data, restarts (in HDF5), and the other HDF5 data are now each controlled through a parameter: io.data.write, io.restarts.write, io.vtk.write In particular, it is thus now possible to write VTK data while opening the data and restart files in read-only mode (or not at all); This makes it possible to, by keeping the restarts file (and the binary!) around, quickly (re-)compute a a small segment from the history. The likeliest scenario: VTK files are huge and for a computation with tens of thousands of timesteps (e.g. 30000), we know already (from the data, or just expect it) that an event occurred during the last 500 timesteps. So we recompute them from a restart in the following mode: -io.restarts.first 29500 -io.restarts.write false -io.data.write false -io.vtk.write true Note that care needs to be taken to use a number for io.restarts.first which is compatible with the old io.restarts.spacingElias Pipping authored- The configuration parameter io.writeVTK is now io.vtk.write - The configuration parameters restarts.first and restarts.spacing were moved into the io namespace - Restarts are stored in a separate file - VTK data, restarts (in HDF5), and the other HDF5 data are now each controlled through a parameter: io.data.write, io.restarts.write, io.vtk.write In particular, it is thus now possible to write VTK data while opening the data and restart files in read-only mode (or not at all); This makes it possible to, by keeping the restarts file (and the binary!) around, quickly (re-)compute a a small segment from the history. The likeliest scenario: VTK files are huge and for a computation with tens of thousands of timesteps (e.g. 30000), we know already (from the data, or just expect it) that an event occurred during the last 500 timesteps. So we recompute them from a restart in the following mode: -io.restarts.first 29500 -io.restarts.write false -io.data.write false -io.vtk.write true Note that care needs to be taken to use a number for io.restarts.first which is compatible with the old io.restarts.spacing
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
one-body-problem.cc 12.48 KiB
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_IPOPT
#undef HAVE_IPOPT
#endif
#include <atomic>
#include <cmath>
#include <csignal>
#include <exception>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <dune/common/bitsetvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/formatstring.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/solvers/solver.hh>
#include <dune/tnnmg/problem-classes/convexproblem.hh>
#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/fufem/hdf5/file.hh>
#include "assemblers.hh"
#include "diameter.hh"
#include "enumparser.hh"
#include "enums.hh"
#include "gridselector.hh"
#include "hdf5-writer.hh"
#include "hdf5/restart-io.hh"
#include "matrices.hh"
#include "program_state.hh"
#include "one-body-problem-data/bc.hh"
#include "one-body-problem-data/mybody.hh"
#include "one-body-problem-data/mygeometry.hh"
#include "one-body-problem-data/myglobalfrictiondata.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"
size_t const dims = MY_DIM;
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
Dune::ParameterTreeParser::readINITree( 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[]) {
try {
Dune::MPIHelper::instance(argc, argv);
auto const parset = getParameters(argc, argv);
MyGeometry::render();
MyGeometry::write();
using GridView = Grid::LeafGridView;
using MyAssembler = MyAssembler<GridView, dims>;
using Matrix = MyAssembler::Matrix;
using Vector = MyAssembler::Vector;
using LocalVector = Vector::block_type;
using ScalarMatrix = MyAssembler::ScalarMatrix;
using ScalarVector = MyAssembler::ScalarVector;
auto const weakPatch =
getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"));
// {{{ Set up grid
GridConstructor<Grid> gridConstructor;
auto grid = gridConstructor.getGrid();
refine(*grid, weakPatch,
parset.get<double>("boundary.friction.smallestDiameter"));
double minDiameter = std::numeric_limits<double>::infinity();
double maxDiameter = 0.0;
for (auto &&e : elements(grid->leafGridView())) {
auto const geometry = e.geometry();
auto const diam = diameter(geometry);
minDiameter = std::min(minDiameter, diam);
maxDiameter = std::max(maxDiameter, diam);
}
std::cout << "min diameter: " << minDiameter << std::endl;
std::cout << "max diameter: " << maxDiameter << std::endl;
auto const leafView = grid->leafGridView();
auto const leafVertexCount = leafView.size(dims);
std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
auto myFaces = gridConstructor.constructFaces(leafView);
BoundaryPatch<GridView> const neumannBoundary(leafView);
BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
BoundaryPatch<GridView> const &surface = myFaces.upper;
// Dirichlet Boundary
Dune::BitSetVector<dims> noNodes(leafVertexCount);
Dune::BitSetVector<dims> dirichletNodes(leafVertexCount);
for (size_t i = 0; i < leafVertexCount; ++i) {
if (myFaces.right.containsVertex(i))
dirichletNodes[i][0] = true;
if (myFaces.lower.containsVertex(i))
dirichletNodes[i][1] = true;
#if MY_DIM == 3
if (myFaces.front.containsVertex(i) || myFaces.back.containsVertex(i))
dirichletNodes[i][2] = true;
#endif }
// Set up functions for time-dependent boundary conditions
using Function = Dune::VirtualFunction<double, double>;
Function const &velocityDirichletFunction = VelocityDirichletCondition();
Function const &neumannFunction = NeumannCondition();
MyAssembler const myAssembler(leafView);
MyBody<dims> const body(parset);
Matrices<Matrix> matrices;
myAssembler.assembleElasticity(body.getYoungModulus(),
body.getPoissonRatio(), matrices.elasticity);
myAssembler.assembleViscosity(body.getShearViscosityField(),
body.getBulkViscosityField(),
matrices.damping);
myAssembler.assembleMass(body.getDensityField(), matrices.mass);
ScalarMatrix relativeFrictionalBoundaryMass;
myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
relativeFrictionalBoundaryMass);
relativeFrictionalBoundaryMass /= frictionalBoundary.area();
EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
relativeFrictionalBoundaryMass);
// Assemble forces
Vector gravityFunctional;
myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
// Problem formulation: right-hand side
std::function<void(double, Vector &)> computeExternalForces =
[&](double _relativeTime, Vector &_ell) {
myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
_relativeTime);
_ell += gravityFunctional;
};
using MyProgramState = ProgramState<Vector, ScalarVector>;
MyProgramState programState(leafVertexCount);
auto const firstRestart = parset.get<size_t>("io.restarts.first");
auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
auto const writeRestarts = parset.get<bool>("io.restarts.write");
auto const writeData = parset.get<bool>("io.data.write");
bool const handleRestarts = writeRestarts or firstRestart > 0;
auto dataFile =
writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
auto restartFile = handleRestarts
? std::make_unique<HDF5::File>(
"restarts.h5",
writeRestarts ? HDF5::Access::READWRITE
: HDF5::Access::READONLY)
: nullptr;
auto restartIO = handleRestarts
? std::make_unique<RestartIO<MyProgramState>>(
*restartFile, leafVertexCount)
: nullptr;
if (firstRestart > 0) // automatically adjusts the time and timestep
restartIO->read(firstRestart, programState);
else
programState.setupInitialConditions(parset, computeExternalForces,
matrices, myAssembler, dirichletNodes,
noNodes, frictionalBoundary, body);
MyGlobalFrictionData<LocalVector> frictionInfo(
parset.sub("boundary.friction"), weakPatch);
auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), frictionalBoundary, frictionInfo, programState.weightedNormalStress);
myGlobalFriction->updateAlpha(programState.alpha);
Vector vertexCoordinates(leafVertexCount);
{
Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
for (auto &&v : vertices(leafView))
vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
using MyVertexBasis = typename MyAssembler::VertexBasis;
auto dataWriter =
writeData ? std::make_unique<
HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
*dataFile, vertexCoordinates, myAssembler.vertexBasis,
surface, frictionalBoundary, weakPatch)
: nullptr;
MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
IterationRegister iterationCount;
auto const report = [&](bool initial = false) {
if (writeData) {
dataWriter->reportSolution(programState, *myGlobalFriction);
if (!initial)
dataWriter->reportIterations(programState, iterationCount);
dataFile->flush();
}
if (writeRestarts and !initial and
programState.timeStep % restartSpacing == 0) {
restartIO->write(programState);
restartFile->flush();
}
if (parset.get<bool>("io.printProgress"))
std::cout << "timeStep = " << std::setw(6) << programState.timeStep
<< ", time = " << std::setw(12) << programState.relativeTime
<< ", tau = " << std::setw(12) << programState.relativeTau
<< std::endl;
if (parset.get<bool>("io.vtk.write")) {
ScalarVector stress;
myAssembler.assembleVonMisesStress(body.getYoungModulus(),
body.getPoissonRatio(),
programState.u, stress);
vtkWriter.write(programState.timeStep, programState.u, programState.v,
programState.alpha, stress);
}
};
report(true);
// Set up TNNMG solver
using NonlinearFactory = SolverFactory<
dims,
MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
Grid>;
NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes);
using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
StateUpdater<ScalarVector, Vector>>;
MyUpdater current(
initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, dirichletNodes, matrices,
programState.u, programState.v, programState.a),
initStateUpdater<ScalarVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
programState.alpha, *frictionalBoundary.getVertices(),
parset.get<double>("boundary.friction.L"), parset.get<double>("boundary.friction.V0")));
auto const refinementTolerance =
parset.get<double>("timeSteps.refinementTolerance");
auto const mustRefine = [&](MyUpdater &coarseUpdater,
MyUpdater &fineUpdater) {
ScalarVector coarseAlpha;
coarseUpdater.state_->extractAlpha(coarseAlpha);
ScalarVector fineAlpha;
fineUpdater.state_->extractAlpha(fineAlpha);
return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
};
std::signal(SIGXCPU, handleSignal);
std::signal(SIGINT, handleSignal);
std::signal(SIGTERM, handleSignal);
AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
EnergyNorm<ScalarMatrix, ScalarVector>>
adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
programState.relativeTime, programState.relativeTau,
computeExternalForces, stateEnergyNorm, mustRefine);
while (!adaptiveTimeStepper.reachedEnd()) {
programState.timeStep++;
iterationCount = adaptiveTimeStepper.advance();
programState.relativeTime = adaptiveTimeStepper.relativeTime_;
programState.relativeTau = adaptiveTimeStepper.relativeTau_;
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;
}
}
} catch (Dune::Exception &e) {
Dune::derr << "Dune reported error: " << e << std::endl;
} catch (std::exception &e) {
std::cerr << "Standard exception: " << e.what() << std::endl;
}
}