# prosperOilwellHtc
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
htc <- seq(2, by = 0.5, to = 10)
df <- data.frame(htc)
df["liquid_rate"]     <- NA
df["wellhead_temp"]   <- NA
df["sol_bhp"]         <- NA
df["sol_dp_friction"] <- NA

# openserver variables
sys_sol_liqrate       <- "PROSPER.OUT.SYS.RESULTS[0][0][0].SOL.LIQRATE"
sys_sol_whtemperature <- "PROSPER.OUT.SYS.RESULTS[0][0][0].SOL.WHTEMPERATURE"
sys_sol_bhp           <- "PROSPER.OUT.SYS.Results[0].Sol.BHP"
sys_prfric            <- "PROSPER.OUT.SYS.Results[0].Sol.PRFRIC"
# Initialize and start OpenServer
pserver <- OpenServerR6$new()     # another way to start an OpenServer instance
cmd = "PROSPER.START"
DoCmd(pserver, cmd)
[1] 0
# open model
model_file <- get_model_filename(model = "oilwell.Out")
open_cmd <- "PROSPER.OPENFILE"
open_cmd <- paste0(open_cmd, "('", model_file, "')")
DoCmd(pserver, open_cmd)
[1] 0
# iterate through heat transfer coefficient vector
i <-  1
for (h in df[["htc"]]) {
    # set HTC
    DoSet(pserver, "PROSPER.SIN.EQP.Geo.Htc", h)
    DoCmd(pserver, "PROSPER.anl.SYS.CALC")    # do calculation
    # store liquid rate, wellhead temperature, solution BHP and dP friction
    df[["liquid_rate"]][i]     <- as.double(DoGet(pserver, sys_sol_liqrate))
    df[["wellhead_temp"]][i]   <- as.double(DoGet(pserver, sys_sol_whtemperature))
    df[["sol_bhp"]][i]         <- as.double(DoGet(pserver, sys_sol_bhp))
    df[["sol_dp_friction"]][i] <- as.double(DoGet(pserver, sys_prfric))
    i <-  i + 1
}
print(df)
    htc liquid_rate wellhead_temp  sol_bhp sol_dp_friction
1   2.0    12207.70      193.8414 2669.929        514.1321
2   2.5    12201.75      190.1239 2671.276        512.6521
3   3.0    12195.78      186.5259 2672.628        511.1868
4   3.5    12189.82      183.0434 2673.978        509.7431
5   4.0    12182.66      179.6687 2675.599        508.3648
6   4.5    12176.46      176.4035 2677.004        506.7877
7   5.0    12170.23      173.2415 2678.415        505.2016
8   5.5    12164.05      170.1792 2679.813        503.6503
9   6.0    12157.92      167.2131 2681.202        502.1317
10  6.5    12151.84      164.3399 2682.580        500.6453
11  7.0    12145.80      161.5562 2683.948        499.1901
12  7.5    12139.76      158.8588 2685.316        497.7242
13  8.0    12133.74      156.2446 2686.678        496.2813
14  8.5    12127.75      153.7108 2688.035        494.8627
15  9.0    12120.55      151.2492 2689.666        493.6298
16  9.5    12102.93      148.8171 2693.657        493.7540
17 10.0    12097.17      146.5076 2694.961        492.4105
# convert dataframe to tidy dataset
df_gather <- gather(df, var, value, liquid_rate:sol_dp_friction)

# plot
g <- ggplot(df_gather, aes(x = htc, y = value)) +
    geom_line() +
    geom_point() +
    facet_wrap(var ~., scales = "free_y")

print(g)

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