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, j in 3:p-2, k in 3:p-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?