one-body-problem.cc

#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"

// for getcwd
#include <unistd.h>

#define USE_OLD_TNNMG

size_t const dims = MY_DIM;

Dune::ParameterTree getParameters(int argc, char *argv[]) {
  Dune::ParameterTree parset;
  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/one-body-problem.cfg", parset);
  Dune::ParameterTreeParser::readINITree(
      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/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);

    char buffer[256];
    char *val = getcwd(buffer, sizeof(buffer));
    if (val) {
        std::cout << buffer << std::endl;
        std::cout << argv[0] << std::endl;
    }

    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, Dune::mcmgVertexLayout());
      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;
  }
}