Adding more user controls

[erny123]

Hi all, sure you've seen my slack messages.

Trixi.jl has a great user interface and it's very flexible. However, anyone who is visiting the Trixi.jl page and is actually using it for research is at least familiar with the basics of Spectral Elements and the Discontinuous Galerkin Method and will probably be extremely familiar with numerical methods (linear and nonlinear solvers). I think the one thing that Trixi.jl is missing is the flexibility and malleability to customize the solving algorithms for each term. For example, if one wants to treat a linear terms in a different manner than nonlinear. We also need more flexibility in how to individually treat the hyperbolic and parabolic terms.
As much as I like the DifferentialEquations.jl , it assumes far too much and it is not flexible for new problems and equations, at least from a surface user view.

Suggestions:

• Julia has LinearSolve.jl and NonlinearSolve.jl . Allowing the flexibility to use these packages will be extremely useful
• For a linear term to be solved explicitly,semidiscretization() should give either the matrix operator and it's inverse, a matrix free operator and its matrix free inverse (this can either be LinearMaps.jl or LinearOperators.jl . For a term that needs to be solved implicitly (usually nonlinear but could be parabolic terms) it should return a semidiscretization() should return a nonlinear function to be used in SciMLBase.NonlinearProblem or SciMLBase.NonlinearFunction

Dr. Doering suggested it be part of rhs! and a rhs_parabolic!. Additionally, matrix formulations are definitely costly and I think more people would be interested in iterative solvers like Krylov.jl

[ranocha]

I may not understand everything correctly. What exactly would you like to achieve? Could you please provide an example that is as simple as possible, e.g., for a 1D linear advection equation?

[JoshuaLampert]

Xref https://discourse.julialang.org/t/using-ordinarydiffeq-integrator-interface-for-split-step/122313

[erny123]

@ranocha I made a simple 1D example of 2 coupled strongly nonlinear advection equations:

https://github.com/erny123/TRIXIExamples/blob/main/TRIXI-Coupled-Nonlinear.ipynb

You can see in the first solution using the regular method, that OrdinaryDiffEq.jl is not good enough for solving. The second solution creates a linear Trixi equation and a nonlinear Trixi equation. I then use the integrator interface to solve the linear part with SSPRK and the nonlinear part with an implicit method. The second solution is much faster and correct.

The point is, there needs to be more flexibility for users to solve the terms of the semidiscretized equations with their own algorithms.

[ranocha]

You mentioned that everything was very slow on Discourse. My first impression is that this is caused by non-constant global variables like c0 used in pumpgaussinit etc. called from BCLeftRight. See https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-untyped-global-variables

[ranocha]

If you want to use use a SplitODEProblem and corresponding solver from OrdinaryDiffEq.jl, you can basically use your semi_lin and semi_nonlin for it. Since the code you have posted is not minimal but quite a lot, I will only give some pseudocode how to set it up:

function rhs_nonlin!(du_ode, u_ode, parameters, t)
    semi = parameters.semi_nonlin
    Trixi.rhs!(du_ode, u_ode, semi, t)
end

function rhs_lin!(du_ode, u_ode, parameters, t)
    semi = parameters.semi_lin
    Trixi.rhs!(du_ode, u_ode, semi, t)
end

parameters = (; semi_nonlin, semi_lin)
ode = SplitODEProblem(rhs_nonlin!, rhs_lin!, u0, tspan, parameters)

[erny123]

Thanks @ranocha . I'll give it a shot. Yes, the code is not minimal. I'm not sure there is a way to minimize it, however. Given the Trixi.jl interface, the issue is that there is no support for complex numbers. Given the equation is just a 2-variable complex advection equation with a nonlinear coupling source term, this was the easiest example I could think of.

[ranocha]

Is the example representative of the kind of problems you need to solve or will they be more complicated (e.g., multiple space dimensions)?

[erny123]

The ultimate goal is to model 3D with a system of N coupled equations with nonlinear coupling source terms along with a nonlinear self-coupling term (intensity for the Kerr effect).

In optics, the complex representation is important because of the phase of the wave and how it effects the diffraction (diffusion) along with the nonlinear coupling terms.

Additionally, in general, Maxwell's equations in a medium are represented by a complex relative permittivity (e.g. Kramers-Kronig relation).

In a 1D setting, I would use the split step method and treat the linear part of the system with an explicit SSP method and the nonlinear part with an implicit method. I believe this is generally the standard way to solve coupled advection wave equations, but I could be wrong

[erny123]

Here's an example of the coupled nonlinear schrodinger equation:

https://w1.mtsu.edu/faculty/wding/files/ShuNLS.pdf

Although these have small couplilng factors.

and this work has done work on IMEX SSPRK methods:

https://arxiv.org/abs/1702.04621

[erny123]

Here's another reason I believe we need more controls:

Say a system has a flux equal to the convolution in time of a response function and a variable at that point in space.
For example, the flux at the current time $t_{n}$ would be something like:

$$ F(x,t) = \int_{-\infty}^{t_{n}} R(t-\tau) \: U(x, \tau) d\tau $$

where $R$ is the response function and $U$ is the variable to be solved at position $x$.

To do this convolution, I would do a Z-transform (discrete Laplace transform) of the $R$ and $U$ independently then multiply the transforms together then inverse transform:

$$ F(x,t) = Z^{-1}\{\: Z\{ R \} \: Z\{ U \} \:\} $$

Then plug it into the hyperbolic equation:

$$ \partial_{t} U(x,t) + \partial_{x} F(x,t) = 0 $$

For this, I would need access the values at each point at previous times.