add python script for convergence test

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)

src/convergenceTest.py

+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")

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

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();