diff --git a/packrat/packrat.lock b/packrat/packrat.lock
index 06abd433d2ba9a704607383ffcc1b135de24601b..b5743cdb85e927893ec150ca789d9377c73d2b3b 100644
--- a/packrat/packrat.lock
+++ b/packrat/packrat.lock
@@ -1,7 +1,7 @@
PackratFormat: 1.4
PackratVersion: 0.5.0
RVersion: 4.0.3
-Repos: CRAN=https://cloud.r-project.org
+Repos: CRAN=https://cran.rstudio.com/
Package: BH
Source: CRAN
@@ -146,6 +146,12 @@ Source: CRAN
Version: 0.6.27
Hash: 8ba156629c1537635b39332c93d8ed33
+Package: doParallel
+Source: CRAN
+Version: 1.0.16
+Hash: 0bfafffe940c8fc4bc4f785018b17dde
+Requires: foreach, iterators
+
Package: dplyr
Source: CRAN
Version: 1.0.2
@@ -169,6 +175,12 @@ Source: CRAN
Version: 0.4.1
Hash: fc0a252b8e427847d13e89f56ab4665e
+Package: foreach
+Source: CRAN
+Version: 1.5.1
+Hash: 6dda99e5cba20abf021b0e6cfb3472df
+Requires: iterators
+
Package: fs
Source: CRAN
Version: 1.5.0
@@ -234,6 +246,11 @@ Source: CRAN
Version: 0.3.1
Hash: 9d6de5178c1cedabfb24e7d2acc9a092
+Package: iterators
+Source: CRAN
+Version: 1.0.13
+Hash: 750adbdf6750e2dbcf2066c443d64892
+
Package: jsonlite
Source: CRAN
Version: 1.7.1
@@ -363,6 +380,11 @@ Version: 1.1
Hash: 080900612c4d677d521131810ce29f26
Requires: Rcpp, askpass, curl
+Package: rappdirs
+Source: CRAN
+Version: 0.3.3
+Hash: fd602edfd7dbf365642e35dd9cc97958
+
Package: rcmdcheck
Source: CRAN
Version: 1.3.3
@@ -388,6 +410,12 @@ Source: CRAN
Version: 2.2.0
Hash: 7f3378b51f10acd43af6814111572889
+Package: reticulate
+Source: CRAN
+Version: 1.18
+Hash: 4837a6417a9d3e7f6f934847e535253c
+Requires: Rcpp, jsonlite, rappdirs
+
Package: rex
Source: CRAN
Version: 1.2.0
diff --git a/vignettes/theis_vignette.Rmd b/vignettes/theis_vignette.Rmd
index 72d7c24ec8c29084246f1bbd654a79966f879186..d5a6e4de03afbcacd784c462ea7bfc90e755f3f5 100644
--- a/vignettes/theis_vignette.Rmd
+++ b/vignettes/theis_vignette.Rmd
@@ -16,6 +16,7 @@ knitr::opts_chunk$set(
```{r setup}
library(r2ogs6)
+library(tidyverse)
```
## Hi there!
@@ -82,23 +83,14 @@ get_values <- function(value,
# demo
paste(get_values(100), collapse = " ")
```
-Let's use this function to set up our parameters.
+We can now use this function for a minimal working example. Let's set up an ensemble with only the `storage` parameter. First we create the corresponding vector...
```{r}
-
-permeability <- get_values(ogs6_obj$media[[1]]$properties[[1]]$value)
-
-porosity <- get_values(ogs6_obj$media[[1]]$properties[[3]]$value,
- zero_to_one = TRUE)
-
+# Create the vector based on the original value
storage <- get_values(ogs6_obj$media[[1]]$properties[[4]]$value)
-
-# check values
-paste(permeability, collapse = " ")
-paste(porosity, collapse = " ")
paste(storage, collapse = " ")
```
-After we've prepared our parameters, we've got everything we need for a minimal working example. Let's set up an ensemble with only the `storage` vector first.
+... then we pass the original reference as well as the vector to an ensemble object.
```{r}
# First we define an ensemble object
@@ -114,9 +106,9 @@ ogs6_ens$run_simulation()
### Plot results
-#### Single plot
+#### Single simulation
-Our simulations (hopefully) ran successfully - great! Now it'd be nice to see some results. For now, we're interested in the `pressure` and `v`elocity data.
+Our simulations (hopefully) ran successfully - great! Now it'd be nice to see some results. For now, we're interested in the `pressure` and `v`elocity data of the first simulation.
```{r}
# This will get a list of tibbles - one for each simulation.
@@ -129,20 +121,163 @@ storage_tbls <-
You may leave out the `point_ids` parameter. It defaults to all points in the dataset - I specify it here because there's about 500 which would slow down building this vignette significantly.
+```{r}
+# Let's look at the first 10 rows of the dataset for the first simulation.
+head(storage_tbls[[1]], 10)
+```
+
+You can see there's one row per timestep for each point ID. Say we want to plot the *average* pressure for each `timestep`, that is, the mean of the `pressure` values where the rows have the same `timestep`.
+
+
+```{r}
+# Get average pressure from first dataframe
+avg_pr_df <- storage_tbls[[1]] %>%
+ group_by(timestep) %>%
+ summarise(avg_pressure = mean(pressure))
+avg_pr_df
+```
+
+
+```{r}
+# Plot pressure over time for first simulation
+ggplot(data = avg_pr_df) +
+ geom_point(mapping = aes(x = as.numeric(row.names(avg_pr_df)),
+ y = avg_pressure)) + xlab("Timestep")
+
+
+```
+#### Multiple simulations
+
+```{r}
+# Get list of dataframes that each look like avg_pr_df
+avg_pr_dfs <- list()
+
+for(i in seq_len(length(storage_tbls))){
+ old_df <- storage_tbls[[i]]
+ new_df <- old_df %>%
+ group_by(timestep) %>%
+ summarise(avg_pressure = mean(pressure))
+
+ # Add column so we can separate the values based on their old df
+ new_df$df_id <- rep(i, nrow(new_df))
+
+ # Add column so we can plot based on old df rownames
+ new_df$orig_rwns <- as.numeric(row.names(new_df))
+ avg_pr_dfs <- c(avg_pr_dfs, list(new_df))
+}
+
+# Combine dataframes into single dataframe for plot
+avg_pr_combined_df <- bind_rows(avg_pr_dfs)
+
+# Make plot
+ggplot(avg_pr_combined_df,
+ aes(x = orig_rwns,
+ y = avg_pressure)) +
+ geom_point(aes(color = df_id)) +
+ geom_line(aes(color = df_id)) +
+ xlab("Timestep") +
+ facet_grid(rows = vars(df_id))
+
+```
+
+Now we have the average `pressure` over time for each of the simulations (rows). Since they're pretty similar, let's put them in the same plot for a better comparison.
+
+```{r}
+ggplot(avg_pr_combined_df, aes(x = orig_rwns,
+ y = avg_pressure,
+ group = df_id)) +
+ geom_point(aes(color = df_id)) +
+ geom_line(aes(color = df_id)) +
+ xlab("Timestep")
+```
+
+
+#### Matrix plot
+
+So far we've only considered the `storage` parameter. Now we want to have a look at how the other two parameters influence the simulation and each other. Our `storage` vector had 7 values. Assuming the same length for `permeability` and `porosity` would lead us to a whopping 343 combinations, so we will use our `get_values` function to create smaller vectors. Say we choose a length of 4 with deviations of -1%, -10% and -50% as well as the original value.
+
+```{r}
+
+new_deviations <- c(-0.5, -0.1, -0.01, 0)
+
+permeability <- get_values(ogs6_obj$media[[1]]$properties[[1]]$value,
+ new_deviations)
+porosity <- get_values(ogs6_obj$media[[1]]$properties[[3]]$value,
+ new_deviations)
+storage <- get_values(ogs6_obj$media[[1]]$properties[[4]]$value,
+ new_deviations)
+
+# Total length of each vector we want to pass to the ensemble
+vec_len <- length(permeability) * length(porosity) * length(storage)
+
+# Replicate elements based on vec_len
+per_long <- rep(permeability, (vec_len / length(permeability)))
+por_long <- rep(porosity, (vec_len / length(porosity)))
+sto_long <- rep(storage, (vec_len / length(storage)))
+```
+
+Now we take these long vectors and feed them to a new ensemble. We can base this on the same `OGS6` object as before since our project file hasn't changed.
+
+```{r}
+ogs6_ens_big <-
+ OGS6_Ensemble$new(
+ ogs6_obj = ogs6_obj,
+ parameters = list(list(ogs6_obj$media[[1]]$properties[[1]]$value, per_long),
+ list(ogs6_obj$media[[1]]$properties[[3]]$value, por_long),
+ list(ogs6_obj$media[[1]]$properties[[4]]$value, sto_long)),
+ use_original_values = FALSE)
+
+#ogs6_ens_big$run_simulation()
+```
+
+This will take a while on most systems since we're creating 64 objects and running the same number of simulations. Go grab a coffee (inclusive-)or two. As soon as the simulations are done, we can extract the point data just as we did before.
+
```{r}
# This will get a list of tibbles - one for each simulation.
-storage_tbls <-
+per_por_sto_tbls <-
lapply(ogs6_ens$ensemble, function(x){
- x$pvds[[1]]$get_point_data(Names = c("pressure", "v"))
+ x$pvds[[1]]$get_point_data(point_ids = c(0, 1, 2),
+ Names = c("pressure", "v"))
})
```
+Plotting time! Well, not quite. This time we have to think a bit more about how we plot. There's 3 parameters now and we replicated their original vectors to get our combinations so our ensemble is all over the place. We'll start by putting the data from the tibbles we extracted into a combined dataframe first like we did in the `storage`-only run.
-#### Matrix plot
+```{r}
+# Get list of dataframes
+avg_pr_dfs <- list()
+
+for(i in seq_len(length(per_por_sto_tbls))){
+ old_df <- per_por_sto_tbls[[i]]
+ new_df <- old_df %>%
+ group_by(timestep) %>%
+ summarise(avg_pressure = mean(pressure))
+
+ # Add column so we can separate the values based on their old df
+ new_df$df_id <- rep(i, nrow(new_df))
+
+ # Add column so we can plot based on old df rownames
+ new_df$orig_rwns <- as.numeric(row.names(new_df))
+ avg_pr_dfs <- c(avg_pr_dfs, list(new_df))
+}
+
+# Combine dataframes into single dataframe for plot
+avg_pr_combined_df <- bind_rows(avg_pr_dfs)
+
+```
-Maybe we want to have a look at how different parameters influence each other. For this, we can use a matrix plot like so:
+```{r}
+# # Make plot
+# ggplot(avg_pr_combined_df,
+# aes(x = orig_rwns,
+# y = avg_pressure)) +
+# geom_point(aes(color = df_id)) +
+# geom_line(aes(color = df_id)) +
+# xlab("Timestep") +
+# facet_grid(rows = vars(df_id))
+```
## Theis solution for well pumping