R/package-gmg.R

Defines functions rmf_write_gmg rmf_read_gmg rmf_create_gmg

Documented in rmf_create_gmg rmf_read_gmg rmf_write_gmg

#' Create an \code{RMODFLOW} gmg object
#' 
#' \code{rmf_create_gmg} creates an \code{RMODFLOW} gmg object.
#' 
#' @param rclose residual convergence criterion for the inner iteration; defaults to 0.001
#' @param iiter maximum number of PCG iterations for each linear solution; defaults to 100
#' @param hclose head-change convergence criterion for nonlinear problems; ignored of mxiter = 1. Defaults to 0.001
#' @param mxiter maximum number of outer iterations; defaults to 50
#' @param damp value of the damping parameter; defaults to 0.5
#' @param iadamp flag controlling adaptive damping. Either 0 (no damping), 1 (Cooley's method for first nonlinear iteration; default) or 2 (relative reduced residual damping)
#' @param ioutgmg flag controlling solver output. Either 0 (only solver input is printed), 1 (detailed solver output; default), 2 (basic solver output), 3 (same as 1 but printed to terminal) or 4 (same as 2 but printed to terminal)
#' @param iunitmhc flag and unit number for writing the maximum head change values. If > 0, these values are written to unit number iunitmhc which should be listed in the Name file. If <= 0 (default), no values are written.
#' @param ism type of smoother in preconditioner. Either 0 (ILU(0) smoothing; default) or 1 (Symmetric Gauss-Seidel smoothing)
#' @param isc flag controlling semi-coarsening in preconditioner. Either 0 (rows, columns & layers; default), 1 (rows & columns), 2 (columns & layers), 3 (rows & layers) or 4 (no coarsening)
#' @param dup maximum damping value that should be applied to any iteration when the solver is not oscillating. Only used when iadamp = 2; defaults to 0.5
#' @param dlowminimum damping value to be generated by the adaptive-damping procedure. Only used when iadamp = 2; defaults to 0.01
#' @param chglimit maximum allowed head change at any cell between outer iterations. Only used when iadamp = 2; defaults to 0.001
#' @param relax relaxation parameter for the ILU PCG method. Defaults to 1. Only used when isc = 4.
#'
#' @return Object of class gmg
#' @export
#' @seealso \code{\link{rmf_read_gmg}}, \code{\link{rmf_write_gmg}} and \url{https://water.usgs.gov/ogw/modflow/MODFLOW-2005-Guide/index.html?gmg.htm}
rmf_create_gmg <- function(rclose = 0.001,
                           iiter = 100,
                           hclose = 0.001,
                           mxiter = 50,
                           damp = 0.5,
                           iadamp = 1,
                           ioutgmg = 1,
                           iunitmhc = 0,
                           ism = 0,
                           isc = 0,
                           dup = 0.5,
                           dlow = 0.01,
                           chglimit = 0.001,
                           relax = 1) {
  gmg <- NULL
  
  # data set 0
  # to provide comments, use ?comment on the resulting gmg object
  
  # data set 1
  gmg$rclose <- rclose
  gmg$iiter <- iiter
  gmg$hclose <- hclose
  gmg$mxiter <- mxiter
  
  # data set 2
  gmg$damp <- damp
  gmg$iadamp <- iadamp
  gmg$ioutgmg <- ioutgmg
  gmg$iunitmhc <- ifelse(iunitmhc < 0, 0, iunitmhc)
  
  # data set 3
  gmg$ism <- ism
  gmg$isc <- isc
  if(gmg$iadamp == 2) {
    gmg$dup <- dup
    gmg$dlow <- dlow
    gmg$chglimit <- chglimit
  }
  
  # data set 4
  if(gmg$isc == 4) gmg$relax <- relax
  
  class(gmg) <- c('gmg','rmf_package')
  return(gmg)
}

