From: Wolfgang Bangerth Date: Tue, 7 Jan 2020 16:45:11 +0000 (-0700) Subject: Update the discussion of HDF5 and R scripting in step-3. X-Git-Tag: v9.2.0-rc1~709^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9250%2Fhead;p=dealii.git Update the discussion of HDF5 and R scripting in step-3. This is a really nice introduction to what can be done by putting data into hdf5 files and then postprocessing them in R. I'm really only updating stylistic issues here. --- diff --git a/examples/step-3/doc/results.dox b/examples/step-3/doc/results.dox index 9b04473491..d67030f1eb 100644 --- a/examples/step-3/doc/results.dox +++ b/examples/step-3/doc/results.dox @@ -222,20 +222,32 @@ suggestions: Again, the difference between two adjacent values goes down by about a factor of four, indicating convergence as ${\cal O}(h^2)$. -

Using the %HDF5 Library to output the solution and additional data

-%HDF5 is a commonly used format that can be read by many scripting languages (e.g. R or Python). Here are some -insights on what is possible. -

Changing the output from .gnuplot to .h5:

+ + + +

Using %HDF5 to output the solution and additional data

+ +%HDF5 is a commonly used format that can be read by many scripting +languages (e.g. R or Python). It is not difficult to get deal.II to +produce some %HDF5 files that can then be used in external scripts to +postprocess some of the data generated by this program. Here are some +ideas on what is possible. + + +

Changing the output to .h5

+ To fully make use of the automation we first need to introduce a private variable for the number of global refinement steps unsigned int n_refinement_steps , which will be used for the output filename. -In make_grid() we replace triangulation.refine_global(5); with +In make_grid() we then replace triangulation.refine_global(5); with @code n_refinement_steps = 5; triangulation.refine_global(n_refinement_steps); @endcode -The deal.II library has two different %HDF5 bindings, one in the HDF5 namespace -and another one in DataOut. -Although the HDF5 deal.II binding supports both, serial and MPI, the %HDF5 DataOut binding +The deal.II library has two different %HDF5 bindings, one in the HDF5 +namespace (for interfacing to general-purpose data files) +and another one in DataOut (specifically for writing files for the +visualization of solutions). +Although the HDF5 deal.II binding supports both serial and MPI, the %HDF5 DataOut binding only supports parallel output. For this reason we need to initialize an MPI communicator with only one processor. This is done by adding the following code. @@ -246,7 +258,8 @@ int main(int argc, char* argv[]) ... } @endcode -Next we change the output routine as described in the namespace documentation. +Next we change the `Step3::output_results()` output routine as +described in the DataOutBase namespace documentation: @code const std::string filename_h5 = "solution_" + std::to_string(n_refinement_steps) + ".h5"; DataOutBase::DataOutFilterFlags flags(true, true); @@ -254,49 +267,76 @@ DataOutBase::DataOutFilter data_filter(flags); data_out.write_filtered_data(data_filter); data_out.write_hdf5_parallel(data_filter, filename_h5, MPI_COMM_WORLD); @endcode +The resulting file can then be visualized just like the VTK file that +the original version of the tutorial produces; but, since %HDF5 is a +more general file format, it can also easily be processed in scripting +languages for other purposes. -

Adding the point value and the mean (see extension above) into the .h5 file:

-After outputing the solution, the file can be opened again to include more datasets. -This allows us to keep all the necessary information of our experiment in a single result file, -which can then be read and processed by some postprocessing script.
-Have a further look at HDF5::Group::write_dataset() for further information on the possible output options.
-First we include the necessary header into our file. +

Adding the point value and the mean (see extension above) into the .h5 file

+ +After outputing the solution, the file can be opened again to include +more datasets. This allows us to keep all the necessary information +of our experiment in a single result file, which can then be read and +processed by some postprocessing script. +(Have a look at HDF5::Group::write_dataset() for further +information on the possible output options.) + +To make this happen, we first include the necessary header into our file: @code -"" +#include @endcode Adding the following lines to the end -of our output routine adds the additional information to our %HDF5 file. +of our output routine adds the information about the value of the +solution at a particular point, as well as the mean value of the +solution, to our %HDF5 file: @code HDF5::File data_file(filename_h5, HDF5::File::FileAccessMode::open, MPI_COMM_WORLD); Vector point_value(1); -point_value[0] = VectorTools::point_value(dof_handler, solution, Point<2>(1./3, 1./3)); +point_value[0] = VectorTools::point_value(dof_handler, solution, + Point<2>(1./3, 1./3)); data_file.write_dataset("point_value", point_value); Vector mean_value(1); -mean_value[0] = VectorTools::compute_mean_value(dof_handler, QGauss<2>(fe.degree + 1), solution, 0); +mean_value[0] = VectorTools::compute_mean_value(dof_handler, + QGauss<2>(fe.degree + 1), + solution, 0); data_file.write_dataset("mean_value",mean_value); @endcode -

