# Efficient way to handle operator splitting (periodically updating solution mid-solve) in JuliaDiffEq

**URL:** https://discourse.julialang.org/t/efficient-way-to-handle-operator-splitting-periodically-updating-solution-mid-solve-in-juliadiffeq/30757
**Category:** Numerics
**Tags:** question, diffeq
**Created:** [November 5, 2019, 4:26pm UTC](https://discourse.julialang.org/t/efficient-way-to-handle-operator-splitting-periodically-updating-solution-mid-solve-in-juliadiffeq/30757 "2019-11-05T16:26:12Z")
**Posts on this page:** 1
**Showing post:** 1

<div class="post-metadata">

### Author: ![jimmielin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jimmielin/32/11066_2.png) [@jimmielin](https://discourse.julialang.org/u/jimmielin)
#### Post date: [November 5, 2019, 4:26pm UTC](https://discourse.julialang.org/t/efficient-way-to-handle-operator-splitting-periodically-updating-solution-mid-solve-in-juliadiffeq/30757/1 "2019-11-05T16:26:12Z")

</div>

I am writing a simple model of atmospheric chemistry in Julia using `DifferentialEquations.jl` to solve the chemical kinetics ODEs.

While chemical kinetics can be modeled using a system of ODEs, some processes in modern atmospheric models cannot be expressed in terms of differential equations, i.e. time-dependent emissions, sub-grid scale processes, etc. They can be thought as simple `if`-`else` clauses that need to be executed to mutate the solution array (species concentrations) at given time steps. Thus the ODE solve needs to be ran step-by-step in conjunction with these if-else clauses, i.e. “operator splitting”.

Thus I am trying to find a performance-efficient way to mutate the solution array in the middle of solving the ODE. e.g. the code below for a simple stratospheric mechanism

```julia
function strato_t!(du, u, p, t)
    O, O1D, O3, NO, NO2, M, O2 = u # species
    param = p(t)
    SUN = param[1]
    du[1] = 2*(2.643E-10)*SUN*SUN*SUN*O2 - 1*(8.018E-17)*O*O2 + 1*(6.120E-04)*SUN*O3 - 1*(1.576E-15)*O*O3 + 1*(7.110E-11)*O1D*M - 1*(1.069E-11)*NO2*O + 1*(1.289E-02)*SUN*NO2 # O
    du[2] = 1*(1.070E-03)*SUN*SUN*O3 - 1*(7.110E-11)*O1D*M - 1*(1.200E-10)*O1D*O3 # O1D
    du[3] = 1*(8.018E-17)*O*O2 - 1*(6.120E-04)*SUN*O3 - 1*(1.576E-15)*O*O3 - 1*(1.070E-03)*SUN*SUN*O3 - 1*(1.200E-10)*O1D*O3 - 1*(6.062E-15)*NO*O3 # O3
    du[4] = (-1)*(6.062E-15)*O3*NO + 1*(1.069E-11)*NO2*O + 1*(1.289E-02)*SUN*NO2 # NO
    du[5] = 1*(6.062E-15)*O3*NO - 1*(1.069E-11)*O*NO2 - 1*(1.289E-02)*SUN*NO2 # NO2
    du[6] = 0 # M (fix)
    du[7] = 0 # O2 (fix)
    nothing
end
function model_onestep()
    u0 = [6.624e+8, 9.906e+1, 5.326e+11, 8.725e+8, 2.240e+8, 8.120e+16, 1.697e+16]
    p = t -> [t/3600.0]
    oprob = ODEProblem(strato_t!, u0, (0.0, 3600.0), p)
    @time sol = solve(oprob, dense=false, save_on=false, calck=false, save_everystep=false)
end

```

When solved in a single call to `solve` spanning from `t = 0.0` to `t = 3600.0`, it takes

```julia
0.454232 seconds (2.65 M allocations: 147.560 MiB, 4.91% gc time)

```

(I believe there are allocation problems here, but irrelevant for now)

However I need to mutate the solution array at given steps, e.g. every `600.0` s. My current approach is to `solve` and `remake` in a loop, which is the prevalent approach in current Fortran-based models:

```julia
function model_chunked()
    u0 = [6.624e+8, 9.906e+1, 5.326e+11, 8.725e+8, 2.240e+8, 8.120e+16, 1.697e+16]
    p = t -> [t/3600.0]
    oprob = ODEProblem(strato_t!, u0, (0.0, 600.0), p)
    for t in 0.0:600.0:3600.0
        @time sol = solve(oprob, dense=false, save_on=false, calck=false, save_everystep=false)
        
        u = sol.u[end]
        # ... need to do something with the solution here
        # to mimic operator splitting
        # e.g. if(t == 1200.0) u[2] += 1.5e+1
        
        oprob = remake(oprob; u0=u, tspan=(t, t+600.0))
    end
end

```

This allocates a lot:

```julia
  0.293247 seconds (1.82 M allocations: 101.559 MiB, 6.51% gc time)
  0.284285 seconds (1.82 M allocations: 101.539 MiB, 5.27% gc time)
  0.289665 seconds (1.82 M allocations: 101.533 MiB, 6.29% gc time)
  0.290912 seconds (1.82 M allocations: 101.529 MiB, 6.90% gc time)
  0.288586 seconds (1.82 M allocations: 101.520 MiB, 5.07% gc time)
  0.294330 seconds (1.82 M allocations: 101.514 MiB, 6.53% gc time)
  0.288314 seconds (1.82 M allocations: 101.508 MiB, 6.46% gc time)
  2.101101 seconds (12.77 M allocations: 710.813 MiB, 5.94% gc time)

```

Of which the `(1.82M, 101.5 MiB)` are what I assume to be setup allocations in the solve routine. There is also a lot of garbage collection going on, which is not good for performance.

**What would be a efficient way to mutate the solution array at given “checkpoints” within the DiffEq `solve`**? Is there a straightforward way to “reuse” the temporary allocations made by `solve` since the problem formulation has not really changed?

Thank you!

---

_[View the full topic](https://discourse.julialang.org/t/efficient-way-to-handle-operator-splitting-periodically-updating-solution-mid-solve-in-juliadiffeq/30757)._