From 3d89cd18a963422bd2c1109aa7f6210d48726449 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Thu, 19 Mar 2020 19:39:03 +0100
Subject: [PATCH] add python script for convergence test
---
src/CMakeLists.txt | 3 +-
src/convergenceTest.py | 111 +++++++++++++++++++++++++++++++++++++++++
src/dg.cc | 21 +++++++-
src/mmdg.cc | 5 ++
4 files changed, 137 insertions(+), 3 deletions(-)
create mode 100644 src/convergenceTest.py
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7f1cc3f..0733f27 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 0000000..f542831
--- /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 4b55a71..4eb8d25 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 4c96269..47c2e83 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();
--
GitLab