Using R and ggplot2 to generate PDFs with plots of the grid, the solution and convergence curves:

+ + +

Using R and ggplot2 to generate plots

+ +The data put into %HDF5 files above can then be used from scripting +languages for further postprocessing. In the following, let us show +how this can, in particular, be done with the +R +programming language, a widely used language in statistical data +analysis. (Similar things can also be done in Python, for example.) If you are unfamiliar with R and ggplot2 you could check out the data carpentry course on R - here.
+here. Furthermore, since most search engines struggle with searches of the form "R + topic", -we recommend using the specializes service RSeek instead. -The most prominent difference between R and other languages is that the assignment operator (a = 5) is interchangable with -(a <- 5). As the latter is considered standard we will use it in our examples as well. -To open the .h5 file in R you have to install the rhdf5 package, which is a part of the Bioconductor package. -
First we will include all necessary packages and have a look at how the data is structured in our file. -@code -library(rhdf5) #library for handling HDF5 files -library(ggplot2) #main plotting library -library(grDevices) #needed for output to PDF -library(viridis) #contains good colormaps for sequential data +we recommend using the specializes service RSeek instead. + +The most prominent difference between R and other languages is that +the assignment operator (`a = 5`) is typically written as +`a <- 5`. As the latter is considered standard we will use it in our examples as well. +To open the `.h5` file in R you have to install the rhdf5 package, which is a part of the Bioconductor package. + +First we will include all necessary packages and have a look at how the data is structured in our file. +@code{.r} +library(rhdf5) # library for handling HDF5 files +library(ggplot2) # main plotting library +library(grDevices) # needed for output to PDF +library(viridis) # contains good colormaps for sequential data + refinement <- 5 h5f <- H5Fopen(paste("solution_",refinement,".h5",sep="")) print(h5f) @endcode This gives the following output -@code +@code{.unparsed} HDF5 FILE name / filename @@ -308,41 +348,55 @@ filename 3 point_value H5I_DATASET FLOAT 1 4 solution H5I_DATASET FLOAT x 1089 @endcode -The datasets can be accessed by h5f\$name. The function dim(h5f\$cells) gives us -the dimensions of our cell matrix. We can see the following three matrices, as well as the two additional -data points we added. +The datasets can be accessed by h5f\$name. The function +dim(h5f\$cells) gives us the dimensions of the matrix +that is used to store our cells. +We can see the following three matrices, as well as the two +additional data points we added.
  • cells: a 4x1024 matrix that stores the (C++) vertex indices for each cell
  • nodes: a 2x1089 matrix storing the position values (x,y) for our cell vertices
  • solution: a 1x1089 matrix storing the values of our solution at each vertex
