Large ManifoldProjection causing memory issues in DifferentialEquations.jl

Specific Domains Modelling & Simulations
1

Hello,

I’m solving a PDE in three spatial dimension where the field takes values in S^3. A convenient way to work with this is to define a four-vector u and project it onto the sphere at each point. So I have u_i(x,y,z), i=1,2,3,4 and want to demand that:

u_1(x,y,z)^2 + u_2(x,y,z)^2 + u_3(x,y,z)^2 + u_4(x,y,z)^2 = 1

To solve the PDE, I discretise space, so that the problem becomes an ODE at every point in space. I’m aiming for a 100^3 grid, so around a million ODEs.

The way to solve this, I think, is to use ManifoldProjection.

The constraint function looks like:

function g(resid,u,p,t)
	
	for i in 3:p[2][1]-2, j in 3:p[2][2]-2, k in 3:p[2][3]-2
		resid[1,i,j,k] = u[1,i,j,k]^2 + u[2,i,j,k]^2 + u[3,i,j,k]^2 + u[4,i,j,k]^2 - 1.0
		resid[2,i,j,k] = 0.0
		resid[3,i,j,k] = 0.0
		resid[4,i,j,k] = 0.0
	end
  
end

I read that resid should be the same size as your field, which is why I’ve got these extra zero field.

And my solver is this:

function flow(sk,p)
    
# latticespacing, latticepoints and some parameter are in p. sk is a Struct.
    p=(sk.ls,sk.lp,1.0)

    tspan = (0.0,0.01)
	
	cb2 = ManifoldProjection(g,save=false)

    prob = ODEProblem(dEdu!,sk.u,tspan,p)
	sol = solve(prob, Tsit5(),reltol=1e-8, abstol=1e-8, save_everystep=false, callback=cb2)

    return sol[end]
    
end

When I run this, the system quickly runs out of memory even for small numbers of lattice points. When the callback is off, it doesn’t (so I think the problem doesn’t like in dEdu!). I don’t think I’m saving anything extra, since I’ve used save=false.

What stupid thing am I doing?

Many thanks,
Chris

2

There’s no need to do a general projection here. Instead just define a DiscreteCallback where the condition always returns true and then divide all values by u[1,i,j,k]^2 + u[2,i,j,k]^2 + u[3,i,j,k]^2 + u[4,i,j,k]^2 - 1.0. That will avoid using Newton iterations on a huge sparse Jacobian.

3

Thanks Chris - works a treat. And thanks for all your work, DifferentialEquations.jl is very fun.

For anyone who’s reading this with a similar problem, I have some more tips: my model with just this projection was still fairly unstable because the projection only happens at the time steps used by the integrator. Increasing tolerance, decreasing your spatial step or setting your own tiny timesteps all help stability but are speed killers.

A better solution for me was to add a Lagrange Multiplier to the system. Numerically, this doesn’t exactly enforce the constraint so you should still apply the projection at each step. Doing this meant that my problem was stable without having to do anything over-the-top with the solver.

-Chris