LakeComo.Rmd

---
title: "Spread of a substance of interest in Lake Como"
output: rmarkdown::html_vignette
description: >
  This document explain the basic usage of femR to solve PDEs in R. 
vignette: >
  %\VignetteIndexEntry{Spread of a substance of interest in Lake Como}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.align= "center",
                                   collapse = TRUE,
                                   comment = "#>", 
                                   fig.width= 7.3, fig.height = 5.3,
                                   message = FALSE,
                                   warning = FALSE)
                                   
library(sf, quietly = T)
library(mapview, quietly = T)
library(femR, quietly = T)
#library(rmapshaper, quietly = T)
```

Lake Como, nestled in the Lombardy region, is recognized as a crucial freshwater ecosystem esteemed for its biodiversity and ecological significance. Environmental studies have suggested the presence and spreading of a substance of interest within the lake's waters. This hypothetical presence prompts a thorough examination to understand its nature, source, and potential implications on Lake Como's delicate ecosystem.<br>
Assume that several monitoring stations have been strategically established around Lake Como to track and measure concentrations of the aforementioned substance of interest in its waters. These stations play a pivotal role in continuous surveillance, providing crucial data to assess levels of the substance of interest and its potential impacts on the lake's ecosystem.<br>
The following interactive figure shows the comprehensive map of Lake Como and includes georeferenced black points denoting the strategic placement of monitoring stations across its expanse. 
The red points denote potential 'hotspots' where, in the hypothetical scenario under investigation, abnormal concentrations of the substance of interest have been recorded. 
```{r, echo=FALSE}
data("lake_como_boundary", package="femR")
st_crs(lake_como_boundary) <- 4326
domain <- Domain(st_geometry(lake_como_boundary))
mesh <- build_mesh(domain, maximum_area = 0.25e-5, minimum_angle = 20)
mesh_sf <- st_as_sfc(mesh)
set.seed(0)
num_stations <- 20
boundary_points <- which(mesh$boundary() == 1)
boundary_id <- sample(x = 1:length(boundary_points),size = num_stations)
stations <- lapply(as.list(as.data.frame(t(mesh$nodes()[boundary_id,]))),
                   st_point)

stations <- st_sfc(stations, crs=4326)
station_name <- as.character(1:num_stations)
stations <- st_as_sf(data.frame(station_name=station_name),
                     geometry = stations)

sources <- st_as_sf(data.frame(station_name=station_name[c(16,9)]),
                    geometry=st_geometry(stations)[c(16,9)])
map.type <- "Esri.WorldTopoMap"
mappa_completa <-  mapview(stations, legend=F, cex=5, col.region="black",
        map.type=map.type) + 
  mapview(sources, legend=F, cex=5, col.region="red") +
  mapview(lake_como_boundary, col.regions = "white", lwd = 2, 
          alpha.regions=0.2,legend=F)
mappa_completa

mappa_mesh <- mapview(mesh_sf, col.region="white",alpha.regions=0.2,legend=F,  alpha=0.2,
        lwd=0.4,
        map.type=map.type) + 
  mapview(lake_como_boundary, col.regions = "white", lwd = 2, 
          alpha.regions=0.2,legend=F)

beta <- lapply(st_geometry(stations)[c(9,8)], as.numeric)
beta <- matrix(c(beta[[1]][1]-beta[[2]][1],
                 beta[[1]][2]-beta[[2]][2]), nrow=2)*10^2

mu = 0.01
Q <- 5
```

In modeling the spread of the substance of interest, it's crucial to consider the presence of a tributary in the northern region of the lake. Consequently, a transport term has been identified to represent the movement of the substance of interest carried by the inflow to other areas within Lake Como. Furthermore, non-homogeneous Dirichlet boundary conditions are applied to boundary regions near the stations that recorded anomalies in the substance of interest concentration. Conversely, homogeneous Dirichlet boundary conditions are applied to boundary regions unaffected by the anomalies.

## Steady-state solution 
In this section, we focus on solving the steady-state problem. Let $\Omega$ represent Lake Como. We assume that the concentration of the substance of interest, denoted by $u$, spreads across $\Omega$ according to the following partial differential equation (PDE):

$$
\begin{cases}
-\mu \Delta u + \boldsymbol{\beta} \cdot \nabla u = 0 \qquad & in \ \Omega \\
 u = g|_{\partial \Omega} \qquad         & on \ \partial \Omega,
\end{cases}
$$
where $\mu \in \mathbb{R}$, $\boldsymbol{\beta} \in \mathbb{R}^2$ are the diffusion coefficient and the transport coefficient, respectively. The function $g$ used to enforce the boundary conditions consists of the sum of two Gaussian-like functions centered at the locations of the stations that detected anomalous concentrations of the substance of interest. 

The upcoming sections will illustrate how to obtain a discrete solution of the problem using `femR` package. 
```{r, eval=FALSE}
library(sf, quietly = T)
library(mapview, quietly = T)
library(femR, quietly = T)

