# Using OrdinaryDiffEQ Integrator interface for split step
**URL:** https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313
**Category:** General Usage
**Tags:** ode, math, ordinarydiffeq
**Created:** [November 6, 2024, 12:50am UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313 "2024-11-06T00:50:10Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 6, 2024, 12:50am UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/1 "2024-11-06T00:50:10Z")
hi all,
Looking to see if someone can help me out.
Currently working with Trixi.jl and it’s great. However, because it’s fairly new, there are some things that have not been implemented. For example, handling nonlinear source terms like those in the Optical Kerr effect.
Trixi.jl will return a semi-discretized Discontinuous Galerkin ODE problem with the hyperbolic part of the equation along with the nonlinear source term. Currently, there is no way to split it up so I can’t use the `SplitODEProblem`.
My next thought was to use the `integrator` interface: [Integrator Interface · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)
Here’s what I’m thinking:
- First set up the linear `Trixi` problem with just the hyperbolic part with the nonlinear source term set to zero. This will return a semi-discretized `ODEProblem`.
- Then set up a nonlinear `Trixi` problem where the hyperbolic part is set to zero but the nonlinear source term is not zero.
- Create separate integrators for the linear and nonlinear problems using `linear_int = init( linear_prob, Explicit_Algorithm )` and `nonlinear_int = init( nonlinear_prob, Implicit_Algorithm)`
**Start loop**
1. step with linear-explicit: `step!(linear_int)`
2. re-initialize nonlinear integrator: `reinit!(nonlinear_int, linear_int.u; t0=linear_int.t)`
3. take nonlinear step: `step!(nonlinear_int)`
4. re-initialize linear-explicit: `reinit!(linear_int, nonlinear_int.u; t0=nonlinear_int.t)`
5. repeat loop
Not sure if this would be any good or if there is an easier way to do this
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 6, 2024, 2:42pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/2 "2024-11-06T14:42:37Z")
So looked like this worked. Using the ` SSPRK53` solver for the linear term and the `ImplicitEuler` for the nonlinear term.
The only issue is that it still seems fairly slow and I’m not sure if there’s a way to optimize this.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [November 6, 2024, 2:48pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/3 "2024-11-06T14:48:56Z")
What happens if you switch the implicit solver to `Rosenbrock23`/ `Rodas5P`? `ImplicitEuler` is very rarely what you want to use.
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 6, 2024, 3:55pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/4 "2024-11-06T15:55:23Z")
Yup. I switched it to `ImplicitEuler` after any implicit Runge-Kutta schemes were taking far too long.
I’ll try other solvers.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [November 6, 2024, 3:57pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/5 "2024-11-06T15:57:12Z")
Also, is there a way to solve this with IMEX solvers?
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 6, 2024, 5:09pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/6 "2024-11-06T17:09:14Z")
I’ve tried IMEX but not sure how to tell it to only solve the nonlinear parts with an implicit method.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [November 6, 2024, 5:17pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/7 "2024-11-06T17:17:30Z")
You want a [Split ODE Problems · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/dev/types/split_ode_types/#split_ode_prob)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 6, 2024, 5:21pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/8 "2024-11-06T17:21:48Z")
> [@erny123](#):
>
> The only issue is that it still seems fairly slow and I’m not sure if there’s a way to optimize this.
Likely the performance of this with a PDE is heavily dominated by the linear algebra, so you need to make sure you’re using a matrix-free Newton-Krylov technique or defining the sparsity pattern of the nonlinear operator.
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 6, 2024, 5:23pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/9 "2024-11-06T17:23:16Z")
The point is that I can’t get the SplitODE functions which is why I’m doing this in the first place
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 6, 2024, 6:05pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/10 "2024-11-06T18:05:14Z")
Ok finally got some code to share:
Here’s my IMEX function
```julia
function IMEXintegration!(solMat::AbstractArray, linear_int, nonlinear_int, semi, steps::Int64, saveevery::Int64, printout::Bool )
println("Beginning IMEX Integration")
println("-------------------------------")
j = 1 #index to which to add to solMat
printint = Int(round(steps/100))
for i in 1:steps
if(printout)
if mod(i, printint) == 0.0
println("Step: ", i)
end
end
if mod(i, saveevery) == 0.0
solMat[:,:,j] .= PlotData1D(linear_int.u, semi).data
j = j+1
end
step!(linear_int)
reinit!(nonlinear_int, linear_int.u; t0=linear_int.t)
step!(nonlinear_int)
reinit!(linear_int, nonlinear_int.u; t0=nonlinear_int.t)
end
end
```
And setting up the integrators:
```julia
integrator_lin = init(ode_lin, SSPRK53(); abstol = time_int_tol, reltol = time_int_tol, dt = dT, save_everystep = false, saveat = visnodes, callback = callbacks)
integrator_nonlin = init(ode_nonlin, Rosenbrock23(autodiff=false,linsolve = KrylovJL_GMRES()); abstol = time_int_tol, reltol = time_int_tol, dt = dT)
```
and so me simple benchmarking:
```julia
@btime IMEXintegration!($Solut, $integrator_lin, $integrator_nonlin,$semi_lin,100, -1560831, false)
```
with output:
```julia
34.888 ms (152945 allocations: 110.60 MiB)
```
And the `Trixi` output with via callbacks
```julia
Step: 154
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearAdvectionEquation4var1D' with DGSEM(polydeg=4)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 20153 run time: 1.87000000e-07 s
Δt: 5.00000000e-05 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 2.01535000e+00 (50.384%) time/DOF/rhs!: 1.07176823e-07 s
PID: Inf s
#DOFs per field: 1280 alloc'd memory: 500.469 MiB
#elements: 256
Variable: Asr Asi A0r A0i
L2 error: 1.44485280e-02 1.17063821e-04 5.39406253e-01 4.37033018e-03
Linf error: 1.00383538e-01 8.13320259e-04 1.00404289e+00 8.13486849e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearAdvectionEquation4var1D' with DGSEM(polydeg=4)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 20154 run time: 1.95000000e-07 s
Δt: 5.00000000e-05 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 2.01540000e+00 (50.385%) time/DOF/rhs!: 1.37684896e-08 s
PID: Inf s
#DOFs per field: 1280 alloc'd memory: 501.597 MiB
#elements: 256
Variable: Asr Asi A0r A0i
L2 error: 1.44485280e-02 1.17063821e-04 5.39406253e-01 4.37033018e-03
Linf error: 1.00383538e-01 8.13320259e-04 1.00404289e+00 8.13486849e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
```
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 6, 2024, 6:21pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/11 "2024-11-06T18:21:46Z")
Some more benchmarking:
```julia
funca(a) = step!(a)
funcb(a,b) = reinit!(a,b.u; t0=b.t)
#compile first
funca(integrator_lin)
funcb(integrator_nonlin,integrator_lin)
```
and benchmark:
`@btime funca($integrator_lin)`
` 94.646 μs (700 allocations: 22.58 KiB)`
`@btime funcb($integrator_nonlin, $integrator_lin)`
`23.629 μs (148 allocations: 5.08 KiB)`
`@btime funca($integrator_nonlin)`
`146.916 μs (674 allocations: 1.08 MiB)`
`@btime funcb($integrator_lin, $integrator_nonlin)`
`1.630 μs (7 allocations: 240 bytes)`
It looks like there’s some significant time taken to reinitialize the initial value `u0` between the linear and nonlinear steps. Obviously the best way would be for them to share the `u` values
EDIT: Also there seems to be a good amount of allocations that could be easily allocated initially since they’re not much larger than a couple of MiB
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 7, 2024, 12:55pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/12 "2024-11-07T12:55:51Z")
So there’s something interesting happening. There is some numerical artifact that starts right in front of the right pulse depending on what `dt` I use. Setting it to a specific fairly high value seems to get rid of it. Anyone have any idea what’s going on?
**Initial Conditions**
**Simulations with Different** `dt`
`dt = 0.0001`
`dt = 0.001`
`dt = 0.002`
`dt = 0.003`
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [November 7, 2024, 11:06pm UTC](https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313/13 "2024-11-07T23:06:20Z")
I made a notebook for anyone that wants to play around with this:
[TRIXIExamples/TRIXI-Coupled-Nonlinear.ipynb at main · erny123/TRIXIExamples](https://github.com/erny123/TRIXIExamples/blob/main/TRIXI-Coupled-Nonlinear.ipynb)