r2ogs5.Rmd

---
title: ""
author: Philipp Schad
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simple Example for r2ogs5}
  %\VignetteKeyword{HTML}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(tibble.print_min = 4L, tibble.print_max = 4L)
options(rmarkdown.html_vignette.check_title = FALSE)
```

# Basic functions of r2ogs5

```{r setup, message=FALSE}
library(r2ogs5)
```

In this vignette, the basic functions of *r2ogs5* are illustrated with a simple 
conservative mass transport simulation. 

1. [Installation](#inst)
2. [Setting up a simulation manually](#sim1)
3. [Modifying a simulation object, and writing input files](#mod) 
4. [Setting up a simulation from input files](#sim2)
5. [Run simulations, and retrieve output](#run)

## Installation <a name="inst"></a>

The first option for installation is directly from the online repository via
```{r, eval=FALSE}
devtools::install_git("https://gitlab.opengeosys.org/ogs5/r2ogs5")
```
or 
```{r, eval=FALSE}
remotes::install_git("https://gitlab.opengeosys.org/ogs5/r2ogs5")
```

in your R console.  

The second option is first cloning the repository by typing in a terminal 
```{bash, eval=FALSE}
git clone https://gitlab.opengeosys.org/ogs5/r2ogs5.git
cd r2ogs5
R --slave -e 'devtools::install(".")'
```

Note that it is necessary to have git https://www.git-scm.com/
and the R-package `devtools` installed on your system. 
  
Finally, the actual simulation program, OpenGeoSys 5 needs to be available on the system.
*r2ogs5* does not come with OGS5 executables, as they usually
need to be compiled specifically for the machine run on. For this tutorial, the default compiled
executables should be sufficient which can be downloaded from https://www.opengeosys.org/releases/ for Linux, macOS, and Windows.  
An OGS5 executable needs to be available. A function to browse the OS "PATH" variable 
is available, to search for an OGS5 executable by name. Alternatively, the name (here ogs_fem), and location of the OGS5 executable have to be set accordingly via `options`. 

```{r}
library(r2ogs5)
```

```{r, eval=TRUE, echo=FALSE}
options(r2ogs5.default_ogs5_bin = search_ogs5_bin_path(ogs_exe = "ogs_fem", return_ogs5_bin = TRUE))
```

```{r, eval=FALSE}
# Note: the next line only makes sense if an executable "ogs_fem" exists
# in the directory "/usr/local/bin/"
options(r2ogs5.default_ogs5_bin = "/usr/local/bin/ogs_fem")
```

The current location of the executable can be shown by 
```{r}
getOption("r2ogs5.default_ogs5_bin")
```


There are more options to set a default simulation directory, benchmark directory, and python environment (needed for some simulation output gathering functionalities). By preference, all options for a  *r2ogs5* project can be stored in a "config.yml" file in the corresponding working directory and updated via `update_options_from_config` (see `?update_options_from_config`). An example for a config.yml file is: 

```{bash, eval = FALSE}
default: 
  r2ogs5.default_ogs5_bin: "/usr/local/bin/ogs_fem"
  r2ogs5.default_ogs5_bin_path: "/usr/local/bin/"

reactive_transport: 
  r2ogs5.default_ogs5_bin: "/usr/local/bin/ogs_ipqc"
  r2ogs5.default_ogs5_bin_path: "/usr/local/bin/"

```
Here, a default setting and an alternative is available. To use the ogs executable 
compiled for reactive transport simulations "ogs_ipqc", the options could be loaded via 
```{r, eval = FALSE}
update_options_from_config(config = "reactive_transport")
```

**Note:** In order to install the package with vignettes, e.g. by calling 

```{r, eval=FALSE}
devtools::install_git("https://gitlab.opengeosys.org/ogs5/r2ogs5", build_vignettes = TRUE)
```
a basic ogs5 executable named "ogs_fem" must be available in the search path. For Linux, the executable should be in one of the directories found when entering `Sys.getenv("PATH")` in the R console. For Windows, it should be located in *"C:/Programme/OpenGeoSys/bin/"* or *"C:/Programs/OpenGeoSys/bin/"* under the current implementation.  


## Setting up a simulation <a name="sim"></a>

By means of a very simple hydraulic transport flow model, the use of the central 
package functions to create input, run a simulation, and retrieve output are demonstrated. 

A relative directory "tmp" to write OGS5 files in will be defined here, it can be replaced by the user with any path to a custom directory of choice. 
```{r}
temp_dir <- "tmp"
```


First, a simulation object of 
class *ogs5* has to be created, where these blocs are added. This is done using
the assign operator ( `<-` ), and the function `create_ogs5()`. 
```{r}
ex1 <- create_ogs5(sim_name = "ex1",
                   sim_id = 1L,
                   sim_path = temp_dir)