-Now we can use this data to generate our various plots. Plotting with ggplot2 usually splits into two steps. -At first the data needs to be manipulated and added to a data.frame. -After that a ggplot object is constructed and manipulated by adding plot elements to it. +Now we can use this data to generate various plots. Plotting with ggplot2 usually splits into two steps. +At first the data needs to be manipulated and added to a data.frame. +After that, a ggplot object is constructed and manipulated by adding plot elements to it. nodes and cells contain all the information we need to plot our grid. -The following code wraps all the data into one dataframe for plotting our grid -@code -#counting in R starts at 1 instead of 0 so we need to increment all vertex indices -cell_ids <- h5f$cells+1 -#store the x and y positions of each vertex in one big vector in a cell by cell fashion -# (every 4 entries belong to one cell) +The following code wraps all the data into one dataframe for plotting our grid: +@code{.r} +# Counting in R starts at 1 instead of 0, so we need to increment all +# vertex indices by one: +cell_ids <- h5f$cells+1 + +# Store the x and y positions of each vertex in one big vector in a +# cell by cell fashion (every 4 entries belong to one cell): cells_x <- h5f$nodes[1,][cell_ids] cells_y <- h5f$nodes[2,][cell_ids] -#construct a vector that stores the matching cell by cell grouping (1,1,1,1,2,2,2,2,...) + +# Construct a vector that stores the matching cell by cell grouping +# (1,1,1,1,2,2,2,2,...): groups <- rep(1:ncol(cell_ids),each=4) -#finally put everything into one dataframe + +# Finally put everything into one dataframe: meshdata <- data.frame(x = cells_x, y = cells_y, id = groups) @endcode -With the finished dataframe we have everything we need to plot our grid -@code -pdf (paste("grid_",refinement,".pdf",sep=""),width = 5,height = 5) #Open new PDF file -plt <- ggplot(meshdata,aes(x=x,y=y,group=id)) #construction of our plot object, at first only data -plt <- plt + geom_polygon(fill="white",colour="black") #actual plotting of the grid as polygons + +With the finished dataframe we have everything we need to plot our grid: +@code{.r} +pdf (paste("grid_",refinement,".pdf",sep=""),width = 5,height = 5) # Open new PDF file +plt <- ggplot(meshdata,aes(x=x,y=y,group=id)) # Construction of our plot + # object, at first only data + +plt <- plt + geom_polygon(fill="white",colour="black") # Actual plotting of the grid as polygons plt <- plt + ggtitle(paste("grid at refinement level #",refinement)) -print(plt) #show the current state of the plot/ add it to the pdf -dev.off() #close PDF file + +print(plt) # Show the current state of the plot/add it to the pdf +dev.off() # Close PDF file @endcode + +The contents of this file then look as follows (not very exciting, but +you get the idea):
@@ -350,22 +404,26 @@ dev.off() #close PDF file
-To make a 2D pseudocolor plot of our solution we will use geom_raster . + +We can also visualize the solution itself, and this is going to look +more interesting. +To make a 2D pseudocolor plot of our solution we will use geom_raster. This function needs a structured grid, i.e. uniform in x and y directions. Luckily our data at this point is structured in the right way. -A description of the data manipulation needed for unstructured grids will be added to step-6. -The following code plots a pseudocolor representation of our surface into a new PDF. -@code -pdf (paste("pseudocolor_",refinement,".pdf",sep=""),width = 5,height = 4.2) #Open new PDF file +The following code plots a pseudocolor representation of our surface into a new PDF: +@code{.r} +pdf (paste("pseudocolor_",refinement,".pdf",sep=""),width = 5,height = 4.2) # Open new PDF file colordata <- data.frame(x = h5f$nodes[1,],y = h5f$nodes[2,] , solution = h5f$solution[1,]) plt <- ggplot(colordata,aes(x=x,y=y,fill=solution)) plt <- plt + geom_raster(interpolate=TRUE) plt <- plt + scale_fill_viridis() plt <- plt + ggtitle(paste("solution at refinement level #",refinement)) -print(plt) #show the current state of the plot/ add it to the pdf -dev.off() #close PDF file -H5Fclose(h5f) #close HDF5 file + +print(plt) +dev.off() +H5Fclose(h5f) # Close the HDF5 file @endcode +This is now going to look as follows:
@@ -374,38 +432,43 @@ H5Fclose(h5f) #close HDF5 file
  • + For plotting the converge 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 -n_ref <- 8 #maximum refinement level for which results are existing -#first we initiate all vectors with the results of the first level -h5f <- H5Fopen("solution_1.h5") -dofs <- dim(h5f$solution)[2] -mean <- h5f$mean_value +@code{.r} +n_ref <- 8 # Maximum refinement level for which results are existing + +# First we initiate all vectors with the results of the first level +h5f <- H5Fopen("solution_1.h5") +dofs <- dim(h5f$solution)[2] +mean <- h5f$mean_value point <- h5f$point_value H5Fclose(h5f) + for (reflevel in 2:n_ref) { - h5f <- H5Fopen(paste("solution_",reflevel,".h5",sep="")) - dofs <- c(dofs,dim(h5f\$solution)[2]) - mean <- c(mean,h5f\$mean_value) + h5f <- H5Fopen(paste("solution_",reflevel,".h5",sep="")) + dofs <- c(dofs,dim(h5f\$solution)[2]) + mean <- c(mean,h5f\$mean_value) point <- c(point,h5f\$point_value) H5Fclose(h5f) } @endcode As we are not interested in the values themselves but rather in the error compared to a "exact" solution we will assume our highest refinement level to be that solution and omit it from the data. -@code -#calculate the error w.r.t. our maximum refinement step -mean_error <- abs(mean[1:n_ref-1]-mean[n_ref]) +@code{.r} +# Calculate the error w.r.t. our maximum refinement step +mean_error <- abs(mean[1:n_ref-1]-mean[n_ref]) point_error <- abs(point[1:n_ref-1]-point[n_ref]) -#remove the highest value from our DoF data -dofs <- dofs[1:n_ref-1] + +# Remove the highest value from our DoF data +dofs <- dofs[1:n_ref-1] convdata <- data.frame(dofs = dofs, mean_value= mean_error, point_value = point_error) @endcode Now we have all the data available to generate our plots. -As is usual we will plot our results on a log-log scale. +It is often useful to plot errors on a log-log scale, which is +accomplished in the following code: @code pdf (paste("convergence.pdf",sep=""),width = 5,height = 4.2) plt <- ggplot(convdata,mapping=aes(x = dofs, y = mean_value)) @@ -413,12 +476,18 @@ plt <- plt+geom_line() plt <- plt+labs(x="#DoFs",y = "mean value error") plt <- plt+scale_x_log10()+scale_y_log10() print(plt) + plt <- ggplot(convdata,mapping=aes(x = dofs, y = point_value)) plt <- plt+geom_line() plt <- plt+labs(x="#DoFs",y = "point value error") plt <- plt+scale_x_log10()+scale_y_log10() print(plt) + dev.off() + +This results in the following plot that shows how the errors in the +mean value and the solution value at the chosen point nicely converge +to zero: @endcode