Hi,
I am using DifferentialEquations.jl to solve a system of SDEs. This is a model of an electrochemical system, where I have reaction-diffusion equations.
My variables are the following (\theta, CO, OH, S_t). All of these are extended variables, so they are arrays, which I concatenate to form the varibles array u.
So far I am using the EM() solver with a fixed time-step (dt = dt_stable).
During the integration I want to ensure that my variables do not go out of some physical bound (0<\theta, CO, OH < 1, S_{t, max} < S_t < S_{t, min}).
I saw on the online documentation that I could use this solution:
isoutofdomain = (u, p, t) -> (any(x -> x < 0, p.view(u, \theta)) || any(x -> x > 1, p.view(u, \theta)) || any(x -> x < p.St_min, p.view(u, St)))
or using the GeneralDomain callback (whihc I still did not implement).
(p.view is a NamedTupleShape object that I use to unpack u )
Reading on the shortcomings of these methods I saw that they simply reject the proposed integration step and adjust the time stepping in order to not go out of bounds in the first place.
This is certainly nice but I am not sure whether the extrapolation of the solution makes sense when dealing with noise in my system.
Additionally, I want to avoid the overhead of having to check my solution at each time step and repeat the integration a number of times to ensure it.
I was thinking of simply clamping the array u after the integration step, in order to ensure my bounds. This way I would also avoid defining the differential equations outside the domain of my variables. Do you think this would be a reasonable solution ?
So I would like to ask how I can manually change the state of integrator.u right after it has been updated, before going to the next integration step and before saving the result to the output object (either the output of the solve() function or by means of a SavingCallback).