diff --git a/lecture/mhs/groundwater/CMakeLists.txt b/lecture/mhs/groundwater/CMakeLists.txt index ff1a2640f5b1ecdd316bf308980d0f3b6d8b635e..0439e9f3b11581c6889e6a9978b77992b59b8f0b 100644 --- a/lecture/mhs/groundwater/CMakeLists.txt +++ b/lecture/mhs/groundwater/CMakeLists.txt @@ -1,4 +1,5 @@ add_input_file_links() +dune_symlink_to_source_files(FILES "contourPlots.py") dumux_add_test(NAME groundwater SOURCES groundwater.cc diff --git a/lecture/mhs/groundwater/contourPlots.py b/lecture/mhs/groundwater/contourPlots.py new file mode 100644 index 0000000000000000000000000000000000000000..63438a53d3f2b81ee3a6c12856e927319cbabade --- /dev/null +++ b/lecture/mhs/groundwater/contourPlots.py @@ -0,0 +1,78 @@ +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() diff --git a/lecture/mhs/groundwater/groundwaterproblem.hh b/lecture/mhs/groundwater/groundwaterproblem.hh index 1b421a96865f3a84fec3ce6a9575c2c3778929be..7396bb8d1a38f3443a9df4fec35a5f8fcb9be549 100644 --- a/lecture/mhs/groundwater/groundwaterproblem.hh +++ b/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: