What
I had a reasoning session with myself here about improving the flow integration in #3039. The conclusion I came to is that it would be very helpful to have a cubic Hermite spline of the flows $q(t)$ over the latest timestep $[t_n, t_{n+1}]$:
$$
\tilde{q}(t) = q(t_n) + \Delta t \left[
\dot{q}(t_n)\left(\frac{t - t_n}{\Delta t}\right)+
\left(3\frac{\Delta q}{\Delta t}-\dot{q}(t_{n+1}) - 2\dot{q}(t_n)\right)\left(\frac{t - t_n}{\Delta t}\right)^2 +
\left(\dot{q}(t_1)+\dot{q}(t_2) - 2\frac{\Delta q}{\Delta t}\right)\left(\frac{t - t_n}{\Delta t}\right)^3
\right].
$$
This is a third degree polynomial which fits the flows (corresponding to the solver solution) exactly at $q(t_n), q(t_{n+1}), \dot{q}(t_n), \dot{q}(t_{n+1})$. For the cumulative flows we integrate this interpolation over the timestep, which yields
$$
\int_{t_n}^{t_{n+1}} \tilde{q}(t)\text{d}t = \frac{\Delta t}{2}[q(t_n) + q(t_{n+1})] + \frac{\Delta t^2}{12}[\dot{q}(t_n) - \dot{q}(t_{n+1})].
$$
For the concentration computations we can compute cumulative flow contributions per sub-timestep by integrating the interpolation over the subtimesteps, say $[t_n + \alpha_1\Delta t, t_n + \alpha_2\Delta t]$:
$$
\int_{t_n + \alpha_1\Delta t}^{t_n + \alpha_2\Delta t} \tilde{q}(t)\text{d}t =
q(t_n)\Delta t (\alpha_2 - \alpha_1) +
\Delta t^2\left[
\frac{1}{2}\dot{q}(t_n)(\alpha_2^2 - \alpha_1^2) +
\frac{1}{3}\left(3\frac{\Delta q}{\Delta t}-\dot{q}(t_{n+1}) - 2\dot{q}(t_n)\right)(\alpha_2^3 - \alpha_1^3) +
\frac{1}{4}\left(\dot{q}(t_1)+\dot{q}(t_2) - 2\frac{\Delta q}{\Delta t}\right)(\alpha_2^4 - \alpha_1^4)
\right].
$$
This isn't as computationally expensive as it might look with the right caching and flop minimization.
How
I propose to (again) formulate the rhs as (ignoring PID control):
$$
f(s, t) = Mq(s, t) + b(t)
$$
That way:
-
$f$ can be passed to the solver as the rhs
-
$q$ can be used for the cumulative flow computations
- The Jacobian computation can be made a bit cheaper; $J(f) = MJ(q)$, so we use AD to compute the Jacobian of $q$ and do the multiplication with $M$ ourselves only once, and $b(t)$ can be ignored.
$\tilde{q}(t)$ above is an approximation of
$$
q(t) = q(s(t), t),
$$
where $s(t)$ is the storage interpolation provided by the solver. Hence, we can use ForwardDiff in combination with DifferentiationInterface.value_and_derivative! to efficiently obtain $q(t_n), \dot{q}(t_n)$ and $q(t_{n+1}), \dot{q}(t_{n+1})$ with one $q(s, t)$ call respectively.
What
I had a reasoning session with myself here about improving the flow integration in #3039. The conclusion I came to is that it would be very helpful to have a cubic Hermite spline of the flows$q(t)$ over the latest timestep $[t_n, t_{n+1}]$ :
This is a third degree polynomial which fits the flows (corresponding to the solver solution) exactly at$q(t_n), q(t_{n+1}), \dot{q}(t_n), \dot{q}(t_{n+1})$ . For the cumulative flows we integrate this interpolation over the timestep, which yields
For the concentration computations we can compute cumulative flow contributions per sub-timestep by integrating the interpolation over the subtimesteps, say$[t_n + \alpha_1\Delta t, t_n + \alpha_2\Delta t]$ :
This isn't as computationally expensive as it might look with the right caching and flop minimization.
How
I propose to (again) formulate the rhs as (ignoring PID control):
That way:
where$s(t)$ is the storage interpolation provided by the solver. Hence, we can use $q(t_n), \dot{q}(t_n)$ and $q(t_{n+1}), \dot{q}(t_{n+1})$ with one $q(s, t)$ call respectively.
ForwardDiffin combination withDifferentiationInterface.value_and_derivative!to efficiently obtain