data("lake_como_boundary", package="femR")
st_crs(lake_como_boundary) <- 4326
```
The `femR` package includes utilities  to read a `sfc` (Simple Features Collection) object, such as the boundary of the Lake, and generating meshes by specifying the `maximum_area` of each triangle and the `minimum_angle` between adjacent sides of the triangles.
```{r}
# create the physical domain
domain <- Domain(st_geometry(lake_como_boundary)) 
# create the mesh
mesh <- build_mesh(domain, maximum_area = 0.25e-5, minimum_angle = 20)
```

Furthermore, the `Mesh` class overloads the `st_as_sfc` method, enabling visualization of the mesh object on an interactive map using the `mapview` package, as demonstrated in the following code snippet.
```{r}
mesh_sf <- st_as_sfc(mesh)

map.type <- "Esri.WorldTopoMap"
mapview(lake_como_boundary, col.regions = "white", lwd = 2, 
          alpha.regions=0.2, legend=F, map.type=map.type) +
mapview(mesh_sf, col.region="white", alpha.regions=0.2, legend=F,  alpha=0.4,
        lwd=0.4, map.type=map.type) 
```

The next step involves defining the `FunctionSpace` where to find the solution of the PDE. By default, first-order finite elements are considered during the construction of `Vh`. Then, the solution of the problem is instantiated as an element of `Vh`.

```{r}
Vh <- FunctionSpace(mesh)
u <- Function(Vh)
```

Then, we define the differential operator in its strong formulation as follows:
```{r}
mu <- 0.001
beta <- c(-4.091349, -4.493403)/20
Lu <- -mu*laplace(u) + dot(beta, grad(u))
```

Both the forcing term of the differential problem and the boundary condition can be defined either as standard R `function`s or as simple numerical values, particularly when representing constant values. 
```{r}
# Dirichlet boundary conditions 
g <- function(points){
  res <- matrix(0, nrow=nrow(points), ncol=1)
  pts_list <- lapply(as.list(as.data.frame(t(points))), st_point)
  dists2 <- st_distance(st_geometry(sources),
                        st_sfc(pts_list, crs=st_crs(sources))) 
  dists2 <- matrix(as.numeric(dists2), 
                   nrow=length(st_geometry(sources)), ncol=nrow(points))/(2.5*10^3)
  for(i in 1:length(sources))
    res <- res + Q*exp(-dists2[i,]) 
  return(res)
}
```
Finally, we build a `Pde` object passing the differential operator `L`, as first parameter, and the forcing term, which is a constant in this specific case, as second parameter:
```{r, results="hide"}
# build PDE object
pde <- Pde(Lu, 0.)
# set boundary conditions
pde$set_boundary_condition(boundary_condition= g, type= "dirichlet")
```
We can compute the discrete solution of the problem calling the `solve` method:
```{r}
pde$solve()
```
It is possible to visualize the computed solution on an interactive map using the `mapview` package as the following code snippet shows.
```{r}
coeff <- apply(mesh$elements(), MARGIN=1, 
               FUN= function(edge){
                     mean(u$coefficients()[edge])
         })
U <- st_sf(data.frame(coeff = coeff), geometry= mesh_sf)
mapview(U, alpha.regions=0.8, alpha = 0, map.type=map.type)
```

```{r,eval=FALSE, include=FALSE}
library(webshot2)
library(leaflet)
cntr_crds <- c(mean(mesh$nodes()[, 1]),
               mean(mesh$nodes()[, 2]))

solution <- mapview(U, alpha.regions=0.8,  legend=F,
                    alpha = 0,                
                    map.type=map.type)
folder.name <- paste0(getwd(),"/../tests/")

html_fl =  paste0(folder.name,"mappa_mesh.html")
png_fl = paste0(folder.name, "mappa_mesh.png")
mappa_mesh@map <- mappa_mesh@map %>% setView(cntr_crds[1], cntr_crds[2], zoom = 10)
mapshot2(mappa_mesh, url = html_fl, file = png_fl, delay=10, 
         cliprect = c(190,135, 600,500), zoom=1.5 )


