cantorfaultnetwork.cc

#include <config.h>

#include <cmath>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>

/*
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include <dune/grid/uggrid/ugwrapper.hh>
#pragma GCC diagnostic pop
*/

// dune-common includes
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/timer.hh>

// dune-istl includes
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>

// dune-fufem includes
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>

// dune-grid includes
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

// dune-solver includes
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/twonorm.hh>

// dune-faultnetworks includes
#include <dune/faultnetworks/utils/debugutils.hh>
#include <dune/faultnetworks/utils/matrixreader.hh>
#include <dune/faultnetworks/utils/vectorreader.hh>
#include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh>

#include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh>
#include <dune/faultnetworks/solvers/osccgsolver.hh>

#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
#include <dune/faultnetworks/assemblers/osclocalassembler.hh>

#include <dune/faultnetworks/oscrhs.hh>
#include <dune/faultnetworks/oscdata.hh>
#include <dune/faultnetworks/oscdatahierarchy.hh>

#include <dune/faultnetworks/faultp1nodalbasis.hh>
#include <dune/faultnetworks/interfacenetwork.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
#include <dune/faultnetworks/dgmgtransfer.hh>

#include <dune/faultnetworks/faultfactories/cantorfaultfactory.hh>

int main(int argc, char** argv) { try
{
    MPIHelper::instance(argc, argv);

    const int dim = 2;
    //typedef double ctype;


    typedef typename Dune::FieldVector<double, dim> WorldVectorType;
    typedef typename Dune::BlockVector<Dune::FieldVector<double,1>> VectorType;
    typedef typename Dune::BitSetVector<1> BitVector;
    typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> MatrixType;


    #if HAVE_UG
    typedef typename Dune::UGGrid<dim> GridType;
    #else
    #error No UG!
    #endif

    typedef typename GridType::LevelGridView GV;
    typedef FaultP1NodalBasis<GV, double> DGBasis;
    typedef typename DGBasis::LocalFiniteElement DGFE;

    typedef OscLocalAssembler<GridType, DGFE, DGFE> LocalOscAssembler;
    typedef IntersectionJumpMassAssembler<GV, DGFE, DGFE> LocalInterfaceAssembler;
 
  // parse parameter file
    ParameterTree parameterSet;
    if (argc==2)
        ParameterTreeParser::readINITree(argv[1], parameterSet);
    else
        ParameterTreeParser::readINITree("cantorfaultnetwork.parset", parameterSet);

    const std::string path = parameterSet.get<std::string>("path");
    const std::string resultPath = parameterSet.get<std::string>("resultPath");

    int problemCount = 0;
    while (parameterSet.hasSub("problem" + std::to_string(problemCount))) {
	
    ParameterTree& problemParameters = parameterSet.sub("problem" + std::to_string(problemCount));

    const std::string oscDataFile = problemParameters.get<std::string>("oscDataFile");

    const int coarseResolution = problemParameters.get<int>("coarseResolution");
    const int fineResolution = problemParameters.get<int>("fineResolution");
    const int exactResolution = problemParameters.get<int>("exactResolution");
    const int minCantorResolution = problemParameters.get<int>("minCantorResolution");

    const double c = problemParameters.get<double>("penaltyFactor");
    const int patchDepth = problemParameters.get<int>("patchDepth");

    const int maxIterations = problemParameters.get<int>("maxIterations");
    const double solverTolerance = problemParameters.get<double>("tol");
    const bool nestedIteration = problemParameters.get<bool>("nestedIteration");

    const int fineGridN = std::pow(2, fineResolution);

    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(coarseResolution) + "_" + std::to_string(fineResolution) + "_" + std::to_string(exactResolution) + ".log");
    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt


    std::cout << std::endl;
    std::cout << "Parameter file read successfully!" << std::endl;

    typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType;
    OscMatrixType matrix;

    // interface integral jump penalty, B:
    OscMatrixType B(2,2);
    B[0][0] = 1;
    B[1][1] = 1;
    B[0][1] = -1;
    B[1][0] = -1;

    MatrixReader<OscMatrixType> matrixReader(path + oscDataFile);
    matrixReader.read(matrix);

    const int oscGridN = matrix.N();
    const int oscGridLevelIdx = std::round(std::log2(oscGridN)) - coarseResolution;

    if (oscGridN>fineGridN)
        DUNE_THROW(Dune::Exception, "Provided oscData too large!");

    std::cout << "OscData file read successfully!" << std::endl << std::endl;

    const int fineLevelIdx = fineResolution - coarseResolution;
    const int exactLevelIdx = exactResolution - coarseResolution;

    size_t maxCantorLevel = exactLevelIdx+minCantorResolution;
    std::map<size_t, size_t> levelToCantorLevel;
    for (int i=0; i<=exactLevelIdx; i++) {
        levelToCantorLevel[i] = i+minCantorResolution;
    }

    // build multilevel cantor fault network
    CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, coarseResolution, exactLevelIdx, maxCantorLevel);
    const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();

    std::vector<int> writeLevels {0, 1, fineLevelIdx};
    for (size_t i=0; i<writeLevels.size(); i++) {
        int writeLevel = writeLevels[i];
        bool writeGrid = !(writeLevel==fineLevelIdx);

        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz");
        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid);
    }

    const GridType& grid = faultFactory.grid();

    OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
    oscData.set(matrix);
    OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);
    
    Timer totalTimer;
    totalTimer.reset();
    totalTimer.start();


    // ----------------------------------------------------------------
    // ---              compute initial iterate                     ---
    // ----------------------------------------------------------------
    VectorType coarseSol;
    std::shared_ptr<DGBasis> coarseBasis;

    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);

    std::cout << std::endl;
    std::cout << "Computing initial iterate!" << std::endl;

    if (!nestedIteration) {
        const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0);
        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork);
        LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper());

        std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size());
        for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) {
            const int k = minCantorResolution + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i);
            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
            coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
        }

        coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler);
        coarseGlobalAssembler.solve();

        coarseSol = coarseGlobalAssembler.getSol();
        coarseBasis = coarseGlobalAssembler.basis();
    } else {
        const LevelInterfaceNetwork<GV>& nLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx-1);
        coarseBasis = std::make_shared<DGBasis>(nLevelInterfaceNetwork);

        coarseSol.resize(coarseBasis->size());

        std::stringstream solFilename;
        solFilename << resultPath << "solutions/level_" << std::to_string(fineResolution-1);

        std::ifstream file(solFilename.str().c_str(), std::ios::in|std::ios::binary);
        if (not(file))
            DUNE_THROW(SolverError, "Couldn't open file " << solFilename.str() << " for reading");

        Dune::MatrixVector::Generic::readBinary(file, coarseSol);

        file.close();
    }

    // ----------------------------------------------------------------
    // ---              compute exact solution                      ---
    // ----------------------------------------------------------------
    std::cout << std::endl;
    std::cout << "Computing exact solution!" << std::endl;

    const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(exactLevelIdx);
    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> exactGlobalAssembler(exactLevelInterfaceNetwork);
    LocalOscAssembler exactLocalAssembler(oscDataHierarchy[exactLevelIdx]->data(), oscDataHierarchy[exactLevelIdx]->mapper());

    std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size());
    for (size_t i=0; i<exactLocalInterfaceAssemblers.size(); i++) {
        const int k = minCantorResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i);
        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
        exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
    }

    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
    exactGlobalAssembler.solve();

    const VectorType& exactSol = exactGlobalAssembler.getSol();
    std::shared_ptr<DGBasis> exactBasis = exactGlobalAssembler.basis();

    EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix());

    // ----------------------------------------------------------------
    // ---                set up cg solver                          ---
    // ----------------------------------------------------------------
    std::cout << std::endl;
    std::cout << "---------------------------------------------" << std::endl; 
    std::cout << "Setting up the fine level CG solver:" << std::endl; 
    std::cout << "---------------------------------------------" << std::endl << std::endl;

    const LevelInterfaceNetwork<GV>& fineLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx);

    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> fineGlobalAssembler(fineLevelInterfaceNetwork);
    LocalOscAssembler fineLocalAssembler(oscDataHierarchy[fineLevelIdx]->data(), oscDataHierarchy[fineLevelIdx]->mapper());

    std::vector<std::shared_ptr<LocalInterfaceAssembler>> fineLocalInterfaceAssemblers(fineLevelInterfaceNetwork.size());
    for (size_t i=0; i<fineLocalInterfaceAssemblers.size(); i++) {
        const int k = minCantorResolution + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i);
        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
        fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
    }

    fineGlobalAssembler.assemble(fineLocalAssembler, fineLocalInterfaceAssemblers, l2FunctionalAssembler);
    fineGlobalAssembler.solve();

    const VectorType& fineSol = fineGlobalAssembler.getSol();
    const VectorType& fineRhs = fineGlobalAssembler.rhs();
    std::shared_ptr<DGBasis> fineBasis = fineGlobalAssembler.basis();

    EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix());

    std::cout << "Setting up the preconditioner!" << std::endl;

    BitVector activeLevels(fineLevelIdx+1, true);

    std::vector<std::shared_ptr<LocalOscAssembler>> localAssemblers;
    std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>> localIntersectionAssemblers;

    localAssemblers.resize(activeLevels.size());
    localIntersectionAssemblers.resize(activeLevels.size());

    for (size_t i=0; i<activeLevels.size(); i++) {
        localAssemblers[i] = std::make_shared<LocalOscAssembler>(oscDataHierarchy[i]->data(), oscDataHierarchy[i]->mapper());

        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(i);
        std::vector<std::shared_ptr<LocalInterfaceAssembler>>& levelLocalIntersectionAssemblers = localIntersectionAssemblers[i];
        levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size());
        for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) {
            const int k = minCantorResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j);
            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
            levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
        }
    }

    typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
    typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
    MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;

    MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
    for (size_t i=0; i<preconditioner.size()-1; i++) {
        LevelPatchPreconditioner& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i);

        levelFaultPreconditioner.setPatchDepth(patchDepth);
        levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
    }
    preconditioner.build();

    std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
    std::cout << std::endl << std::endl;

    // set initial iterate
    VectorType initialX;
    DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis);
    coarseFineTransfer.prolong(coarseSol, initialX);

    VectorType rhs(fineRhs);

    DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis);

    // solve
    OscCGSolver<MatrixType, VectorType, DGBasis>
    solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner,
                    maxIterations, solverTolerance, &exactEnergyNorm,
                    Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN)


    solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineResolution) + ".vec";

    solver.check();
    solver.preprocess();
    solver.solve();

    std::cout << std::endl;
    std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl; 

    std::cout.rdbuf(coutbuf); //reset to standard output again

    std::cout << "Problem " << problemCount << " done!" << std::endl;
    problemCount++;
    }
} catch (Dune::Exception e) {

    std::cout << e << std::endl;

}
}