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