html_fl = paste0(folder.name, "mappa_completa.html")
png_fl = paste0(folder.name, "mappa_completa.png")
mappa_completa@map <- mappa_completa@map %>% setView(cntr_crds[1], cntr_crds[2], zoom = 10) 
mapshot2(mappa_completa, url = html_fl, file = png_fl, delay=10, 
         cliprect = c(190,135, 600,500), zoom=1.5)

rect_coords <- data.frame(lng1=9.24706, lat1=46.17883,
                          lng2=9.42630, lat2=46.06229)

mappa_rect <- mappa_completa
mappa_rect@map <- mappa_rect@map %>% addRectangles(lng1=rect_coords$lng1,
                                     lat1=rect_coords$lat1,
                                     lng2=rect_coords$lng2,
                                     lat2=rect_coords$lat2, weight=2,
                                     opacity=1, group="rect",
                                     color = "red", 
                                     fillColor = "transparent")
html_fl = paste0(folder.name, "mappa_rect.html")
png_fl = paste0(folder.name, "mappa_rect.png")
mapshot2(mappa_rect, url = html_fl, file = png_fl, delay=10, 
         cliprect = c(190,135, 600,500), zoom=1.5)


html_fl = paste0(folder.name, "solution_zoom.html")
png_fl = paste0(folder.name, "solution_zoom.png")
solution@map <- solution@map %>% setView(as.numeric(st_geometry(stations)[9][[1]])[1], 
                                         as.numeric(st_geometry(stations)[9][[1]])[2], zoom = 12)
mapshot2(solution, url = html_fl, file = png_fl, delay=10, 
         cliprect = c(190,135, 600,500), zoom=1.5)

html_fl = paste0(folder.name, "solution.html")
png_fl = paste0(folder.name, "solution.png")
solution@map <- solution@map %>% setView(as.numeric(st_geometry(stations)[17][[1]])[1], 
                                         as.numeric(st_geometry(stations)[17][[1]])[2], zoom = 10)
mapshot2(solution, url = html_fl, file = png_fl, delay=10, 
         cliprect = c(190,145, 600,500), zoom=1.5)

```

## Time-dependent solution

Now, we aim to estimate the evolution of substance of interest concentrations in Lake Como. We consider a Gaussian-like function $g$ from the previous section, appropriately modified to account for the presence of time. This function will be used to define both the initial condition for the substance of interest concentration and to impose the boundary conditions. Therefore, the partial differential equation modeling the evolution of the substance of interest concentration in Lake Como is as follows:

$$
\begin{cases}
\frac{\partial u}{\partial t} -\mu \Delta u + \boldsymbol{\beta} \cdot \nabla u = 0 \qquad & in \ \Omega \times (0,T) \\
 u = g            & in \ \Omega, \ t=0  \\
 u = g|_{\partial \Omega} \qquad         & on \ \partial \Omega \times (0,T),
\end{cases}
$$

The following window wraps all the steps needed to estimate the evolution of the substance of interest concentration relying on `femR` package.
```{r, results='hide'}
time_interval <- c(0., 0.3)
times <- seq(time_interval[1], time_interval[2], length.out=50)
# function space
Vh <- FunctionSpace(mesh %X% times)
# PDE solution
u <- Function(Vh)
# differential operator
Lu <- dt(u) - mu*laplace(u) + dot(beta, grad(u))
# dirichlet bc
g <- function(points, times){
  res <- matrix(0, nrow=nrow(points), ncol=length(times))
  pts_list <- lapply(as.list(as.data.frame(t(points))), st_point)
  dists2 <- st_distance(st_geometry(sources),
                        st_sfc(pts_list, crs=st_crs(sources))) 
  dists2 <- matrix(as.numeric(dists2), 
                   nrow=length(st_geometry(sources)), ncol=nrow(points))/(2.5*10^3)
  for(i in 1:length(sources)){
      res[,1:length(times)] <- res[,1:length(times)] + Q*exp(-dists2[i,])
  }
  return(res)
}
# initial condition
u0 <- g(mesh$nodes(), times[1])
# build PDE object
pde <- Pde(Lu, 0.)
# set initial conditions
pde$set_initial_condition(u0)
# set boundary conditions
pde$set_boundary_condition(boundary_condition= g, type= "dirichlet")
# compute the discrete solution
pde$solve()
```

```{r}
# plot 
plot(u)
```