```

The created *ogs5* object contains two sections; one named `output`, where results can be stored after simulation, and another called `input`, to store the different input objects of specific classes that corresponding to OpenGeoSys5 input files. The different input classes contain one or more blocs of a specific class that define a simulation parameter.  

The functions to set up a simulation with *r2ogs5* follow the pattern 
`input_add_<file extension>bloc()`, and can be used repeatedly to add new input blocs to the simulation. First, a process-bloc, which is the
central object for any OGS5 simulation, is defined for a liquid flow process.
```{r}
ex1 <- input_add_pcs_bloc(x=ex1, pcs_name = "waterflow", PCS_TYPE = "LIQUID_FLOW",
                          PRIMARY_VARIABLE = "PRESSURE1")
```
The overall structure of *ogs5* objects are lists, therefore blocs can be extracted via the list operator `$`.

```{r}
str(ex1)
ex1$input$pcs$waterflow
```
Depending on the bloc class, and the process that is simulated, there are different arguments that need to be
specified in the functions to create input blocs. In general, all uppercase arguments are
equivalent to **ogs5** keywords, and can be looked up at https://ogs5-keywords.netlify.app/ogs/wiki/public/doc-auto. A short overview of the input classes needed for this example, and their function is given in \Cref{tab:inputClass}.

In the following, another `pcs` bloc for a mass transport process will be added to `ex1`, and other blocs in 
the same fashion. 
```{r}
ex1 <- input_add_pcs_bloc(x=ex1, pcs_name = "tracer",
                          PCS_TYPE = "MASS_TRANSPORT",
                          PRIMARY_VARIABLE = "Tracer1",
                          RELOAD = c("1", "1"))

ex1 <- input_add_bc_bloc(x=ex1, bc_name = "pressure_out",
                         PRIMARY_VARIABLE = "PRESSURE1",
                         DIS_TYPE = "CONSTANT 9810",
                         PCS_TYPE = "LIQUID_FLOW",
                         GEO_TYPE = "POINT point1")

ex1 <- input_add_bc_bloc(x=ex1, bc_name = "tracer_in",
                         PRIMARY_VARIABLE = "Tracer",
                         DIS_TYPE = "CONSTANT 1.19438",
                         PCS_TYPE = "MASS_TRANSPORT",
                         GEO_TYPE = "POINT point0",
                         TIM_TYPE = "CURVE 1")

ex1 <- input_add_ic_bloc(x = ex1, ic_name = "pressure",
                         PCS_TYPE = "LIQUID_FLOW",
                         PRIMARY_VARIABLE = "PRESSURE1",
                         GEO_TYPE = "DOMAIN",
                         DIS_TYPE = "CONSTANT 9810")

ex1 <- input_add_ic_bloc(x = ex1, ic_name = "tracer",
                         PCS_TYPE = "MASS_TRANSPORT",
                         PRIMARY_VARIABLE = "Tracer",
                         GEO_TYPE = "DOMAIN",
                         DIS_TYPE = "CONSTANT 0.0")

ex1 <- input_add_mcp_bloc(x = ex1,
                          NAME = "Tracer",
                          MOBILE = 1)

ex1 <- input_add_mfp_bloc(x = ex1, FLUID_NAME = "water",
                          FLUID_TYPE = "LIQUID",
                          DENSITY = "1 1.0E+3",
                          VISCOSITY = "1 1.0E-3")

ex1 <- input_add_mmp_bloc(x = ex1, GEOMETRY_DIMENSION = "1",
                          GEOMETRY_AREA = "1.2",
                          POROSITY = "1 0.38",
                          PERMEABILITY_TENSOR = "ISOTROPIC 1.0E-8",
                          MASS_DISPERSION = "1 1.0",
                          TORTUOSITY = "1 1.0",
                          NAME = "gravel")

ex1 <- input_add_msp_bloc(x = ex1, NAME = "gravel" ,
                          DENSITY = "1 1.8E+3")

ex1 <- input_add_num_bloc(x = ex1, num_name = "waterflow",
                          PCS_TYPE = "LIQUID_FLOW",
                          LINEAR_SOLVER = "2 6 1.0E-14 100 1.0 100 4")

ex1 <- input_add_num_bloc(x = ex1, num_name = "tracer",
                          PCS_TYPE = "MASS_TRANSPORT",
                          LINEAR_SOLVER = "2 6 1.0E-14 100 1.0 100 4")

ex1 <- input_add_out_bloc(x = ex1, out_name = "waterflow",
                          PCS_TYPE = "LIQUID_FLOW",
                          NOD_VALUES = "PRESSURE1\n VELOCITY_X1",
                          GEO_TYPE = "DOMAIN",
                          DAT_TYPE = "PVD",
                          TIM_TYPE = "STEPS 10")

