diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7f1cc3f048a06cc669ac8b46b789986ed1318198..0733f2789cfe13704721b3793ae5163ddf4c2cc5 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -14,4 +14,5 @@ add_executable("mmdg-2d" mmdg.cc)
 target_compile_definitions("mmdg-2d" PRIVATE GRIDDIM=2)
 target_link_dune_default_libraries("mmdg-2d")
 
-dune_symlink_to_source_files(FILES grids parameterMMDG.ini parameterDG.ini)
+dune_symlink_to_source_files(FILES grids parameterMMDG.ini parameterDG.ini
+  convergenceTest.py)
diff --git a/src/convergenceTest.py b/src/convergenceTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..f5428315ee2475dd899a36594b568b6e8ac27d0d
--- /dev/null
+++ b/src/convergenceTest.py
@@ -0,0 +1,111 @@
+import numpy as np
+from os.path import join, exists
+from subprocess import call
+from math import log
+
+
+# writes a dgf file for a grid with N elements per directions and dimension dim
+def writeDGF(N, dim):
+     # open dgf file
+    file = open(join("grids", "cube" + str(dim) + "d.dgf"), "w")
+
+    file.write("DGF\n\nInterval\n")
+
+    # write lower left corner
+    for i in range(0, dim):
+        file.write(str(0) + "\t")
+    file.write("\n")
+
+    # write upper right corner
+    for i in range(0, dim):
+        file.write(str(1) + "\t")
+    file.write("\n")
+
+    # write number of grid elements
+    for i in range(0, dim):
+        file.write(str(N) + "\t")
+    file.write("\n#\n\nSimplex\n#")
+
+    # close dgf file
+    file.close()
+
+
+# reads error, run time and number of elements, i.e. the output from the
+# excuted discontinuous Galerkin scheme
+def readData():
+    file = open("convergenceData", "r")
+
+    # read data
+    error = float(file.readline())
+    timeTotal = float(file.readline())
+    numberOfElements = float(file.readline())
+
+    file.close()
+
+    # return data
+    return np.array([error, timeTotal, numberOfElements])
+
+
+# === main ===
+
+# create directory for the plots if necessary
+if not exists("plots"):
+    os.mkdir("plots")
+
+# number of runs from which the run times are taken as mean value
+numberOfMeans = 3
+
+# numbers of grid elements per direction
+numElemsPerDir = [1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
+
+dimensions = [1, 2]
+
+for dim in dimensions:
+    # storage for errors, run times and largestEdgeLength
+    # column 0: error
+    # column 1: total run time
+    # column 2: largest edge length
+    # different rows correspond to different numbers of grid elements
+    data = np.zeros((len(numElemsPerDir), 3))
+
+    for i in range(len(numElemsPerDir)):
+        # write dgf file
+        writeDGF(numElemsPerDir[i], dim)
+
+        # take mean value over numberOfMeans runs
+        for j in range(numberOfMeans):
+            # run executable
+            call(join(".", "dg-" + str(dim) + "d"))
+
+            # store data
+            data[i, :] += readData()
+        # end for j
+    # end for i
+
+    # divide by numberOfMeans
+    data /= numberOfMeans
+
+    # determine experimental order of convergence (EOC)
+    EOC = np.empty(len(numElemsPerDir) - 1)
+
+    for i in range(len(numElemsPerDir) - 1):
+        if data[i, 0] != 0:
+            EOC[i] = log(data[i, 0] / data[i+1, 0]) / \
+                log(data[i+1, 2] / data[i, 2])
+        else:
+            EOC[i] = -1.0*float('inf')
+    # end for i
+
+    print("\n\n-------------------------------------------------------------\n")
+    # print EOCs
+    print("EOCs (" + str(dim) + "d):\t")
+    print(EOC)
+    print("\n\n")
+
+    # write data to file
+    np.savetxt(join("plots", "dgAnalysis_" + str(dim) + "d"), data)
+    np.savetxt(join("plots", "dgEOC_" + str(dim) + "d"), EOC)
+# end for dim
+
+# remove output files
+os.remove("convergenceData")
diff --git a/src/dg.cc b/src/dg.cc
index 4b55a715a6c225af670d92c8deba0fa348978b58..4eb8d25c238fcc38bc045c1e8be15d2518448730 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -5,6 +5,7 @@
 #include <sys/stat.h>
 #include <sys/types.h>
 #include <type_traits>
+#include <iostream>
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/parametertree.hh>
@@ -87,10 +88,26 @@ int main(int argc, char** argv)
 
     dg(mu);
 
-    //print total run time
+    //error and total run time
+    const double error = dg.computeL2error();
     const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
+
+    //write error, runtime and number of grid elements to file
+    std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
+
+    if (errorFile.is_open())
+    {
+      errorFile << error << "\n" << timeTotal << "\n" << gridView.size(0);
+      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: " << dg.computeL2error() << ".\n";
+      " Seconds.\nL2-error: " << error << ".\n";
 
     return 0;
   }
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 4c96269004925504ff7e23707f9db571ffb5d6c1..47c2e830efe66b890e75adcb269603fdfa4bb92f 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -76,6 +76,11 @@ int main(int argc, char** argv)
         "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();