cantorconvergence.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::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("/home/mi/podlesny/software/dune/faultnetworks/src/cantorfaultnetworks/cantorconvergence.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 minResolution = problemParameters.get<int>("minResolution");
    const int maxResolution = problemParameters.get<int>("maxResolution");
    const int minCantorResolution = problemParameters.get<int>("minCantorResolution");

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

    const int maxGridN = std::pow(2, maxResolution);


    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minResolution) + "_" + std::to_string(maxResolution) + "_" + std::to_string(minCantorResolution) +  "_" + std::to_string((int) penaltyFactor) + ".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)) - minResolution;

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

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

    const int maxLevelIdx = maxResolution - minResolution;

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

    // build multilevel cantor fault network
    CantorFaultFactory<GridType> faultFactory(levelToCantorLevel, minResolution, maxLevelIdx, maxCantorLevel);
    const InterfaceNetwork<GridType>& interfaceNetwork = faultFactory.interfaceNetwork();
    for (size_t i=0; i<interfaceNetwork.size(); i++) {
        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
        networkWriter.write(interfaceNetwork.levelInterfaceNetwork(i), true);
    }

    std::cout << "Interface network built successfully!" << std::endl << std::endl;

    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 direct solutions                    ---
    // ----------------------------------------------------------------
    std::cout << std::endl;
    std::cout << "Computing direct solutions!" << std::endl;

    std::vector<VectorType> solutions;
    solutions.resize(maxLevelIdx+1);
    std::vector<double> discErrors;
    discErrors.resize(maxLevelIdx+1);

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

    const LevelInterfaceNetwork<GV>& maxLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(maxLevelIdx);
    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> maxGlobalAssembler(maxLevelInterfaceNetwork);
    LocalOscAssembler maxLocalAssembler(oscDataHierarchy[maxLevelIdx]->data(), oscDataHierarchy[maxLevelIdx]->mapper());

    std::vector<std::shared_ptr<LocalInterfaceAssembler>> maxLocalInterfaceAssemblers(maxLevelInterfaceNetwork.size());
    for (size_t j=0; j<maxLocalInterfaceAssemblers.size(); j++) {
        //const double pen = penaltyFactor;
        const double pen = penaltyFactor*std::pow(2.0, minResolution + maxLevelInterfaceNetwork.getIntersectionsLevels().at(j));
        maxLocalInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
    }

    maxGlobalAssembler.assemble(maxLocalAssembler, maxLocalInterfaceAssemblers, l2FunctionalAssembler);
    maxGlobalAssembler.solve();

    solutions[maxLevelIdx] = maxGlobalAssembler.getSol();
    discErrors[maxLevelIdx] = 0;

    const DGBasis& maxBasis = maxGlobalAssembler.basis();

    for (int i=maxLevelIdx-1; i>=0; i--) {
        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(i);
        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> globalAssembler(levelInterfaceNetwork);
        LocalOscAssembler localAssembler(oscDataHierarchy[i]->data(), oscDataHierarchy[i]->mapper());

        std::vector<std::shared_ptr<LocalInterfaceAssembler>> localInterfaceAssemblers(levelInterfaceNetwork.size());
        for (size_t j=0; j<localInterfaceAssemblers.size(); j++) {
            //const double pen = penaltyFactor;
            const double pen = penaltyFactor*std::pow(2.0, minResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j));
            localInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
        }

        globalAssembler.assemble(localAssembler, localInterfaceAssemblers, l2FunctionalAssembler);
        globalAssembler.solve();

        const VectorType& sol = globalAssembler.getSol();
        const DGBasis& basis = globalAssembler.basis();

        DGMGTransfer<DGBasis> levelTransfer(basis, maxBasis);
        levelTransfer.prolong(sol, solutions[i]);

        VectorType discError(solutions[maxLevelIdx]);
        discError -= solutions[i];
        discErrors[i] = discError.two_norm();
    }

    // print results
    std::cout << std::endl << std::endl;
    std::cout << "level  discError"<< std::endl;
    std::cout << "----------------" << std::endl;
    for (size_t i=0; i<discErrors.size(); i++) {
        std::cout << i << "      " << discErrors[i] << std::endl;
    }

    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;
}
}