ex1 <- input_add_out_bloc(x = ex1, out_name = "tracer",
                          PCS_TYPE = "MASS_TRANSPORT",
                          NOD_VALUES = "Tracer",
                          GEO_TYPE = "DOMAIN",
                          DAT_TYPE = "PVD",
                          TIM_TYPE = "STEPS 10")

ex1 <- input_add_out_bloc(x = ex1, out_name = "tracer_tec",
                          PCS_TYPE = "MASS_TRANSPORT",
                          NOD_VALUES = "Tracer",
                          GEO_TYPE = "DOMAIN",
                          DAT_TYPE = "TECPLOT", TIM_TYPE = "STEPS 1")

ex1 <- input_add_st_bloc(x = ex1, st_name = "waterflow",
                         PCS_TYPE = "LIQUID_FLOW",
                         PRIMARY_VARIABLE = "PRESSURE1",
                         GEO_TYPE = "POINT point0",
                         DIS_TYPE = "CONSTANT_NEUMANN 7.8611E-06")

ex1 <- input_add_tim_bloc(x = ex1,  tim_name = "waterflow",
                          PCS_TYPE = "LIQUID_FLOW",
                          TIME_START = "0",
                          TIME_END = "36000000",
                          TIME_STEPS = "229 3600")

ex1 <- input_add_tim_bloc(x = ex1,  tim_name = "tracer",
                          PCS_TYPE = "MASS_TRANSPORT",
                          TIME_START = "0",
                          TIME_END = "36000000",
                          TIME_STEPS = "229 3600")

```

\newpage
```{r inputClass, echo=FALSE, results='markup'}

df <- tibble::tibble("Input class" = c(
  "pcs",
  "gli",
  "msh",
  "num",
  "tim",
  "ic",
  "bc",
  "st",
  "mcp",
  "mfp",
  "mmp",
  "msp",
  "out"),
  Desctiption = c(
    "Process definition", 
    "Geometry",
    "Finite element mesh",
    "Numerical properties",
    "Time discretization",
    "Initial condition",
    "Boundary condition",
    "Source/Sink term",
    "Component properties",
    "Fluid properties",
    "Medium properties",
    "Solid properties",
    "Output configuration")
  )

knitr::kable(df, caption = "Most important ogs5 input classes")
```

Next, an `rfd` bloc with a time series of tracer concentrations at one point 
will be generated. 

```{r}
ex1 <- input_add_rfd_bloc(x = ex1, rfd_name = "tracer",
                          mkey = "CURVES",
                          data = tibble::tibble(
                            time = c(0, 3600, 3600.1, 720, 36000000),
                            conc = c(1,1,0,0,0)))
```

Finally, two geometry blocs are created, and added to the `ex1`. First, a `gli` 
bloc with three-dimensional coordinates $x, y, z$ for two points is added. Then a `msh` bloc is created 
via two helper functions. Here, only a one-dimensional mesh (on $x$) is created. In order to create a $3d$ mesh, for example, the arguments `ly`, `lz`, `ny`, and `nz` have to be specified accordingly. For further options please refer to the function documentation. 

```{r}
## add gli
ex1 <- input_add_gli_points(x = ex1,
                            ogs5_points = tibble::tibble(
                              x = c(0, 4.7),
                              y = c(0,0),
                              z = c(0,0),
                              name = c("point0", "point1")))
# ad msh bloc
mesh_lst <- create_structured_mesh_nodes_ele(lx = 4.7, nx = 94)
ex1 <- input_add_msh_bloc(x = ex1, msh_name = "base_mesh",
                          NODES = mesh_lst$NODES,
                          ELEMENTS = mesh_lst$ELEMENTS) 
```
The names of the points link to the boundary condition blocs (`bc`), which have been added before:
```{r}
ex1$input$bc$tracer_in$GEO_TYPE
ex1$input$bc$pressure_out$GEO_TYPE
```

The location of the points, lines, or surfaces on the finite element mesh, for which boundary conditions, sink/source terms, etc are defined in the `gli` section.
```{r}
ex1$input$gli
```


## Modifying a simulation object, and writing input files <a name="mod"></a> 

In order to overwrite an existing input bloc, the entry has to be deleted first.
Then a new bloc with the same name as before but different values can be added. 
```{r}
ex1$input$msp <- NULL
ex1 <- input_add_msp_bloc(ex1,
                          NAME = "gravel",
                          DENSITY = "1 1.79e+3")
```
Alternatively, an entry can also be changed directly by assigning a new character 
string. Here, care has to be taken to include all necessary specifications, and
follow the format requirements such as whitespaces, and line breaks. 
```{r}
ex1$input$msp$gravel$DENSITY <- "1 1.8E+3"
str(ex1$input$msp)
```

The helper function `cal_change_parameters` to use for calibration can also be used for this purpose. For more see `vignette("cal_bayesOpt"). 

