Skip to content
Snippets Groups Projects
Select Git revision
  • 2e8fd2849be7a50e401296882dfe42d40f1acdae
  • master default protected
  • release/1.0
3 results

mmdg.cc

Blame
  • mmdg.cc 4.45 KiB
    #ifdef HAVE_CONFIG_H
    # include "config.h"
    #endif
    
    #include <sys/stat.h>
    #include <sys/types.h>
    #include <type_traits>
    #include <string>
    
    #include <dune/common/exceptions.hh>
    #include <dune/common/parametertree.hh>
    #include <dune/common/parametertreeparser.hh>
    
    #include <dune/grid/common/mcmgmapper.hh>
    #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
    #include <dune/grid/io/file/dgfparser/dgfug.hh>
    
    #include <dune/mmesh/mmesh.hh>
    
    #include <dune/mmdg/mmdg.hh>
    #include <dune/mmdg/clockhelper.hh>
    #include <dune/mmdg/problems/coupleddgproblem.hh>
    #include <dune/mmdg/problems/mmdgproblem1.hh>
    #include <dune/mmdg/problems/mmdgproblem2.hh>
    #include <dune/mmdg/problems/mmdgproblem3.hh>
    #include <dune/mmdg/problems/mmdgproblem4.hh>
    
    int main(int argc, char** argv)
    {
      try
      {
        using Clock = std::chrono::steady_clock;
    
        //start run time measurement
        const Clock::time_point timeInitial = Clock::now();
    
        static constexpr int dim = GRIDDIM;
    
        using Grid = Dune::MovingMesh<dim>;
        using GridFactory = Dune::GmshGridFactory<Grid, /*implicit=*/false>;
        using GridView = typename Grid::LeafGridView;
        using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
        using IGrid = typename Grid::InterfaceGrid;
        using IGridView = typename IGrid::LeafGridView;
        using IMapper =
          typename Dune::MultipleCodimMultipleGeomTypeMapper<IGridView>;
        using Coordinate = Dune::FieldVector<double, dim>;
        using Problem = CoupledDGProblem<Coordinate>;
        using MMDG = MMDG<GridView, Mapper, IGridView, IMapper, Problem>;
    
        //get parameters from file "parameterDG.ini"
        Dune::ParameterTree pt;
        Dune::ParameterTreeParser::readINITree("parameterMMDG.ini", pt);
    
        //assume constant penalty parameter //TODO
        const double mu = std::stod(pt["mu"]);
        const double xi = std::stod(pt["xi"]); //coupling parameter
        const double aperture = std::stod(pt["d"]);
    
        Problem* problem; //problem to be solved
    
        std::string gridType;
    
        //determine problem type
        if (pt["problem"] == "mmdg1")
        {
          problem = new MMDGProblem1<Coordinate>(aperture);
          gridType = "mmdg2_5e-2";
        }
        else if (pt["problem"] == "mmdg2")
        {
          problem = new MMDGProblem2<Coordinate>();
          gridType = "mmdg2_5e-2";
        }
        else if (pt["problem"] == "mmdg3")
        {
          problem = new MMDGProblem3<Coordinate>(aperture);
          gridType = "mmdg3_5e-2";
        }
        else if (pt["problem"] == "mmdg4")
        {
          problem = new MMDGProblem4<Coordinate>(aperture);
          gridType = "mmdg2_5e-2";
        }
        else
        {
          DUNE_THROW(Dune::Exception,
            "Invalid problem type or parameter 'problem' not specified in file "
            "parameterDG.ini (use 'mmdg2' or TODO)");
        }
    
        if (pt.hasKey("gridfile")) //non-default grid type
        {
          gridType = pt["gridfile"];
        }
    
        //create a grid from .dgf file
        GridFactory gridFactory( "grids/" + gridType + ".msh" );
        const Grid& grid = *gridFactory.grid();
        const GridView& gridView = grid.leafGridView();
        const IGrid& iGrid = grid.interfaceGrid();
        const IGridView& iGridView = iGrid.leafGridView();
    
        //element mapper for the leaf grid views
        const Mapper mapper(gridView, Dune::mcmgElementLayout());
        const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
    
        MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
        mmdg(mu, xi);
    
        //error and total run time
        const double error = mmdg.computeL2error();
        const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
    
        //determine maximum edge length
        double largestEdgeLength = 0.;
        for (const auto& edge : edges(gridView))
        {
          largestEdgeLength = std::max(largestEdgeLength, edge.geometry().volume());
        }
    
        //write error, runtime and maximum edge length
        std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
    
        if (errorFile.is_open())
        {
          errorFile << error << "\n" << timeTotal << "\n" << largestEdgeLength;
          errorFile.close();
        }
        else
        {
          std::cout << "Unable to write error data to file.\n";
        }
    
        //print total run time and error
        std::cout << "Finished with a total run time of " << timeTotal <<
          " Seconds.\nL2-error: " << error << ".\n";
    
        return 0;
      }
      catch (Dune::Exception& e)
      {
        std::cerr << "Dune reported error: " << e.what() << std::endl;
      }
      catch (std::exception& e)
      {
        std::cerr << "Error: " << e.what() << std::endl;
      }
      catch (...)
      {
        std::cerr << "Unknown exception thrown!" << std::endl;
      }
    }