deSolve: differential equations with two consecutive dynamics

Question

I am simulating a ring tube with flowing water and a temperature gradient using deSolve::ode(). The ring is modelled as a vector where each element has a temperature value and position.

I am modelling the heat diffusion formula:

1)

But I'm struggling with also moving the water along the ring. In theory, it's just about substituting the temperature at the element i in the tube vector with that at the element s places earlier. Since s may not be an integer, it can be separated into the integer part (n) and the fractional part (p): s=n+p. Consequently, the change in temperature due to the water moving becomes:

2)

The problem is that s equals to the water velocity v by the dt evaluated at each iteration of the ode solver.

My idea is to treat the phenomenons as additive, that is first computing (1), then (2) and finally adding them together. I'm afraid though about the effect of time. The ode solver with implicit methods decides the time step automatically and scales down linearly the unitary change delta.

My question is whether just returning (1) + (2) in the derivative function is correct or if I should break the two processes apart and compute the derivatives separately. In the second case, what would be the suggested approach?

EDIT: As by suggestion by @tpetzoldt I tried to implement the water flow using ReacTran::advection.1D(). My model has multiple sources of variation of temperature: the spontaneous symmetric heat diffusion; the water flow; a source of heat that is turned on if the temperature near a sensor (placed before the heat source) drops below a lower threshold and is turned off if raises above an upper threshold; a constant heat dispersion determined by a cyclical external temperature.

Below the "Moving water" section there is still my previous version of the code, now substituted by ReacTran::advection.1D(). The plot_type argument allows visualizing either a time sequence of the temperature in the water tube ("pipe"), or the temperature sequence at the sensors (before and after the heater).

library(deSolve)
library(dplyr)
library(ggplot2)
library(tidyr)
library(ReacTran)

test <- function(simTime = 5000, vel = 1, L = 500, thresh = c(16, 25), heatT = 25,
                                    heatDisp = .0025,   baseTemp = 15, alpha = .025,
                                    adv_method = 'up', plot_type = c('pipe', 'sensors')) {
    
    plot_type <- match.arg(plot_type)

    thresh <- c(16, 25)

    sensorP <- round(L/2)

    vec <- c(rep(baseTemp, L), 0)

    eventfun <- function(t, y, pars) {

        heat <- y[L + 1] > 0

        if (y[sensorP] < thresh[1] & heat == FALSE) { # if heat is FALSE -> T was above the threshold
            #browser()
            y[L + 1] <- heatT
        }

        if (y[sensorP] > thresh[2] & heat == TRUE) { # if heat is TRUE -> T was below the threshold
            #browser()
            y[L + 1] <- 0
        }

        return(y)
    }

    rootfun <- function (t, y, pars) {

        heat <- y[L + 1] > 0

        trigger_root <- 1

        if (y[sensorP] < thresh[1] & heat == FALSE & t > 1) { # if heat is FALSE -> T was above the threshold
            #browser()
            trigger_root <- 0
        }

        if (y[sensorP] > thresh[2] & heat == TRUE & t > 1) { # if heat is TRUE -> T was below the threshold
            #browser()
            trigger_root <- 0
        }


        return(trigger_root)
    }

    roll <- function(x, n) {
        x[((1:length(x)) - (n + 1)) %% length(x) + 1]
    }

    fun <- function(t, y, pars) {

        v <- y[1:L]

        # Heat diffusion: dT/dt = alpha * d2T/d2X
        d2Td2X <- c(v[2:L], v[1]) + c(v[L], v[1:(L - 1)]) - 2 * v

        dT_diff <- pars * d2Td2X

        # Moving water
        # nS <- floor(vel)
        # pS <- vel - nS
        #
        # v_shifted <- roll(v, nS)
        # nS1 <- nS + 1
        # v_shifted1 <- roll(v, nS + 1)
        #
        # dT_flow <- v_shifted + pS * (v_shifted1 - v_shifted) - v
        dT_flow <- advection.1D(v, v = vel, dx = 1, C.up = v[L], C.down = v[1],
                                                        adv.method = adv_method)$dC

        dT <- dT_flow + dT_diff

        # heating of the ring after the sensor
        dT[sensorP + 1] <- dT[sensorP  + 1] + y[L + 1]

        # heat dispersion
        dT <- dT - heatDisp * (v - baseTemp + 2.5 * sin(t/(60*24) * pi * 2))

        return(list(c(dT, 0)))
    }

    out <- ode.1D(y = vec, times = 1:simTime, func = fun, parms = alpha, nspec = 1,
                                    events = list(func = eventfun, root = T),
                                    rootfunc = rootfun)


    if (plot_type == 'sensors') {

        ## Trend of the temperature at the sensors levels
        out %>%
            {.[,c(1, sensorP + 1, sensorP + 3, L + 2)]} %>%
            as.data.frame() %>%
            setNames(c('time', 'pre', 'post', 'heat')) %>%
            mutate(Amb = baseTemp + 2.5 * sin(time/(60*24) * pi * 2)) %>%
            pivot_longer(-time, values_to = "val", names_to = "trend") %>%
            ggplot(aes(time, val)) +
            geom_hline(yintercept = thresh) +
            geom_line(aes(color = trend)) +
            theme_minimal() +
            theme(panel.spacing=unit(0, "lines")) +
            labs(x = 'time', y = 'T°', color = 'sensor')
    } else {

    ## Trend of the temperature in the whole pipe
    out %>%
        as.data.frame() %>%
        pivot_longer(-time, values_to = "val", names_to = "x") %>%
        filter(time %in% round(seq.int(1, simTime, length.out = 40))) %>%
        ggplot(aes(as.numeric(x), val)) +
        geom_hline(yintercept = thresh) +
        geom_line(alpha = .5, show.legend = FALSE) +
        geom_point(aes(color = val)) +
        scale_color_gradient(low = "#56B1F7", high = "red") +
        facet_wrap(~ time) +
        theme_minimal() +
        theme(panel.spacing=unit(0, "lines")) +
        labs(x = 'x', y = 'T°', color = 'T°')
    }
}

It's interesting that setting an higher number of segment (L = 500) and high speed (vel = 2) it's possible to observe a spiking sequence in the post heating sensor. Also, the processing time drastically increases, but more as an effect of increased velocity than due to increased pipe resolution.

My biggest doubt now is whether ReacTran::advection.1D() does make sense in my context since I'm modeling water temperature, while this function seems more related to the concentration of a solute in flowing water.

Answer 1

The problem looks like a PDE example with a mobile and a fixed phase. A good introduction about the "method of lines" (MOL) approach with R/deSolve can be be found in the paper about ReachTran from Soetaert and Meysman (2012) doi.org/10.1016/j.envsoft.2011.08.011.

An example PDE can be found at slide 55 of some workshop slides, more in the teaching package RTM.

R/deSolve/ReacTran tries to make ODEs/PDEs easy, but pitfalls remain. If numerical dispersion or oscillations occur, it can be caused by violating the Courant–Friedrichs–Lewy condition.