Mesh generation in CFD problem: An approach using Object Oriented Programming in R

 

 

Background

I am currently working in understanding the water evaporation dynamics from soil. I read the paper of Kondo et al. The paper develop a model for calculating volumetric water content, specific humidity and temperature by mean of the solution of a system of PDE.

In order to solve the PDE, a solution approximation in finite differences is implemented. The algorith implementation involves a pre-processor module the which is conformed by fluid dynamics, boundary conditions and grid generation.

Provided boundary conditions, fluid dynamics equations and geometry, the grid can be generated. The objective was to develop a procedure for generating a mesh taking into account changes in geometry and boundary conditions

Solution

I developed a class of type mesh, consisting of four attributes ( number of elements (ix), the compact matrix (mat), the number of nodes (nodes), and dimension order (dimension) ) and one method (print).

Several functions are required to set the class constructor: calculate the increment in every space dimension (x, y, or z) (xyz_dist); calculate the time increment to guarantee numerical stability (t_dis); set the value of boundary conditions (bound_cond), calculate the total number of element in mess (sum1toi)

Mesh generator requires an list as input parameter, consisting of dimension, xyz_dist, t_dist, boun_cond, and Courant sigma value. An instance of mesh class can be obtaining by passing the inparam list as argument, and result in matrix format can be obtained by mean of print.mesh method.

The mesh constructor and print method was based upon the class writing code read in the book The Art of R Programming by Norman Matloff

Mesh generation is tested for the 1D case. In the code shown below, the mesh for time intervals (column 1), 1D space (column 2), and boundary condition value (column 3). The mesh can be extended up to 3D by configuring the len_nod matrix in inparam list.

xyz_dist <- function(len_nod=matrix(data = c(c(rep(0.1,3)), c(rep(10,3))),nrow =2,ncol=3, byrow = TRUE) ){
  len <- len_nod[1,]
  innod <- len_nod[2,]
  rtrn <- vector(mode="numeric", length = 3)
  for (i in 1:length(rtrn)){
    rtrn[i]<-len[i]/(innod[i]-1)
  }
  rtrn
}



t_dist <- function(xyz_d=vector(mode = "numeric",length = 3), sigma=0.5){
  rtrn <- vector(mode = "numeric", length = 1)
  for (i in 1:length(rtrn)){
    rtrn[i] <- sigma * max(xyz_d)
  }
  rtrn
}

boun_cond <- function(theta=0.1, q=0.1, t=298.15){
  return(1) #Dirichlet
}

sum1toi <- function(i,j) return(sum(i)*(j+2)) #i is mesh nodes vector, j is dimension

inparam <- list(di=1
                ,len_nod=matrix(data = c(c(rep(0.1,3)), c(rep(10,3))),nrow =2,ncol=3, byrow = TRUE)
                ,geom_dist = function(len_nod=matrix(data=0.0, nrow= 2,ncol=3) ){
                  len <- len_nod[1,]
                  innod <- len_nod[2,]
                  rtrn <- vector(mode="numeric", length = 3)
                  for (i in 1:length(rtrn)){
                    rtrn[i]<-len[i]/(innod[i]-1)
                  }
                  rtrn
                }
                ,time_dist = function(xyz_d=vector(mode = "numeric",length = 3), sigma=0.5){
                  rtrn <- vector(mode = "numeric", length = 1)
                  for (i in 1:length(rtrn)){
                    rtrn[i] <- sigma * max(xyz_d)
                  }
                  rtrn
                }
                ,bc = function(theta=0.1, q=0.1, t=298.15){
                  return(3.5e-05*(theta/0.4)**15) #Dirichlet
                }
                ,sigma= 0.5)

mesh <- function(inparam){
  # debug(inparam[[4]])
  rtrn <- list()
  class(rtrn) <- "cdf_mesh"
  spa_di <- inparam[[1]]
  lnods <- (inparam[[2]])
  n <- sum1toi(max(lnods),spa_di)
  rtrn$mat <- vector(length= n)
  rtrn$ix <- (1:n)
  rtrn$dimensions <- spa_di
  rtrn$nodes <- lnods[2,]
  rtrn$n <- sum1toi(max(lnods[2,]),spa_di)
  fun_spa_dist <- inparam[[3]]
  fun_t_dist <- inparam[[4]]
  fun_bc <- inparam[[5]]
  sigma <- inparam[[6]]
  t_i <- 0
  for (i in seq(1,n,spa_di+2)){
    t_i <- t_i + fun_t_dist(fun_spa_dist(lnods),sigma)
    rtrn$mat[i]<- t_i
  }
  xyz_step <- 0
  for (j in 1:spa_di){
    for (i in seq(j+1,n,spa_di+2)){
      xyz_step <- xyz_step +fun_spa_dist(lnods)[j]
      rtrn$mat[i]<- xyz_step
    }
    for (i in seq(spa_di+2,n,spa_di+2)){
      xyz_step <- xyz_step +fun_spa_dist(lnods)[j]
      theta <- 0.1/xyz_step
      rtrn$mat[i]<- fun_bc(theta=theta)
    }
  }
  rtrn
}

print.mesh <- function(mesh){
  m_n <- length(mesh$mat)
  geomdim <- mesh$dimension
  n_col <- mesh$dimension+2
  nodes <- max(mesh$nodes)
  mat_vector <- mesh$mat
  n <- mesh$n
  full_mesh <- matrix(nrow=nodes, ncol = n_col)
  for (i in 1:(geomdim+2)) {
    full_mesh[1:nodes,i] <- mat_vector[seq(i,m_n,n_col)]
  }
  return(full_mesh)
}

# Let's test the clas
mesh1 <- mesh(inparam)

# Print mesh
print.mesh(mesh1)
##              [,1]       [,2]         [,3]
##  [1,] 0.005555556 0.01111111 1.6066287663
##  [2,] 0.011111111 0.02222222 0.4355992520
##  [3,] 0.016666667 0.03333333 0.1311159516
##  [4,] 0.022222222 0.04444444 0.0431405106
##  [5,] 0.027777778 0.05555556 0.0153262862
##  [6,] 0.033333333 0.06666667 0.0058211136
##  [7,] 0.038888889 0.07777778 0.0023446175
##  [8,] 0.044444444 0.08888889 0.0009947598
##  [9,] 0.050000000 0.10000000 0.0004420810
## [10,] 0.055555556 0.11111111 0.0002048122
# Plot the boundary condition against x
mesh_1D <- print.mesh(mesh1)
plot(mesh_1D[,2:3])


 

Use of ggplot and knitr in a rmarkdown to report results of activated sludge simulation

I am currently attending an online specialization course on water sanitation modelling at Unesco-IHE.
On the first assignment, I consider writing the report straight away in R using Rmarkdown PDF template, rather than jumping from Excel to Word by modelling, formatting and writing.
I found the experience of using Rmarkdown time saving and enjoyable, as I was capable of performing a regression on the fly as I was writing the report and produce neat and nice visualizations and tables with ggplot and knitr::kable.
Rmarkdown stand out with several capabilities:

  • Writing nice formulas like chemical reaction kinetics,
  • representing the chemical processes (endogenous respiration and air mass transfer in Petersen matrix), by using the kable function
  • Embedding R code for data import, visualization, transformation and modelling using readxl, ggplot, and linear regression.
  • Embedding jpg images in the repor using the command ![]+image filepath.

You can see the rmd and html sample files in the following urls:

Enjoy!

Database dependencies graph visualization



Database dependencies graph visualization

The journey from JSON file to graph format, and visualization with Gephi

Rafael Ventura

2017-05-06

Goal

To enable Business Intelligence team visualising legacy database table dependencies on designing a migration process

Scope

As this an exploratory study, R was chosen to read, process and transform a JSON file with database table dependencies into a file with a graph format.

Procedure

1.- Read the JSON file

To read the JSON file, I used the rjon library to read from a text file having the data in a JSON structure.

library("rjson")
json_file <- "jsonTV.txt"
json_data <- fromJSON(file=json_file)

2.- Transform JSON structure into a data frame

For transforming JSON structure into a data frame, I developed a function populateData.

populateDFData <- function(json_data,df_data) {
  rowIndex <- 0
  for (i in 1:length(json_data)) {

  lengthImport <- length(json_data[[i]]$imports)
  if (lengthImport>0) {
    for (j in 1:lengthImport) {
      rowIndex <- rowIndex+1
      df_data[rowIndex,"name"]<-json_data[[i]]$name
      df_data[rowIndex,"size"]<-json_data[[i]]$size
      import <- json_data[[i]]$imports
      df_data[rowIndex,"imports"]<-import[[j]]
    } } else{
      rowIndex <- rowIndex + 1
      df_data[rowIndex,"name"]<-json_data[[i]]$name
      df_data[rowIndex,"size"]<-json_data[[i]]$size
      df_data[rowIndex,"imports"]<-"<NA>"
    } 
    
  }
}

I declare a data frame whose number of rows were calculated using a function calculateLength, that defines a record for each occurrence of the tuple (name, import), where name is the table and import its dependencies.

calculateLength <- function(jsonObject) {
  dfLength <- 0
  for (i in 1:length(jsonObject)) {
    importSize <- length(jsonObject[[i]]$import)
    if (importSize == 0) {incrm <- 1} else {incrm<-importSize}
    dfLength <- dfLength + incrm
  }
  return(dfLength)
}
dfLength <- calculateLength(json_data) 
df_data <- data.frame(name=rep("",dfLength), 
                      size=rep("",dfLength), 
                      imports=rep("",dfLength))
df_data$name<- as.character(df_data$name)
df_data$size<- as.character(df_data$size)
df_data$imports<- as.character(df_data$imports)

Data frame is populated

df_data <- populateDFData(json_data = json_data, 
                          df_data = df_data)

3.- Export data to file with graph format

As we wanted to visualise the graph with Gephi, data input needs to contain headers named Source and Target. A data frame declaration is necessary. Besides, we need to filter only tables containing the expression “AMS” and dependencies with one or more tables. Hence, it was necessary to define two filter conditions. A graph data frame was defined using the function graph.data.frame, and this object was exported to a graphml file using the function write.graph.

