Closed ChrisRackauckas closed 4 years ago
For those curious, we already have partial linear methods implemented:
so actually the infrastructure for doing all of this is already in there. Just the specific perform_step!
methods for specific algorithms need updates to handle AbstractDiffEqOperator
s.
Does this mean we need to implement expmv
in ApproxFun?
Hi,
Looks very nice, indeed. A couple of things:
For operator splitting, it is important that the exponentiations are placed symmetrically, like so: When time-stepping, this will then give a rate of convergence of dt^2. This is intuitively extended to multiple operators. Maybe this is already that way?
Post perform_step
, every cycle:
Multiple RHS for Crank--Nicolson (or other schemes): Usually, the wavefunction is written as a linear combination of angular basis vectors (spherical harmonics) and the Hamiltonian will then also be partitioned according to spherical symmetry: For a second-order finite-differences discretization of the kinetic operator, -laplace/2, this is a tridiagonal matrix. It turns out that , which we can exploit by
?gtsv
), with multiple right-hand sides (if an optimized implementation exists)It would be silly for me not to use this obvious fact, but I am unsure how easy this is to accomplish (edit: in a general framework, I mean).
Regarding my unclear comments on analytical asymptotics: If sufficiently far from the 'interaction region' (i.e. outside a certain radius) instead of propagating waves numerically, one can map the solution onto different plane waves exp(ik.r) which are trivial to propagate. I think this can be implemented as a custom step, along with calculating observables, etc.
Does this mean we need to implement expmv in ApproxFun?
I'm developing the exponential integrators to default to expm
but have a choice for expmv
. From this:
https://github.com/JuliaApproximation/SpectralTimeStepping.jl/blob/master/src/timeevolution.jl#L55
it looks like a dense exponential is fine because the operator is only NxN where N is the number of coefficients, right?
(the expmv
variant will take a bit of time though, mostly because we are missing the relevant pieces to do this on matrices efficiently in the Julia ecosystem, at least beyond Phi0)
For operator splitting, it is important that the exponentiations are placed symmetrically
Yes, methods which do splitting have this stuff as part of the method. Essentially the splitting stuff will be like ODEProblem((A,B),u0,tspan)
where I am given the operators from the user where they determine the splits, and the methods will then be defined appropriately symmetrically (all of the algorithms I've found use that Baker–Campbell–Hausdorff result.
Post perform_step, every cycle:
Those look like they can be handled by the callbacks. You'll need to save to observables to your own array though, but it wouldn't be difficult. That may be tutorial-worthy though.
For a second-order finite-differences discretization of the kinetic operator, -laplace/2, this is a tridiagonal matrix. Not saving more than one copy of the (sub-)Hamiltonian
Or save zero copies with a matrix-free version. Or alias the underlying matrices. I think this could be something left to the user.
tridiagonal matrix--matrix product tridiagonal matrix--vector solve (?gtsv), with multiple right-hand sides (if an optimized implementation exists)
Julia will actually do this automatically if the operator is set as Tridiagonal
. The lufact!
will specialize on Tridiagonal
matrices to do a faster factorization. Then after that's created, you can repeatedly \
the multiple RHS. Or you can use a block Krylov method which IterativeSolvers.jl is working on.
This means I should probably make sure that the DiffEqArrayOperator
works with special matrix types.
Implemented in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/1061 and continued in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/1064
It seems many physicists have problems of the form
u' = A*u
andu' = A(t)*u
. ThisA
can be from a PDE discretization, oru' = H(t)*u
naturally shows up from Schrodinger's equation. In these cases, we need to exploit the linearity of the problem in order to do well.This can all use the
DiffEqOperator
path where a user defines aDiffEqOperator
as their function. This contains anupdate_coefficients!(A,t,u)
function as well for handling time and state dependence. Then there are different things to tackle here. For one, we can make the stiff solvers "linear-aware". This can be done by, whenever we would be solving a nonlinear equation using a Newton method, instead do anIn the
#Linear Solve
part, it will just callupdate_coefficients!(A,t,u)
and then\
, so whatever overrides the user gives for how the matrix updates in time and whatever override for how the solving is supposed to take place will be used. (@shivin9 this is why the interface matching is important here)This will make things like
Trapezoid
automatically beCrankNicholson
, and will be very useful with things like BDF methods and ApproxFun. In addition, we can pull some inspiration from @jagot for Magnus integrators.http://iopscience.iop.org/article/10.1088/1367-2630/14/10/105008
Using the split operators interface we can also do things on
A + B(t)
both linear operators. But I don't know of algorithms directly on that, but technically the definition already exists.The
DiffEqOperator
will have a matrix-free pretty soon which closely matches that of LinearMaps.jl, but theDiffEqOperator
interface definesexpm
andexpmv
on the operator itself:https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/diffeq_operator.jl#L27
This means that the matrix-free versions will allow an optional method for directly defining
expmv
(and the array versions will get a quick upgrade to allow for direct definitions ofexpm
andexpmv
if someone has a use for it). This solves @jagot's problem since he wanted to use things liketo speed up the integrator by putting that definition in manually, and that interface would allow doing so.
DiffEqApproxFun.jl will get a special
AbstractDiffEqOperator
defined from ApproxFuns, which will define a linear operator onFun
types via ApproxFun operators. Even nonlinear ApproxFun operators can take this route because then it'll make it use the specialized ApproxFun Newton solver. What's interesting here is that DiffEq will be calling\
on the operator, so we can specialize the operator to hold the boundary conditions as well to do things likechop(Operator\[bcs;u],tol)
. With the type check above, this will convert the nonlinear parts on ApproxFuns to be just like SpectralTimeStepping.jl which is exactly what we needed @dlfivefifty.As for making it easy to solve finite difference PDEs, we have a library we are building which will automatically generate the linear operators from arbitrary order and derivative order PDE discretizations via Fornburg's algorithm:
https://github.com/JuliaDiffEq/PDEOperators.jl
This is all almost compatible with the full interface and will include things like time dependent BCs to make really good test cases for
u' = A(t)*u
.Things I don't know about: @jagot
What do you mean by this?