#' Read a MODFLOW Geometric Multigrid Solver file
#' 
#' \code{rmf_read_gmg} reads in a MODFLOW Geometric Multigrid Solver file and returns it as an \code{\link{RMODFLOW}} gmg object.
#' @param file filename; typically '*.gmg'
#' @param ... ignored
#' @return object of class gmg
#' @export
#' @seealso \code{\link{rmf_write_gmg}}, \code{\link{rmf_create_gmg}} and \url{https://water.usgs.gov/ogw/modflow/MODFLOW-2005-Guide/index.html?gmg.htm}
rmf_read_gmg <- function(file = {cat('Please select gmg file ...\n'); file.choose()}, ...) {
  
  gmg_lines <- readr::read_lines(file, lazy = FALSE)
  gmg <- list()
  
  # data set 0
  data_set_0 <- rmfi_parse_comments(gmg_lines)
  comment(gmg) <- data_set_0$comments
  gmg_lines <- data_set_0$remaining_lines
  rm(data_set_0)
  
  # data set 1
  data_set_1 <- rmfi_parse_variables(gmg_lines)
  gmg$rclose <- as.numeric(data_set_1$variables[1])
  gmg$iiter <- as.numeric(data_set_1$variables[2])
  gmg$hclose <- as.numeric(data_set_1$variables[3])
  gmg$mxiter <- as.numeric(data_set_1$variables[4])
  gmg_lines <- data_set_1$remaining_lines
  rm(data_set_1)
  
  # data set 2
  data_set_2 <- rmfi_parse_variables(gmg_lines)
  gmg$damp <- as.numeric(data_set_2$variables[1])
  gmg$iadamp <- as.numeric(data_set_2$variables[2])
  gmg$ioutgmg <- as.numeric(data_set_2$variables[3])
  gmg$iunitmhc <- ifelse(is.na(as.numeric(data_set_2$variables[4])), 0, as.numeric(data_set_2$variables[4]))
  if(gmg$iunitmhc < 0) gmg$iunitmhc <- 0
  gmg_lines <- data_set_2$remaining_lines
  rm(data_set_2)
  
  # data set 3
  data_set_3 <- rmfi_parse_variables(gmg_lines)
  gmg$ism <- as.numeric(data_set_3$variables[1])
  gmg$isc <- as.numeric(data_set_3$variables[2])
  if(gmg$iadamp == 2) {
    gmg$dup <- as.numeric(data_set_3$variables[3])
    gmg$dlow <- as.numeric(data_set_3$variables[4])
    gmg$chglimit <- as.numeric(data_set_3$variables[5])
  }
  gmg_lines <- data_set_3$remaining_lines
  rm(data_set_3)
  
  # data set 4
  if(gmg$isc == 4) {
    data_set_4 <- rmfi_parse_variables(gmg_lines)
    gmg$relax <- as.numeric(data_set_4$variables[1])
    gmg_lines <- data_set_4$remaining_lines
    rm(data_set_4)
  }
  
  class(gmg) <- c('gmg','rmf_package')
  return(gmg)
}

#' Write a MODFLOW Geometric Multigrid Solver file
#' 
#' @param gmg an \code{\link{RMODFLOW}} gmg object
#' @param file filename to write to; typically '*.gmg'
#' @param ... ignored
#' @return \code{NULL}
#' @export
#' @seealso \code{\link{rmf_read_gmg}}, \code{\link{rmf_create_gmg}} and \url{https://water.usgs.gov/ogw/modflow/MODFLOW-2005-Guide/index.html?gmg.htm}
rmf_write_gmg <- function(gmg,
                          file = {cat('Please select gmg file to overwrite or provide new filename ...\n'); file.choose()}, ...) {
  
  # data set 0
  v <- packageDescription("RMODFLOW")$Version
  cat(paste('# MODFLOW Geometric Multigrid Solver file created by RMODFLOW, version',v,'\n'), file=file)
  cat(paste('#', comment(gmg)), sep='\n', file=file, append=TRUE)
  
  # data set 1
  rmfi_write_variables(gmg$rclose, as.integer(gmg$iiter), gmg$hclose, as.integer(gmg$mxiter), file=file)
  
  # data set 2
  rmfi_write_variables(gmg$damp, as.integer(gmg$iadamp), as.integer(gmg$ioutgmg), as.integer(gmg$iunitmhc), file = file)
  
  # data set 3
  rmfi_write_variables(as.integer(gmg$ism), as.integer(gmg$isc), rmfi_ifelse0(gmg$iadamp == 2, c(gmg$dup, gmg$dflow, gmg$chglimit), ''), file = file)
  
  # data set 4
  if(gmg$isc == 4) {
    rmfi_write_variables(gmg$relax, file=file)
  }
  
}
rogiersbart/RMODFLOW documentation built on Jan. 14, 2023, 4:21 a.m.