Postprocess and visualize CHT flow over plate results

flow-over-heated-plate/README.md

@@ -57,6 +57,23 @@ You may use additional filters, such as the Calculator and the Plot Over Line, t
 
 ![graph](images/tutorials-flow-over-heated-plate-graph-result.png)
 
+### Compare results
+
+First generate the output for each run by adding export to the participant `Solid` in `precice-config.xml`:
+
+```
+<participant name="Solid">
+  <export:vtk directory="preCICE-output" />
+  <use-mesh name="Fluid-Mesh" from="Fluid" />
+  <use-mesh name="Solid-Mesh" provide="yes" />
+  ...
+</participant>
+```
+
+After that running a case from this tutorial will export data into `solid-*/preCICE-output`. To visualize and compare these results run `python3 plot-final-interface-temperature.py` (You can install the required python packages by running `pip3 install -r plot-final-interface-temperature-requirements.txt`). This will plot the dimensionless temperature `theta = (T-300)/(310-300)` (with `T` being the temperature) across the coupling interface, i.e. where the solid and the fluid meet and exchange heat. The x-axis shows the x coordinate and the y-axis the dimensionless temperature `theta` at the interface. If you want to exclude certain cases, simply comment out the corresponding lines in the script. For reference see below:
+
+![](images/tutorials-flow-over-heated-plate-results-comparison.png)
+
 ## References
 
 [1]  M. Vynnycky, S. Kimura, K. Kanev, and I. Pop. Forced convection heat transfer from a flat plate: the conjugate problem. International Journal of Heat and Mass Transfer, 41(1):45 – 59, 1998.

flow-over-heated-plate/plot-final-interface-temperature-requirements.txt

@@ -0,0 +1,3 @@
+vtk
+numpy
+matplotlib

flow-over-heated-plate/plot-final-interface-temperature.py

@@ -0,0 +1,63 @@
+import vtk
+from matplotlib import pyplot as plt
+import numpy as np
+
+
+def vtk_to_dict(case):
+    vtkFileName = "solid-{}/preCICE-output/Fluid-Mesh-Solid.final.vtk".format(case)
+    # read the vtk file as an unstructured grid
+    reader = vtk.vtkUnstructuredGridReader()
+    reader.SetFileName(vtkFileName)
+    reader.ReadAllVectorsOn()
+    reader.ReadAllScalarsOn()
+    reader.Update()
+
+    # obtain the data
+    data = reader.GetOutput()
+    n_data = data.GetPointData().GetNumberOfTuples()
+
+    name = "Temperature"
+    data_names = []
+    i = 0
+    max_i = data.GetPointData().GetNumberOfArrays()
+    while i < max_i:
+        this_data_name = data.GetPointData().GetArray(i).GetName()
+        data_names.append(this_data_name)
+        if(this_data_name == name):
+            data_id = i
+            break
+        i += 1
+
+    data_dict = {}
+
+    if not data_id:
+        raise Exception(
+            "For file {} name {} not found. Only the following names are available: {}. "
+            "Aborting!".format(vtkFileName, name, data_names))
+    for i in range(n_data):
+        data_dict[data.GetPoint(i)] = data.GetPointData().GetArray(data_id).GetValue(i)
+
+    return data_dict
+
+
+cases = []
+cases.append('fenics')
+cases.append('openfoam')
+cases.append('nutils')
+
+case_labels = {'fenics': 'OpenFOAM-FEniCS', 'openfoam': 'OpenFOAM-OpenFOAM', 'nutils': 'OpenFOAM-Nutils', }
+styles = [':', '-', '--']
+colors = ['r', 'b', 'g', 'k']
+i = 0
+
+for case in cases:
+    case_data = vtk_to_dict(case)
+    x, t = [p[0] for p in case_data.keys()], np.array(list(case_data.values()))
+    theta = (t - 300) / (310 - 300)
+    plt.plot(x, theta, colors[i % 4] + styles[i % 3], label=case_labels[case])
+    i += 1
+
+plt.ylabel("Theta")
+plt.xlabel("x-coordinate along coupling interface")
+plt.legend()
+plt.show()