All blocs can be written out by setting `type = "all"`, (default)
or only a selection of bloc types. Note that r2ogs5 will always create directories
where none exist, and overwrite existing directories. 
```{r,eval=TRUE}
ogs5_write_inputfiles(ex1, type = "all")

ogs5_write_inputfiles(ex1, type = "bc")

ogs5_write_inputfiles(ex1, type = c("bc", "mmp"))
```
If no `folderpath` is specified, the files will be written into the folder specified
at the initialization of `ex1` with `sim_path`. 

## Setting up a simulation from input files <a name="sim2"></a>

The function `input_add_blocs_from_file()` can be used to create *r2ogs5* 
simulation objects directly from the input files in an automatized fashion. As
before, an *ogs5*-object needs to be initialized first.

For example, all files from the directory *temp/* can be read in, and *r2ogs* will print the reading progress.  
```{r, eval=TRUE, results="hide", echo=TRUE}
ex2 <- create_ogs5(
  sim_name = "ex2",
  sim_id = 1L,
  sim_path = paste0(temp_dir, "/two")
  )

ex2 <- input_add_blocs_from_file(
  ogs5_obj = ex2,
  sim_basename = "ex1",
  filename = "all",
  file_dir = temp_dir
  )
```
The argument `sim_basename` should correspond to the common file-basename in the directory. The function will automatically assign the blocs with the corresponding classes, and place them into the correct location in `ex2`.  

Except for \*.tec, \*.vtk, \*.pvd all kinds of files that are in the directory are read in. Files that are not among the basic *OGS5* files are added into `input$additional`. 
It is important to have 'clean' input files in the directory, best without any comments (e.g. \\\\). Furthermore, a single input file or a list of inputfiles can be specified explicitly, to avoid reading any unnecessary files from the directory.
```{r, results='hide'}
ex2 <- input_add_blocs_from_file(
  ogs5_obj = ex2,
  sim_basename = "ex1",
  filename = list("ex1.bc", "ex1.mmp"),
  file_dir = temp_dir
  )
```

## Run simulations and retrieve output <a name="run"></a>  

If the simulation object is ready, and written to a directory, the simulation can finally be handed over to *OGS5*, as shown below with all arguments set to default. 
```{r, eval=TRUE}
ogs5_run(ogs5_obj = ex1,
         ogs_exe = NULL, # default path
         run_path = attributes(ex1)$sim_path,
         log_output = TRUE,
         log_path = paste0(attributes(ex1)$sim_path, "/log"), 
         use_mpi = FALSE,
         wait = TRUE)
```
In order to run, an appropriate *OGS5* executable should be provided, as explained at the beginning. Depending on the type of simulation, special compiling options need to be set. For a download and compilation guide of *OGS5* please refer to https://github.com/ufz/ogs5. 
The argument `log_output` allows writing the *OGS5* output into a file specified under `log_path` or display the output in the console by setting `log_output` to `FALSE`. 
For Parallelized simulation, `use_mpi` must be set to true, the `number_of_cores` specified, and an appropriate *OGS5* executable provided, of course.


The output files specified in `input$out` are written into the same directory where
the input files lie. 
```{r}
str(ex1$input$out)
```
Here we have two *"PVD"*, and one "*TECPLOT*" as output files. The first stores every
10th time step of the variable `PRESSURE`, the second every 10th time step of `Tracer`, and the last 
every time step of `Tracer` over the whole geometrical domain of the model. 

The "tecplot" output can be added to the `ex1$output` via
```{r}
ex1$output <- ogs5_read_many_tecplots(filepath = attributes(ex1)$sim_path,
                  geo_object = "domain")
ex1$output %>% dplyr::select(!filepath)
```
which yields a regular *tibble* with the observed variable in time, and space.
Depending on their dimensionality, results can then be easily plotted.

```{r, message=FALSE, warning=FALSE, results='hide', fig.width=5, fig.asp=0.6}
library(dplyr)
library(ggplot2)
breakthrough <- ex1$output %>%
  select(TIME, X, Tracer) %>% 
  filter(X == 4.7)

ggplot(breakthrough, aes(x = TIME, y = Tracer)) +
  geom_point() +
  labs(title = "End point") +
  theme_classic()
```
```{r, message=FALSE, warning=FALSE, results='hide', fig.width=7, fig.asp=0.6}
ggplot(ex1$output, aes(x = X, y = TIME, fill = Tracer)) +
  geom_tile() +
  labs(title = "Model domain") +
  scale_fill_viridis_c(trans = "log") +
  theme_classic()
  
```