Using OrdinaryDiffEQ Integrator interface for split step

erny123 · November 6, 2024, 12:50am

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

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

step with linear-explicit: step!(linear_int)
re-initialize nonlinear integrator: reinit!(nonlinear_int, linear_int.u; t0=linear_int.t)
take nonlinear step: step!(nonlinear_int)
re-initialize linear-explicit: reinit!(linear_int, nonlinear_int.u; t0=nonlinear_int.t)
repeat loop

Not sure if this would be any good or if there is an easier way to do this

erny123 · November 6, 2024, 2:42pm

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.

Oscar_Smith · November 6, 2024, 2:48pm

What happens if you switch the implicit solver to Rosenbrock23/ Rodas5P? ImplicitEuler is very rarely what you want to use.

erny123 · November 6, 2024, 3:55pm

Yup. I switched it to ImplicitEuler after any implicit Runge-Kutta schemes were taking far too long.

I’ll try other solvers.

Oscar_Smith · November 6, 2024, 3:57pm

Also, is there a way to solve this with IMEX solvers?

erny123 · November 6, 2024, 5:09pm

I’ve tried IMEX but not sure how to tell it to only solve the nonlinear parts with an implicit method.

Oscar_Smith · November 6, 2024, 5:17pm

You want a Split ODE Problems · DifferentialEquations.jl

ChrisRackauckas · November 6, 2024, 5:21pm

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.

erny123 · November 6, 2024, 5:23pm

The point is that I can’t get the SplitODE functions which is why I’m doing this in the first place

erny123 · November 6, 2024, 6:05pm

Ok finally got some code to share:

Here’s my IMEX function

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:

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:

@btime IMEXintegration!($Solut, $integrator_lin, $integrator_nonlin,$semi_lin,100, -1560831, false)

with output:

34.888 ms (152945 allocations: 110.60 MiB)

And the Trixi output with via callbacks


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
────────────────────────────────────────────────────────────────────────────────────────────────────

erny123 · November 6, 2024, 6:21pm

Some more benchmarking:

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

erny123 · November 7, 2024, 12:55pm

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