Merge branch 'feature/add-contour-plot' into 'master'

Add contour plot script for groundwater simulation

See merge request !206

lecture/mhs/groundwater/CMakeLists.txt

 add_input_file_links()
+dune_symlink_to_source_files(FILES "contourPlots.py")
 
 dumux_add_test(NAME groundwater
               SOURCES groundwater.cc

lecture/mhs/groundwater/contourPlots.py

+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+#Instructions: First execute the `groundwater` test, a file called "contourFile.txt" will be written to the respective build-cmake folder. Within the build-cmake/.../groundwater folder you can then run the python script via e.g. `python3 contourPlots.py` to plot the contour and streamlines of the given problem.
+
+plt.rcParams.update({'font.size': 10})
+
+#extra function for reading in the Dumux parameter file
+def read_parameters(filename):
+    parameters = {}
+    with open(filename, 'r') as file:
+        for line in file:
+            line = line.strip()
+            # Ignore empty lines and comments and only consider lines which contain the character "="
+            if line and not line.startswith("#") and not line.startswith("[") and "=" in line:
+                name, value = line.split('=')
+                #parameters[name.strip()] = float(value.strip())
+                parameters[name.strip()] = value.strip()
+    return parameters
+
+#-------READING IN DATA FROM FILE-------
+#read in the necessary parameters from the Dumux input file
+filename = "groundwater.input" #choose the correct name for the dumux parameter file
+parameters = read_parameters(filename)
+cells = parameters.get("Cells").split()
+nx = int(cells[0])
+ny = int(cells[1])
+#read in the data which should be plotted
+myData = pd.read_csv('contourFile.txt') #choose the correct name of the data file
+
+#-------PREPARING DATA FOR CONTOUR PLOT-------
+xArray = myData['x'].values #read in a pandas data column and immediately convert it to a numpy array
+yArray = myData['y'].values
+hArray = myData['h'].values
+#reshape input from a vector to a matrix
+x = xArray.reshape(ny,nx)
+xLinspace = x[0] # each row is the same, take the first
+y = yArray.reshape(ny,nx)
+yLinspace = np.asarray([i[0] for i in y]) # only take first entries of each vector as the structure of y is transposed compared to that of x
+#reshape output from a vector to a matrix
+h = hArray.reshape(ny,nx)
+#create a meshgrid from x and y
+X,Y = np.meshgrid(xLinspace,yLinspace)
+
+#-------CONTOUR PLOT-------
+#contour plot
+numberLevels=20
+fig, ax = plt.subplots()
+contourPlot = ax.contour(X, Y, h, levels=numberLevels, colors='black', linewidths=0.5)
+
+#contourPlot = ax.contour(X, Y, h, levels=[4,12,20,28,36,40,41,42,43,44,45,46,47,48,49,49.5,50])
+ax.clabel(contourPlot, inline=True, fontsize=7)
+ax.set_xlabel('x')
+ax.set_ylabel('y')
+fig.savefig('isolineplot.png', dpi=300)
+
+#-------PREPARING DATA FOR STREAMLINE PLOT-------
+#read in the data for the streamlines
+vxArray = myData['v_x'].values #read in a pandas data column and immediately convert it to a numpy array
+vyArray = myData['v_y'].values
+#reshape velocity fields from a vector to a matrix
+vx = vxArray.reshape(ny,nx)
+vy = vyArray.reshape(ny,nx)
+
+#-------STREAMLINE PLOT-------
+#streamline plot over contour plot
+fig2, ax2 = plt.subplots()
+#contour plot
+contourPlot = ax2.contour(X, Y, h, levels=numberLevels, colors='black')
+ax2.clabel(contourPlot, inline=True, fontsize=10)
+#streamline plot
+streamlinePlot = ax2.streamplot(xLinspace, yLinspace, vx, vy, color='red', density=1, linewidth=0.5)
+ax2.set_xlabel('x')
+ax2.set_ylabel('y')
+fig2.savefig('streamlineplot.png', dpi=300)
+
+plt.show()

lecture/mhs/groundwater/groundwaterproblem.hh

@@ -415,6 +415,10 @@ public:
         std::cout << "         x          y          h           v_x           v_y"<<std::endl;
         std::cout << "------------------------------------------------------------"<<std::endl;
 
+        std::ofstream contourFile;
+        contourFile.open("contourFile.txt");
+        contourFile << "x,y,h,v_x,v_y" << std::endl;
+
         for (const auto& element : elements(this->gridView()))
         {
             int cellIndex = this->variables().index(element);
@@ -432,8 +436,11 @@ public:
             y = element.geometry().center()[1];
 
             printf("%10.4g %10.4g %10.4g %13.4g %13.4g\n", x, y, piezo, v_x, v_y);
+
+            contourFile << std::setprecision(15) << x << "," << y << "," << piezo << "," << v_x << "," << v_y << std::endl;
         }
 
+        contourFile.close();
     }
 
 private: