R/extend_sweep_mat.R

In laurin-f/pkg.WWM: What the Package Does (One Line, Title Case)

# library(dplyr)
# library(stringr)
# library(ggplot2)
# filename <- "freeSoil_anisotropy_sweep_2DS.txt"

#' function to expand sweep matrix generated by COMSOL
#'
#' @param filename filename to COMSOL sweep txt
#' @param byout resolution of extended matrix
#'
#' @return
#' @export
#'
#' @examples
extend_sweep_mat <- function(filename = "freeSoil_anisotropy_sweep_2DS.txt",
                             byout = 5e-8) {
  n_DS <- stringr::str_extract(filename,pattern = "\\d+(?=DS)") %>% as.numeric()
  # #Parameter die in der Datei gesweept wurden
  pars <- c(paste0("DS_",1:n_DS),"injection_rate")
  
  sweep_long <- read_sweep(filename, format = "long")
  
  par_vals <- lapply(sweep_long[-7],unique)
  
  
  
  
  
  extend_long <- data.frame()
  
  for (i in unique(sweep_long$z)) {
    print(paste("z =", i))
    for (j in par_vals$injection_rate) {
      print(paste("injection_rate =", j))
      
      sweep_i <- subset(sweep_long, z == i & injection_rate == j)
      cols <-  paste0("DS_",1:n_DS)
      formula <- as.formula(paste(cols, collapse = "~"))
      value <- "CO2_mol_per_m3"
      sweep_mat <- reshape2::acast(sweep_i, formula, value.var = value)
      
      
      xyz <- lapply(dimnames(sweep_mat), as.numeric)
      names(xyz) <- cols
      
      xyz_seq <-
        (lapply(xyz, function(x)
          seq(min(x), max(x), by = byout)))
      
      xyz_out <- expand.grid(xyz_seq)
      
      
      
      approx_vec <- e1071::interpolate(xyz_out, sweep_mat)
      
      
      approx_list <-
        cbind(
          z = i,
          injection_rate = j,
          xyz_out,
          CO2_mol_per_m3 = approx_vec
        )
      extend_long <- rbind(extend_long, approx_list)
      
    }
  }
  
  
  sweep_extend <- extend_long
  for(i in paste0("DS_",1:n_DS)){
    sweep_extend[,i] <- paste0(i,"=",extend_long[,i])
  }
  extend_wide <- tidyr::pivot_wider(sweep_extend,names_from = matches("DS"),names_sep=", ",values_from=CO2_mol_per_m3) %>% as.data.frame()
  rm(sweep_extend)
  
  
  mod_inj_rates <- unique(extend_wide$injection_rate)
  sweep_colnames <- colnames(extend_wide[-(1:2)])
  tiefen <- 150-unique(extend_wide$z)
  
  extend_arr <- array(dim = c(nrow(extend_wide)/2,ncol(extend_wide)-2,2))
  for(i in 1:2){
    extend_arr[,,i] <- as.matrix(subset(extend_wide[,-(1)],injection_rate == unique(extend_wide$injection_rate)[i]))[,-1]
  }
  
  dimnames(extend_arr) <- list(tiefen,sweep_colnames,mod_inj_rates)
  
  ###############
  #save
  ##############
  save(extend_long,file=paste0(comsolpfad,str_remove(filename,".txt"),"_extend_long.RData"))
  #save(extend_wide,extend_arr,file=paste0(comsolpfad,str_remove(filename,".txt"),"_extend_wide.RData"))
  save(extend_arr,file=paste0(comsolpfad,str_remove(filename,".txt"),"_extend_arr.RData"))
}


# unique(sweep_long$z)
# ggplot(subset(sweep_long,z==143&injection_rate == 0.1))+
#   geom_raster(aes(DS_1,DS_2,fill=CO2_mol_per_m3))+
#   scale_fill_distiller(palette="Spectral")+
#   labs(title="sweep")
# ggplot(subset(extend_long,z==143&injection_rate == 0.1))+
#   geom_raster(aes(DS_1,DS_2,fill=CO2_mol_per_m3))+
#   scale_fill_distiller(palette="Spectral")+
#   labs(title="extend")