# prosperHorizSys.R
library(rOpenserver)
library(tidyr)
library(ggplot2)

# function to get the filename for the model
get_model_filename <- function(model) {
    models_dir <- system.file("models", package = "rOpenserver")
    model_file <- file.path(models_dir, model)
    if (!file.exists(model_file)) stop("Model not found ...") else
        return(model_file)
}
# well lengths and vertical anisotropy vectors
well_length <- seq(500, by = 500, length.out = 17)
kv_kh <- c(0.001, 0.005, 0.01, 0.1, 1.0)
df <- data.frame(well_length)
# Initialize and start OpenServer
pserver <- .OpenServer$new()       # this uses the R6 class directly
cmd = "PROSPER.START"
DoCmd(pserver, cmd)
[1] 0
# open model
model_file <- get_model_filename(model = "HORWELLDP.OUT")
open_cmd <- "PROSPER.OPENFILE"
open_cmd <- paste0(open_cmd, "('", model_file, "')")
DoCmd(pserver, open_cmd)
[1] 0
# iterate through anisotropy values
for (k in kv_kh) {
    DoSet(pserver, "PROSPER.SIN.IPR.Single.Vans", k)
    # iterate through all well length values of interest
    i <-  1
    for (wlen in df[["well_length"]]) {
        # set well length
        DoSet(pserver, "PROSPER.SIN.IPR.Single.WellLen", wlen)
        # set length in zone 1
        DoSet(pserver, "PROSPER.SIN.IPR.Single.HorizdP[0].ZONLEN", wlen)
        DoCmd(pserver, "PROSPER.anl.SYS.CALC")    # do calculation
        # store liquid rate result in dataframe for each anisotropy scenario
        df[[as.character(k)]][i] <-
            as.double(DoGet(pserver, "PROSPER.OUT.SYS.RESULTS[0][0][0].SOL.LIQRATE"))
        i <-  i + 1
    }
}
print(df)
   well_length     0.001     0.005      0.01       0.1        1
1          500  5122.723  6654.130  7376.862  9225.225 10074.80
2         1000  6502.892  8158.698  8763.273 10010.193 10500.44
3         1500  7442.298  8947.091  9422.393 10368.828 10711.88
4         2000  8081.572  9416.705  9815.624 10576.695 10838.98
5         2500  8565.691  9734.755 10077.367 10710.160 10924.37
6         3000  8909.852  9963.207 10263.006 10803.312 10984.46
7         3500  9177.887 10134.045 10400.360 10870.907 11027.95
8         4000  9391.362 10265.520 10505.026 10921.186 11059.91
9         4500  9564.346 10368.811 10585.193 10959.103 11083.47
10        5000  9706.521 10451.304 10647.387 10987.987 11100.84
11        5500  9824.676 10517.964 10696.969 11009.371 11112.55
12        6000  9923.822 10571.524 10736.874 11025.970 11121.27
13        6500 10007.648 10614.573 10768.802 11038.446 11124.79
14        7000 10079.021 10650.133 10794.942 11044.998 11127.86
15        7500 10140.119 10679.603 10813.414 11051.086 11129.07
16        8000 10192.705 10701.022 10830.166 11055.076 11128.84
17        8500 10238.145 10720.846 10843.626 11057.347 11127.42
# convert dataframe to tidy dataset
df_gather <- gather(df, kv_kh, liquid_rate, '0.001':'1')

# plot
g <- ggplot(df_gather, aes(x = well_length, y = liquid_rate, color = kv_kh)) +
    geom_line() +
    geom_point()

print(g)

# shutdown Prosper
Sys.sleep(3)
command = "PROSPER.SHUTDOWN"
pserver$DoCmd(command)
[1] 0
pserver <- NULL