From 7f31a52e5bcbc0b3e81f56c4de8a43605eb25c56 Mon Sep 17 00:00:00 2001 From: Abdullah Mujahid Date: Sun, 25 Feb 2024 21:57:38 +0100 Subject: [PATCH] Add python code in results.dox --- examples/step-3/doc/results.dox | 236 ++++++++++++++++++++++++++++++- examples/step-3/plotting.py | 239 -------------------------------- 2 files changed, 234 insertions(+), 241 deletions(-) delete mode 100644 examples/step-3/plotting.py diff --git a/examples/step-3/doc/results.dox b/examples/step-3/doc/results.dox index 8a9460ce1b..f904c64eb6 100644 --- a/examples/step-3/doc/results.dox +++ b/examples/step-3/doc/results.dox @@ -308,7 +308,7 @@ data_file.write_dataset("mean_value",mean_value);

Using R and ggplot2 to generate plots

-@note Alternatively, one could use the associated python script `plotting.py`. +@note Alternatively, one could use the python code in the next subsection. The data put into %HDF5 files above can then be used from scripting languages for further postprocessing. In the following, let us show @@ -435,7 +435,7 @@ This is now going to look as follows: -For plotting the converge curves we need to re-run the C++ code multiple times with different values for n_refinement_steps +For plotting the convergence curves we need to re-run the C++ code multiple times with different values for n_refinement_steps starting from 1. Since every file only contains a single data point we need to loop over them and concatenate the results into a single vector. @code{.r} @@ -496,3 +496,235 @@ to zero: + +

Using python to generate plots

