From e04bf6737435642d0f9e7add7c43ae465f449859 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Fri, 3 Apr 2020 15:34:38 +0200
Subject: [PATCH] add separate EOCs for bulk and interface
---
dune/mmdg/mmdg.hh | 13 +++++------
src/mmdg.cc | 30 +++++++++++++++++--------
src/mmdgAnalysis.py | 54 +++++++++++++++++++++++++++++++--------------
3 files changed, 65 insertions(+), 32 deletions(-)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 32c5953..b9a9741 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -48,9 +48,11 @@ public:
writeVTKoutput();
}
- //compute L2 error using quadrature rules
+ //compute L2 errors using quadrature rules,
+ //returns {error(bulk), error(interface), error(total)},
+ //the total error is calculated using the norm
// norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) )
- Scalar computeL2error(const int quadratureOrder = 5)
+ const std::array<Scalar, 3> computeL2error(const int quadratureOrder = 5)
{
const Scalar bulkError = Base::computeL2error(quadratureOrder);
Scalar interfaceError(0.0);
@@ -87,11 +89,8 @@ public:
Base::applyQuadratureRule(iGeo, qrule, qlambda);
}
- //print bulk and interface error
- std::cout << "ERROR bulk: " << bulkError << "\t ERROR interface: " <<
- sqrt(interfaceError) << "\n";
-
- return sqrt(bulkError * bulkError + interfaceError);
+ return {bulkError, sqrt(interfaceError),
+ sqrt(bulkError * bulkError + interfaceError)};
}
const InterfaceGridView& iGridView_;
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 35ba464..706813e 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -131,23 +131,33 @@ int main(int argc, char** argv)
MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
mmdg(mu0, xi);
- //error and total run time
- const double error = mmdg.computeL2error();
+ //errors and total run time
+ const std::array<double, 3> errors = mmdg.computeL2error();
const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
- //determine maximum edge length
- double largestEdgeLength = 0.;
+ //determine maximum edge lengths
+ double hMax = 0.0; //maximum edge length (bulk)
+ double hMaxInterface = 0.0; //maximum edge length (interface)
+
for (const auto& edge : edges(gridView))
{
- largestEdgeLength = std::max(largestEdgeLength, edge.geometry().volume());
+ if (grid.isInterface(edge))
+ {
+ hMaxInterface = std::max(hMaxInterface, edge.geometry().volume());
+ }
+ else
+ {
+ hMax = std::max(hMax, edge.geometry().volume());
+ }
}
- //write error, runtime and maximum edge length
+ //write errors, runtime and maximum edge lengths
std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
if (errorFile.is_open())
{
- errorFile << error << "\n" << timeTotal << "\n" << largestEdgeLength;
+ errorFile << errors[0] << "\n" << errors[1] << "\n" << errors[2] << "\n"
+ << hMax << "\n" << hMaxInterface << "\n" << timeTotal;
errorFile.close();
}
else
@@ -155,9 +165,11 @@ int main(int argc, char** argv)
std::cout << "Unable to write error data to file.\n";
}
- //print total run time and error
+ //print total run time and errors
std::cout << "Finished with a total run time of " << timeTotal <<
- " Seconds.\nL2-error: " << error << ".\n\n";
+ " Seconds.\nL2-error (bulk): \t" << errors[0] <<
+ "\nL2-error (interface): \t" << errors[1] << "\nL2-error (total): \t"
+ << errors[2] << ".\n\n";
return 0;
}
diff --git a/src/mmdgAnalysis.py b/src/mmdgAnalysis.py
index c5b972e..a78c888 100644
--- a/src/mmdgAnalysis.py
+++ b/src/mmdgAnalysis.py
@@ -41,14 +41,36 @@ def readData():
file = open("convergenceData", "r")
# read data
- error = float(file.readline())
+ errorBulk = float(file.readline())
+ errorInterface = float(file.readline())
+ errorTotal = float(file.readline())
+ hMax = float(file.readline())
+ hMaxInterface = float(file.readline())
timeTotal = float(file.readline())
- numberOfElements = float(file.readline())
file.close()
# return data
- return np.array([error, timeTotal, numberOfElements])
+ return np.array([errorBulk, errorInterface, errorTotal, hMax, \
+ hMaxInterface, timeTotal])
+
+
+# determines experimental order of convergence (EOC),
+# data is assumed to be a numpy array with
+# first column: error,
+# second column: maximum edge length
+def EOCloop (data):
+ EOC = np.empty(data.shape[0] - 1)
+
+ for i in range(len(EOC)):
+ if data[i, 0] != 0:
+ EOC[i] = log(data[i, 0] / data[i+1, 0]) / \
+ log(data[i, 1] / data[i+1, 1])
+ else:
+ EOC[i] = -1.0*float('inf')
+ # end for i
+
+ return EOC
# === main ===
@@ -72,7 +94,7 @@ dim = 2
# column 1: total run time
# column 2: largest edge length
# different rows correspond to different numbers of grid elements
-data = np.zeros((len(gridfiles), 3))
+data = np.zeros((len(gridfiles), 6))
for i in range(len(gridfiles)):
# write dgf file
@@ -92,25 +114,25 @@ for i in range(len(gridfiles)):
data /= numberOfMeans
# determine experimental order of convergence (EOC)
-EOC = np.empty(len(gridfiles) - 1)
-
-for i in range(len(gridfiles) - 1):
- if data[i, 0] != 0:
- EOC[i] = log(data[i, 0] / data[i+1, 0]) / \
- log(data[i, 2] / data[i+1, 2])
- else:
- EOC[i] = -1.0*float('inf')
-# end for i
+EOCbulk = EOCloop(data[:, [0,3]])
+EOCinterface = EOCloop(data[:, [1,4]])
+EOCtotal = \
+ EOCloop(np.column_stack([data[:,2], np.maximum(data[:,3], data[:,4])]))
# print EOCs
print("\n-------------------------------------------------------------\n")
-print("EOCs (" + str(dim) + "d):\t")
-print(EOC)
+print("EOCs (" + str(dim) + "d):")
+print("\nbulk:")
+print(EOCbulk)
+print("\ninterface:")
+print(EOCinterface)
+print("\ntotal:")
+print(EOCtotal)
print("\n=============================================================\n\n")
# write data to file
np.savetxt(join("plots", "mmdgAnalysis_" + str(dim) + "d"), data)
-np.savetxt(join("plots", "mmdgEOC_" + str(dim) + "d"), EOC)
+np.savetxt(join("plots", "mmdgEOC_" + str(dim) + "d"), EOCtotal)
# remove output files
remove("convergenceData")
--
GitLab