#################################################################
##                          Section 1                          ##
##                     Grid initialization                     ##
#################################################################

n <- 5
m <- 5

types <- c("scratch", "phreeqc", "rds")

# initsim <- c("SELECTED_OUTPUT", "-high_precision true",
#              "-reset false",
#              "USER_PUNCH",
#              "-head C Ca Cl Mg pH pe O2g Calcite Dolomite",
#              "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
#              "SOLUTION 1",
#              "units mol/kgw",
#              "temp 25.0", "water 1",
#              "pH 9.91 charge",
#              "pe 4.0",
#              "C   1.2279E-04",
#              "Ca  1.2279E-04",
#              "Cl 1E-12",
#              "Mg 1E-12",
#              "PURE 1",
#              "O2g     -0.6788 10.0",
#              "Calcite  0.0 2.07E-3",
#              "Dolomite 0.0 0.0",
#              "END")

# needed if init type is set to "scratch" 
# prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")

init_cell <- list(
  "C" = 1.2279E-4,
  "Ca" = 1.2279E-4,
  "Cl" = 0,
  "Mg" = 0,
  "pH" = 9.91,
  "pe" = 4,
  "O2g" = 10,
  "Calcite" = 2.07e-4,
  "Dolomite" = 0
)

grid <- list(
  n_cells = c(n, m),
  s_cells = c(1,1),
  type = types[1],
  init_cell = as.data.frame(init_cell),
  props = names(init_cell),
  init_script = NULL
)


##################################################################
##                          Section 2                           ##
##         Diffusion parameters and boundary conditions         ##
##################################################################

init_diffu <- c(
  "C" = 1.2279E-4,
  "Ca" = 1.2279E-4,
  "Cl" = 0,
  "Mg" = 0,
  "pe" = 4,
  "pH" = 7
)

alpha_diffu <- c(
  "C" = 1E-4,
  "Ca" = 1E-4,
  "Cl" = 1E-4,
  "Mg" = 1E-4,
  "pe" = 1E-4,
  "pH" = 1E-4
)

vecinj_diffu <- list(
  list(
    "C" = 0,
    "Ca" = 0,
    "Cl" = 0.002,
    "Mg" = 0.001,
    "pe" = 4,
    "pH" = 9.91
  )
)

inner_index <- c(5, 15, 25)
inner_vecinj_index <- rep(1, 3)

vecinj_inner <- cbind(inner_index, inner_vecinj_index)

boundary <- list(
  "N" = rep(1, n),
  "E" = rep(0, n),
  "S" = rep(0, n),
  "W" = rep(0, n)
)

diffu_list <- names(alpha_diffu)

diffusion <- list(
  init = init_diffu,
  vecinj = do.call(rbind.data.frame, vecinj_diffu),
  #vecinj_inner = vecinj_inner,
  vecinj_index = boundary,
  alpha = alpha_diffu
)

##################################################################
##                          Section 3                           ##
##                      Phreeqc simulation                      ##
##################################################################

db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
  package = "RedModRphree"
), is.db = TRUE)

phreeqc::phrLoadDatabaseString(db)

# NOTE: This won't be needed in the future either. Could also be done in a. pqi
# file
base <- c(
  "SOLUTION 1",
  "units mol/kgw",
  "temp 25.0",
  "water 1",
  "pH 9.91 charge",
  "pe 4.0",
  "C   1.2279E-04",
  "Ca  1.2279E-04",
  "Mg 0.001",
  "Cl 0.002",
  "PURE 1",
  "O2g -0.1675 10",
  "KINETICS 1",
  "-steps 100",
  "-step_divide 100",
  "-bad_step_max 2000",
  "Calcite", "-m      0.000207",
  "-parms  0.0032",
  "Dolomite",
  "-m      0.0",
  "-parms  0.00032",
  "END"
)

selout <- c(
  "SELECTED_OUTPUT", "-high_precision true", "-reset false",
  "-time", "-soln", "-temperature true", "-water true",
  "-pH", "-pe", "-totals C Ca Cl Mg",
  "-kinetic_reactants Calcite Dolomite", "-equilibrium O2g"
)

# TODO: dt and iterations

iterations <- 10

setup <- list(
  # bound = myboundmat,
  base = base,
  first = selout,
  # initsim = initsim,
  # Cf = 1,
  grid = grid,
  diffusion = diffusion,
  prop = names(init_cell),
  immobile = c(7, 8, 9),
  kin = c(8, 9),
  ann = list(O2g = -0.1675),
  # phase = "FLUX1",
  # density = "DEN1",
  reduce = FALSE,
  # snapshots = demodir, ## directory where we will read MUFITS SUM files
  # gridfile = paste0(demodir, "/d2ascii.run.GRID.SUM")
  # init = init,
  # vecinj = vecinj,
  # cinj = c(0,1),
  # boundary = boundary,
  # injections = FALSE,
  s_grid = c(n, m),
  n_grid = c(n, m),
  dt = 1,
  iterations = iterations,
  timesteps = rep(1, iterations)
)

# not needed yet
# signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5)
# prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act")