+ +In this section we discuss the postprocessing of the data stored in %HDF5 files using the "python" programming language. +The necessary packages to import are +@code{.py} +import numpy as np # to work with multidimensional arrays +import h5py # to work with %HDF5 files + +import pandas as pd # for data frames +import matplotlib.pyplot as plt # plotting +from matplotlib.patches import Polygon + +from scipy.interpolate import griddata # interpolation function +from matplotlib import cm # for colormaps + +@endcode +The %HDF5 solution file corresponding to `refinement = 5` can be opened as +@code{.py} +refinement = 5 +filename = "solution_%d.h5" % refinement +file = h5py.File(filename, "r") +@endcode +The following prints out the tables in the %HDF5 file +@code{.py} +for key, value in file.items(): + print(key, " : ", value) +@endcode +which prints out +@code{.unparsed} +cells : +mean_value : +nodes : +point_value : +solution : +@endcode +There are $(32+1)\times(32+1) = 1089$ nodes. +The coordinates of these nodes $(x,y)$ are stored in the table `nodes` in the %HDF5 file. +There are a total of $32\times 32 = 1024$ cells. The nodes which make up each cell are +marked in the table `cells` in the %HDF5 file. + +Next, we extract the data into multidimensional arrays +@code{.py} +nodes = np.array(file["/nodes"]) +cells = np.array(file["/cells"]) +solution = np.array(file["/solution"]) + +x, y = nodes.T +@endcode + +The following stores the $x$ and $y$ coordinates of each node of each cell in one flat array. +@code{.py} +cell_x = x[cells.flatten()] +cell_y = y[cells.flatten()] +@endcode +The following tags the cell ids. Each four entries correspond to one cell. +Then we collect the coordinates and ids into a data frame +@code{.py} +n_cells = cells.shape[0] +cell_ids = np.repeat(np.arange(n_cells), 4) +meshdata = pd.DataFrame({"x": cell_x, "y": cell_y, "ids": cell_ids}) +@endcode +The data frame looks +@code{.unparsed} +print(meshdata) + + x y ids +0 0.00000 0.00000 0 +1 0.03125 0.00000 0 +2 0.03125 0.03125 0 +3 0.00000 0.03125 0 +4 0.03125 0.00000 1 +... ... ... ... +4091 0.93750 1.00000 1022 +4092 0.96875 0.96875 1023 +4093 1.00000 0.96875 1023 +4094 1.00000 1.00000 1023 +4095 0.96875 1.00000 1023 + +4096 rows × 3 columns +@endcode + +To plot the mesh, we loop over all cells and connect the first and last node to complete the polygon +@code{.py} +fig, ax = plt.subplots() +ax.set_aspect("equal", "box") +ax.set_title("grid at refinement level #%d" % refinement) + +for cell_id, cell in meshdata.groupby(["ids"]): + cell = pd.concat([cell, cell.head(1)]) + plt.plot(cell["x"], cell["y"], c="k") +@endcode +Alternatively one could use the matplotlib built-in Polygon class +@code{.py} +fig, ax = plt.subplots() +ax.set_aspect("equal", "box") +ax.set_title("grid at refinement level #%d" % refinement) +for cell_id, cell in meshdata.groupby(["ids"]): + patch = Polygon(cell[["x", "y"]], facecolor="w", edgecolor="k") + ax.add_patch(patch) +@endcode + +To plot the solution, we first create a finer grid and then interpolate the solution values +into the grid and then plot it. +@code{.py} +nx = int(np.sqrt(n_cells)) + 1 +nx *= 10 +xg = np.linspace(x.min(), x.max(), nx) +yg = np.linspace(y.min(), y.max(), nx) + +xgrid, ygrid = np.meshgrid(xg, yg) +solution_grid = griddata((x, y), solution.flatten(), (xgrid, ygrid), method="linear") + +fig = plt.figure() +ax = fig.add_subplot(1, 1, 1) +ax.set_title("solution at refinement level #%d" % refinement) +c = ax.pcolor(xgrid, ygrid, solution_grid, cmap=cm.viridis) +fig.colorbar(c, ax=ax) + +plt.show() +@endcode + +To check the convergence of `mean_value` and `point_value` +we loop over data of all refinements and store into vectors mean_values and mean_values +@code{.py} +mean_values = np.zeros((8,)) +point_values = np.zeros((8,)) +dofs = np.zeros((8,)) + +for refinement in range(1, 9): + filename = "solution_%d.h5" % refinement + file = h5py.File(filename, "r") + mean_values[refinement - 1] = np.array(file["/mean_value"])[0] + point_values[refinement - 1] = np.array(file["/point_value"])[0] + dofs[refinement - 1] = np.array(file["/solution"]).shape[0] +@endcode + +Following is the plot of mean_values on `log-log` scale +@code{.py} +mean_error = np.abs(mean_values[1:] - mean_values[:-1]) +plt.loglog(dofs[:-1], mean_error) +plt.grid() +plt.xlabel("#DoFs") +plt.ylabel("mean value error") +plt.show() +@endcode + +Following is the plot of point_values on `log-log` scale +@code{.py} +point_error = np.abs(point_values[1:] - point_values[:-1]) +plt.loglog(dofs[:-1], point_error) +plt.grid() +plt.xlabel("#DoFs") +plt.ylabel("point value error") +plt.show() +@endcode + +A python package which mimicks the `R` ggplot2 (which is based on specifying the grammar of the graphics) is +plotnine. +@code{.py} +We need to import the following from the plotnine package +from plotnine import ( + ggplot, + aes, + geom_raster, + geom_polygon, + geom_line, + labs, + scale_x_log10, + scale_y_log10, + ggtitle, +) +@endcode +Then plot the mesh meshdata dataframe +@code{.py} +plot = ( + ggplot(meshdata, aes(x="x", y="y", group="ids")) + + geom_polygon(fill="white", colour="black") + + ggtitle("grid at refinement level #%d" % refinement) +) + +print(plot) +@endcode +Collect the solution into a dataframe +@code{.py} +colordata = pd.DataFrame({"x": x, "y": y, "solution": solution.flatten()}) +@endcode +Plot of the solution +@code{.py} +plot = ( + ggplot(colordata, aes(x="x", y="y", fill="solution")) + + geom_raster(interpolate=True) + + ggtitle("solution at refinement level #%d" % refinement) +) + +print(plot) +@endcode + +Collect the convergence data into a dataframe +@code{.py} +convdata = pd.DataFrame( + {"dofs": dofs[:-1], "mean_value": mean_error, "point_value": point_error} +) + +@endcode +Following is the plot of mean_values on `log-log` scale +@code{.py} +plot = ( + ggplot(convdata, mapping=aes(x="dofs", y="mean_value")) + + geom_line() + + labs(x="#DoFs", y="mean value error") + + scale_x_log10() + + scale_y_log10() +) + +plot.save("mean_error.pdf", dpi=60) +print(plot) +@endcode + +Following is the plot of point_values on `log-log` scale +@code{.py} +plot = ( + ggplot(convdata, mapping=aes(x="dofs", y="point_value")) + + geom_line() + + labs(x="#DoFs", y="point value error") + + scale_x_log10() + + scale_y_log10() +) + +plot.save("point_error.pdf", dpi=60) +print(plot) +@endcode diff --git a/examples/step-3/plotting.py b/examples/step-3/plotting.py deleted file mode 100644 index 0001b14f9a..0000000000 --- a/examples/step-3/plotting.py +++ /dev/null @@ -1,239 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# jupytext_version: 1.16.1 -# kernelspec: -# display_name: Python 3 (ipykernel) -# language: python -# name: python3 -# --- - -# %% [markdown] -# # H5 files and plotting - -# %% -import numpy as np -import h5py - -import pandas as pd -import matplotlib.pyplot as plt -from matplotlib.patches import Polygon - -from scipy.interpolate import griddata -from matplotlib import cm - -# %% -refinement = 5 -filename = "solution_%d.h5" % refinement -file = h5py.File(filename, "r") - -# %% -for key, value in file.items(): - print(key, " : ", value) - -# %% -nodes = np.array(file["/nodes"]) -cells = np.array(file["/cells"]) -solution = np.array(file["/solution"]) - -x, y = nodes.T - -# %% -nodes - -# %% -cells - -# %% -solution - -# %% -x - -# %% -y - -# %% -cell_x = x[cells.flatten()] -cell_y = y[cells.flatten()] - -# %% -cell_x - -# %% -cell_y - -# %% -cell_ids = np.repeat(np.arange(cells.shape[0]), 4) -cell_ids - -# %% -n_cells = cells.shape[0] -n_cells - -# %% -meshdata = pd.DataFrame({"x": cell_x, "y": cell_y, "ids": cell_ids}) - -# %% -meshdata - -# %% -fig, ax = plt.subplots() -ax.set_aspect("equal", "box") -ax.set_title("grid at refinement level #%d" % refinement) - -for cell_id, cell in meshdata.groupby(["ids"]): - cell = pd.concat([cell, cell.head(1)]) - plt.plot(cell["x"], cell["y"], c="k") - -# %% [markdown] -# An alternative is to use the `matplotlib` built-in `Polygon` class - -# %% -fig, ax = plt.subplots() -ax.set_aspect("equal", "box") -ax.set_title("grid at refinement level #%d" % refinement) -for cell_id, cell in meshdata.groupby(["ids"]): - patch = Polygon(cell[["x", "y"]], facecolor="w", edgecolor="k") - ax.add_patch(patch) - - -# %% [markdown] -# ## A color plot of the solution - -# %% -nx = int(np.sqrt(n_cells)) + 1 -nx *= 10 -xg = np.linspace(x.min(), x.max(), nx) -yg = np.linspace(y.min(), y.max(), nx) - -xgrid, ygrid = np.meshgrid(xg, yg) -solution_grid = griddata((x, y), solution.flatten(), (xgrid, ygrid), method="linear") - -fig = plt.figure() -ax = fig.add_subplot(1, 1, 1) -ax.set_title("solution at refinement level #%d" % refinement) -c = ax.pcolor(xgrid, ygrid, solution_grid, cmap=cm.viridis) -fig.colorbar(c, ax=ax) - -plt.show() - -# %% [markdown] -# ## Convergence - -# %% -mean_values = np.zeros((8,)) -point_values = np.zeros((8,)) -dofs = np.zeros((8,)) - -for refinement in range(1, 9): - filename = "solution_%d.h5" % refinement - file = h5py.File(filename, "r") - mean_values[refinement - 1] = np.array(file["/mean_value"])[0] - point_values[refinement - 1] = np.array(file["/point_value"])[0] - dofs[refinement - 1] = np.array(file["/solution"]).shape[0] - - -# %% -mean_values - -# %% -point_values - -# %% -dofs - -# %% -mean_error = np.abs(mean_values[1:] - mean_values[:-1]) -plt.loglog(dofs[:-1], mean_error) -plt.grid() -plt.xlabel("#DoFs") -plt.ylabel("mean value error") -plt.show() - -# %% -point_error = np.abs(point_values[1:] - point_values[:-1]) -plt.loglog(dofs[:-1], point_error) -plt.grid() -plt.xlabel("#DoFs") -plt.ylabel("point value error") -plt.show() - -# %% [markdown] -# # Using *python* equivalent of *ggplot* of R -# -# **[plotnine](https://plotnine.readthedocs.io/en/v0.12.4/)** - -# %% -# !pip install plotnine - -# %% -from plotnine import ( - ggplot, - aes, - geom_raster, - geom_polygon, - geom_line, - labs, - scale_x_log10, - scale_y_log10, - ggtitle, -) # noqa: E402 - -# %% -plot = ( - ggplot(meshdata, aes(x="x", y="y", group="ids")) - + geom_polygon(fill="white", colour="black") - + ggtitle("grid at refinement level #%d" % refinement) -) - -print(plot) - -# %% -colordata = pd.DataFrame({"x": x, "y": y, "solution": solution.flatten()}) -colordata - -# %% -plot = ( - ggplot(colordata, aes(x="x", y="y", fill="solution")) - + geom_raster(interpolate=True) - + ggtitle("solution at refinement level #%d" % refinement) -) - -print(plot) - -# %% -convdata = pd.DataFrame( - {"dofs": dofs[:-1], "mean_value": mean_error, "point_value": point_error} -) - -convdata - -# %% -plot = ( - ggplot(convdata, mapping=aes(x="dofs", y="mean_value")) - + geom_line() - + labs(x="#DoFs", y="mean value error") - + scale_x_log10() - + scale_y_log10() -) - -plot.save("mean_error.pdf", dpi=60) -print(plot) - -# %% -plot = ( - ggplot(convdata, mapping=aes(x="dofs", y="point_value")) - + geom_line() - + labs(x="#DoFs", y="point value error") - + scale_x_log10() - + scale_y_log10() -) - -plot.save("point_error.pdf", dpi=60) -print(plot) -- 2.39.5