yaml {r} library("igraph") gephi_df <- data.frame(Source=rep("",nrow(df_data)), Target=rep("",nrow(df_data)), Label=rep("",nrow(df_data))) gephi_df$Source <- df_data$imports gephi_df$Target <- df_data$name gephi_df$Label <- df_data$name gephi_df.filter1 <- grepl("AMS",gephi_df$Target) gephi_df.filter2 <- !grepl("<NA>",gephi_df$Source) gephi_df.export <- gephi_df[gephi_df.filter1&gephi_df.filter2,] gephi_df.g <- graph.data.frame(gephi_df.export[,c("Source", "Target")]) write.graph(gephi_df.g, file="tvia-objects-dependencies.graphml", format = "graphml")



Indian Water Quality Visualization Using Heatmaps


Recently I download a file with Indian data quality data in Kaggle (https://www.kaggle.com/venkatramakrishnan/india-water-quality-data). I decided to explore data visualization package to learn how to make neat and tidy visualization using heat maps.File has the following structure:

water.quality<-readRDS("water.quality.rds")
head(water.quality)
##       State.Name     District.Name     Block.Name  Panchayat.Name
## 1 ANDHRA PRADESH EAST GODAVARI(04) PRATHIPADU(10)   GOKAVARAM(04)
## 2 ANDHRA PRADESH EAST GODAVARI(04) PRATHIPADU(10)   GOKAVARAM(04)
## 3 ANDHRA PRADESH EAST GODAVARI(04) PRATHIPADU(10) GAJJANAPUDI(06)
## 4 ANDHRA PRADESH EAST GODAVARI(04) PRATHIPADU(10) GAJJANAPUDI(06)
## 5 ANDHRA PRADESH EAST GODAVARI(04) PRATHIPADU(10)  CHINTALURU(10)
## 6 ANDHRA PRADESH EAST GODAVARI(04) PRATHIPADU(10)       ELURU(16)
##               Village.Name                      Habitation.Name
## 1           VANTHADA(014 )           VANTHADA(0404410014010400)
## 2     PANDAVULAPALEM(022 )     PANDAVULAPALEM(0404410022010400)
## 3         G. KOTHURU(023 )         G. KOTHURU(0404410023010600)
## 4        GAJJANAPUDI(029 )        GAJJANAPUDI(0404410029010600)
## 5         CHINTALURU(028 )         CHINTALURU(0404410028011000)
## 6 P. JAGANNADHAPURAM(035 ) P. JAGANNADHAPURAM(0404410035011600)
##   Quality.Parameter     Year
## 1          Salinity 1/4/2009
## 2          Fluoride 1/4/2009
## 3          Salinity 1/4/2009
## 4          Salinity 1/4/2009
## 5          Salinity 1/4/2009
## 6          Fluoride 1/4/2009

I decided to organize the work in two R scripts: Data preparation and data visualization My first attempt was to use an advance heat maps function (heatmap.2) I found in blog post by Joseph Rickert (http://www.r-bloggers.com/r-for-more-powerful-clustering). As heatmap uses clustering, a matrix of at least 2 rows by 2 columns was required. I added five custom fields, each one with the name of chemical existing in each sample. As heatmaps requires data to be of class numeric matrix, I had to prepare date in order to ensure that a chemical ocurrence in matrix was marked as 1 (existing) or 0 (absent).I gave up using all data as object generated by scaled function was 116GB in size, and went through using samples of state names. Graphics obtained were cluttered and I was not able to add the year in visualization

source("water-quality-data-preparation-v1.R")
source("water-quality-exploration-v1.R")
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess

I searched then in my collection of R-Blogger post and found a recipe of R code for visualization posted by Dikesh Jariwala (https://www.r-bloggers.com/7-visualizations-you-should-learn-in-r/) and found a code for plotting heatmaps using ggplot. I adapted the code to a discrete variable and obtained a more neat and tidy heatmap.

library(devtools)
library(ggplot2)
water.quality1 <- readRDS("water.quality1.rds")
# Create a function for a data frame, from which the plot is drawn, and uploaded to Github.
source("https://github.com/rventuradiaz/indian-water-quality/raw/master/water-quality-data-preparation-v4.R")
wq_df <- wq1_df(water.quality1) #use the function and obtain the wrangled data frame.
# Plot the heatmap
ggplot(wq_df, aes(Year, VillageName))+
  geom_raster(aes(fill = QualityParameter))+
  labs(title ="Heat Map", x = "Year", y = "Village Name" )+
  scale_fill_discrete(name = "QualityParameter")+
  facet_grid(facets = . ~ PanchayatName  )+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size= 5), strip.text.x = element_text(angle = 90, vjust = 0.5, size = 8))

Source data: Indian Water Quality see https://www.kaggle.com/venkatramakrishnan/india-water-quality-data Reference: R for more powerful clustering (see http://www.r-bloggers.com/r-for-more-powerful-clustering/)