Skip to content
Snippets Groups Projects
Commit 2cfdb544 authored by podlesny's avatar podlesny
Browse files

switch to new io

parent c5a80415
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 33 additions and 117 deletions
src/foam/foam.cc
+ 33
− 117
View file @ 2cfdb544
......@@ -56,6 +56,7 @@
#include <dune/tectonic/factories/twoblocksfactory.hh>
#include <dune/tectonic/io/io-handler.hh>
#include <dune/tectonic/io/hdf5-writer.hh>
#include <dune/tectonic/io/hdf5/restart-io.hh>
#include <dune/tectonic/io/vtk.hh>
......@@ -88,12 +89,13 @@
std::vector<std::vector<double>> allReductionFactors;
size_t const dims = MY_DIM;
const std::string sourcePath = "/home/joscha/software/dune/dune-tectonic/src/foam/";
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/foam/foam.cfg", parset);
Dune::ParameterTreeParser::readINITree(sourcePath + "foam.cfg", parset);
Dune::ParameterTreeParser::readINITree(
Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/foam/foam-%dD.cfg", dims), parset);
Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
return parset;
}
......@@ -112,12 +114,22 @@ int main(int argc, char *argv[]) {
std::cout << argv[0] << std::endl;
}
auto const parset = getParameters(argc, argv);
auto outPath = std::filesystem::current_path();
outPath += "/output/" + parset.get<std::string>("outPath");
if (!std::filesystem::is_directory(outPath))
std::filesystem::create_directories(outPath);
const auto copyOptions = std::filesystem::copy_options::overwrite_existing;
std::filesystem::copy(sourcePath + "foam.cfg", outPath, copyOptions);
std::filesystem::copy(Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), outPath, copyOptions);
std::filesystem::current_path(outPath);
std::ofstream out("foam.log");
std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
auto const parset = getParameters(argc, argv);
using Assembler = MyAssembler<DefLeafGridView, dims>;
using field_type = Matrix::field_type;
......@@ -131,8 +143,6 @@ int main(int argc, char *argv[]) {
ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
const size_t bodyCount = contactNetwork.nBodies();
for (size_t i=0; i<contactNetwork.nLevels(); i++) {
......@@ -145,12 +155,12 @@ int main(int argc, char *argv[]) {
const auto& level = *contactNetwork.level(i);
for (size_t j=0; j<level.nBodies(); j++) {
writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
//writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
}
}
for (size_t i=0; i<bodyCount; i++) {
writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
//writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
}
// ----------------------------
......@@ -170,33 +180,13 @@ int main(int argc, char *argv[]) {
using MyProgramState = ProgramState<Vector, ScalarVector>;
MyProgramState programState(nVertices);
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;
IOHandler<Assembler, ContactNetwork> ioHandler(parset.sub("io"), contactNetwork);
auto restartIO = handleRestarts
? std::make_unique<RestartIO<MyProgramState>>(
*restartFile, nVertices)
: nullptr;
if (firstRestart > 0) // automatically adjusts the time and timestep
restartIO->read(firstRestart, programState);
else
programState.setupInitialConditions(parset, contactNetwork);
bool restartRead = ioHandler.read(programState);
if (!restartRead) {
programState.setupInitialConditions(parset, contactNetwork);
}
auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
......@@ -213,87 +203,9 @@ int main(int argc, char *argv[]) {
auto& globalFriction = contactNetwork.globalFriction();
globalFriction.updateAlpha(programState.alpha);
using MyVertexBasis = typename Assembler::VertexBasis;
using MyCellBasis = typename Assembler::CellBasis;
std::vector<Vector> vertexCoordinates(bodyCount);
std::vector<const MyVertexBasis* > vertexBases(bodyCount);
std::vector<const MyCellBasis* > cellBases(bodyCount);
auto& wPatches = blocksFactory.weakPatches();
std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
for (size_t i=0; i<bodyCount; i++) {
const auto& body = contactNetwork.body(i);
vertexBases[i] = &(body->assembler()->vertexBasis);
cellBases[i] = &(body->assembler()->cellBasis);
weakPatches[i].resize(1);
weakPatches[i][0] = wPatches[i].get();
auto& vertexCoords = vertexCoordinates[i];
vertexCoords.resize(nVertices[i]);
auto hostLeafView = body->grid()->hostGrid().leafGridView();
Dune::MultipleCodimMultipleGeomTypeMapper<
LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
for (auto &&v : vertices(hostLeafView))
vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
}
typename ContactNetwork::BoundaryPatches frictionBoundaries;
contactNetwork.boundaryPatches("friction", frictionBoundaries);
std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
contactNetwork.frictionPatches(frictionPatches);
auto dataWriter =
writeData ? std::make_unique<
HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
*dataFile, vertexCoordinates, vertexBases,
frictionPatches, weakPatches)
: nullptr;
const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/home/joscha/software/dune/dune-tectonic/body");
IterationRegister iterationCount;
auto const report = [&](bool initial = false) {
if (writeData) {
dataWriter->reportSolution(programState, contactNetwork, globalFriction);
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")) {
std::vector<ScalarVector> stress(bodyCount);
for (size_t i=0; i<bodyCount; i++) {
const auto& body = contactNetwork.body(i);
body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
body->data()->getPoissonRatio(),
programState.u[i], stress[i]);
}
vtkWriter.write(programState.timeStep, programState.u, programState.v,
programState.alpha, stress);
}
};
report(true);
ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, true);
// -------------------
// Set up TNNMG solver
......@@ -398,8 +310,8 @@ int main(int argc, char *argv[]) {
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;
//std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance << std::endl;
//std::cout << "must refine: " << (energyNorm > refinementTolerance) << std::endl;
return energyNorm > refinementTolerance;
};
......@@ -432,10 +344,14 @@ int main(int argc, char *argv[]) {
typename ContactNetwork::ExternalForces externalForces;
contactNetwork.externalForces(externalForces);
StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
externalForces, stateEnergyNorms);
AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
adaptiveTimeStepper(stepBase, contactNetwork, current,
programState.relativeTime, programState.relativeTau,
externalForces, stateEnergyNorms, mustRefine);
mustRefine);
size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
......@@ -460,7 +376,7 @@ int main(int argc, char *argv[]) {
contactNetwork.setDeformation(programState.u);
report();
ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, false);
if (programState.timeStep==timeSteps) {
std::cout << "limit of timeSteps reached!" << std::endl;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment