Select Git revision
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;
}
}