Out of memory for 1D periodic non-linear system #109
Comments
|
A first remark: why do you use saveToFile does not help here. I would say that julia> @time Field(sol0,Pars);
0.000577 seconds (181 allocations: 94.312 KiB)Maybe https://docs.sciml.ai/Overview/stable/showcase/gpu_spde/#Some-Optimizations But even though, on Mac, the julia 1.9 garbage collector kicks in and the memory never exceed 9Go |
|
Did you solve it? |
|
Hello, yes indeed it was exactly the EigArpack solver. I might have another question, I am trying to find the branch originated from the first point of the following system (pretty similar also to the previous one), but with the automatic branch switching either it doesn't find another solution, or everytime come back to the already known HSS solution. I also tried the deflation operator and followed the same kind of reasoning than the 2d Swift-Hohenberg tutorial because if my parameter λ=1, I expect a slanted snaking (and for λ=0, at least a subcritical branch with a localized pattern). using Revise, Parameters, Setfield, Plots, LinearAlgebra, Random, Roots, NumericalIntegration, UnPack, DiffEqOperators, SparseArrays, JLD2, BifurcationKit
const BK = BifurcationKit
fL=1/2
fN=1.
λ = 0. #or 1
L = fL*10*π # System size
Nx = Int64(fN*512) # Number of nodes
Δx = L / Nx # Grid spacing
x = [-L/2 + i*Δx for i = 0:Nx-1] # Space vector
Pars=(Dc = 1e-2, Da= 1e-1, Di= 1, A = 1, Kd = 1, Ω0 = 0.5, Ω = 9., Ωd = 10., Z = 10., β = 4.,
Dx = CenteredDifference(1 ,3 , Δx, Nx) * PeriodicBC(Float64),
Dxx = CenteredDifference(2 ,3 , Δx, Nx) * PeriodicBC(Float64))
norminf(f) = norm(f, Inf)
########
function Field!(f, f0, Pars)
@unpack Dc, Da, Di, A, Kd, Ω0, Ω, Ωd, Z, β, Dx, Dxx = Pars
C = @view f0[1:Nx]
Na = @view f0[Nx+1:2Nx]
Ni = @view f0[2Nx+1:3Nx]
V = @view f0[3Nx+1:4Nx]
f[1:Nx] .= Dc .* (Dxx * C) .- (Dx * V) .* C .- (Dx * C) .* V .+ A .* Na .- Kd .* C
f[Nx+1] = Nx - sum(Na .+ Ni)
f[Nx+2:2Nx] .= Da * (Dxx * Na)[2:end] .- ((Dx * V).* Na)[2:end] .- ((Dx *Na).* V )[2:end] .+ Ω0* ((1 .+ λ .* Ω .* ( Na .^2)).* Ni)[2:end] .- λ .* Ωd.*( C.* Na)[2:end] .- (1-λ)*β .* Na[2:end]
f[2Nx+1:3Nx] .= Di * (Dxx * Ni) .- ((Dx * V).* Ni) .- ((Dx *Ni).* V ) .- Ω0* (1 .+ λ .* Ω .* ( Na .^2)).* Ni .+ λ .* Ωd.* C.* Na .+ (1-λ)*β .* Na
f[3Nx+1:4Nx] .= Z .* (Dx * C) .* ( 2*C - 3*C .^2) .+ Dxx * V .-V
return f
end
Field(f0,Pars)=Field!(similar(f0),f0,Pars)
######Guess 0
V0 = 0.
Eqn0(Na0) = Pars.Ω0 * (1 - Na0) * (1 + λ * Pars.Ω * Na0^2) - λ * (Pars.Ωd * Pars.A / Pars.Kd) * Na0^2 - (1 - λ) * Pars.β * Na0
Na0 = find_zero(Eqn0, (0, 1))
Ni0 = 1-Na0
C0 = (Pars.A / Pars.Kd) * Na0
sol0 = vcat(C0.* ones(Nx), Na0.* ones(Nx), Ni0.* ones(Nx), V0.* ones(Nx))
######
Pars_short=@set Pars.Z=42.
probBif = BK.BifurcationProblem(Field, sol0, Pars_short, (@lens _.Z) ;
plotSolution = (field, Pars; kwargs...) -> (plot!(field; label="", kwargs... )),
recordFromSolution = (field, Pars) -> (Cmean = integrate(x,field[1:Nx])/L , C=field[1:Nx], Na=field[Nx+1:2*Nx], Ni=field[2Nx+1:3Nx],V=field[3Nx+1:4Nx]))
opt_newton = NewtonPar(tol = 1e-10, maxIter = 20)
opts_br_eq = ContinuationPar(dsmin = 1e-15, dsmax = 0.1, ds = 0.05,
pMin = 0. , pMax = 42.5 , detectBifurcation = 3,
nev = 50,
newtonOptions = opt_newton,
maxSteps = 5000,
nInversion = 10,
detectFold = true,
maxBisectionSteps = 50,
plotEveryStep = 20,
saveEigEveryStep = 1,
saveEigenvectors = 1,
saveSolEveryStep = 0,
saveToFile = false)
br = continuation(probBif, PALC(tangent = BK.Secant(), θ=0.45), opts_br_eq, verbosity=0, plot = true, normC = norminf)
br1 = continuation(br, 1, setproperties(opts_br_eq; ds = -0.005, detectBifurcation = 3, plotEveryStep = 5, maxSteps = 2000, saveEigenvectors=true); nev = 30, plot = true, verbosity = 2)
#I tried other combinaisons of parameters, also changing theta... |
|
For the first part of your question
For BK, I had some issues with Arpack and needed to freeze it to version Arpack = "=0.5.3" (see the project.toml) |
|
You are right, it seems the random initial guess is not appropriate for your equations. In this case, we can use an approach similar to pde2path. We have to modify the current scheme for We first compute the normal form julia> bp = getNormalForm(br, 1, scaleζ = norminf)
Non simple bifurcation point at Z ≈ 42.133386352488884.
Kernel dimension = 2
Normal form:
+ 0.0131 * x1 ⋅ p + 24.1922 ⋅ x1³ + -0.9618 ⋅ x1² ⋅ x2 + 24.192 ⋅ x1 ⋅ x2² + 0.0122 ⋅ x2³
+ 0.0131 * x2 ⋅ p + -0.0005 ⋅ x1³ + 24.1935 ⋅ x1² ⋅ x2 + -0.9747 ⋅ x1 ⋅ x2² + 24.1917 ⋅ x2³The goal is then to find the zeros of this polynomial and this is where the current predictor fails. Using this method
function predictor(bp::BK.NdBranchPoint, ::Val{:exhaustive}, δp::T;
verbose::Bool = false,
ampfactor = T(1),
nbfailures = 30,
maxiter = 100,
perturb = identity,
J = nothing,
normN = x -> norm(x, Inf),
optn = NewtonPar(maxIter = maxiter, verbose = verbose)) where T
# dimension of the kernel
n = length(bp.ζ)
# initial guesses for newton
cis = Iterators.product((-1:1 for _= 1:n)...)
# find zeros of the normal on each side of the bifurcation point
function getRootsNf(_ds)
deflationOp = BK.DeflationOperator(2, 1.0, [zeros(n)]; autodiff = true)
prob = BifurcationProblem((z, p) -> perturb(bp(Val(:reducedForm), z, p)),
zeros(n), _ds, @lens _)
if ~isnothing(J)
@set! prob.VF.J = J
end
failures = 0
# we allow for 30 failures of nonlinear deflation
for ci in cis
prob.u0 .= [ci...] * ampfactor
outdef1 = BK.newton(prob, optn; normN = normN)
@debug _ds ci outdef1.converged prob.u0 outdef1.u
if converged(outdef1)
push!(deflationOp, outdef1.u)
else
failures += 1
end
end
return deflationOp.roots
end
rootsNFm = getRootsNf(-abs(δp))
rootsNFp = getRootsNf(abs(δp))
println("\n──> BS from Non simple branch point")
printstyled(color=:green, "──> we find $(length(rootsNFm)) (resp. $(length(rootsNFp))) roots before (resp. after) the bifurcation point counting the trivial solution (Reduced equation).\n")
return (before = rootsNFm, after = rootsNFp)
endI managed to find: julia> δp0 = -.0025
julia> res = predictor(bp, Val(:exhaustive), δp0; verbose = true, nbfailures=50, ampfactor = 4e-2, optn = NewtonPar(verbose = true, tol = 1e-11, maxIter = 34))
──> BS from Non simple branch point
──> we find 10 (resp. 1) roots before (resp. after) the bifurcation point counting the trivial solution (Reduced equation).
(before = [[0.0, 0.0], [-0.00018674208900862715, -0.0011523382233580843], [0.0009313929572914302, -0.0006791163572741009], [0.0009072142399203207, -0.0007106973667922095], [-0.0011636888712859906, -7.860516301052393e-6], [0.0, 0.0], [0.0011634364099430197, -4.8339594135759534e-6], [-0.0008269646828445703, 0.0008021977527069476], [-0.0006581778938610164, 0.0009463992266879826], [0.0004907732161974406, 0.0010650197930304243]], after = 0)We now need to convert to solutions of your PDE (using Using the version of BK on master, you can do pts = BK.getFirstPointsOnBranch(br, bp, res, opts; δp = δp0, usedeflation = false,)
opts_cont2 = @set opts.newtonOptions.tol=1e-9
brbp = @time BK.multicontinuation(br, bp, pts.before, pts.after, setproperties(opts_cont2; dsmax = 0.035, ds = 3e-3, dsmin = 1e-5, maxSteps = 250, detectBifurcation = 0, nInversion = 2, detectLoop = false, plotEveryStep = 2, saveSolEveryStep = 1, a =0.8);
plot = true, verbosity = 3,
δp = δp0,
verbosedeflation = true,
normC = norm,
)and get
|
|
Hello okay I see. Thank you already for all your help! However, regarding that I have several more questions then. The first one is : is it okay to use the normal form for a system like mine, in the sens that it represents in fact three chemical species evolving in time (C, Na, Ni) but the last equation is a force balance equation for the velocity V which does not contain a time derivative (V= The second, if it is okay to use this method, I wounder how you found these normal form equations because when I use the same formalism as yours, this is what I get : julia> bp = getNormalForm(br, 1, scaleζ = norminf)
Non simple bifurcation point at Z ≈ 42.133386352488884.
Kernel dimension = 2
Normal form:
+ 0.0131 * x1 ⋅ p + 24.1925 ⋅ x1³ + 0.96 ⋅ x1² ⋅ x2 + 24.1929 ⋅ x1 ⋅ x2² + -0.0002 ⋅ x2³
+ 0.0131 * x2 ⋅ p + -0.0014 ⋅ x1³ + 24.1935 ⋅ x1² ⋅ x2 + 0.9588 ⋅ x1 ⋅ x2² + 24.1917 ⋅ x2³They are different from yours and when I pass the following predictor, I don't obtain a new solution : julia> δp0 = -.0025
julia> res = BK.predictor(bp, δp0; verbose = true, nbfailures=50, ampfactor = 4e-2, optn = NewtonPar(verbose = true, tol = 1e-11, maxIter = 34))
──>BS from Non simple branch point
──> we find 1 (resp. 1) roots before (resp. after) the bifurcation point counting the trivial solution (Reduced equation). (before = [[0.0, 0.0]], after = [[0.0, 0.0]])Did you changed something before to obtain your normal form equations? And so the last interrogation is about your predictor function, I can't directly pass it like the one you wrote (because of the argument :exhaustive), I am not sure to understand how it works. Thank you again for your help! |
You mean that V is solution of a linear PDE? If the solution is unique and differentiable wrt to the Z, C, you should be able to convert your equations to a Cauchy problem, so it should be fine (but the devil is in the details). Even the principle of linear stability is not guaranteed then ;) Well, is is a kind of DAE so stability is given by the general eigenvalue problem (see tutorial) where
The equations depend on the expression of the eigenvectors. I forgot to tell you that I used
|
|
Hi! Should ew close this? |
Hello,
I use a 1D periodic non-linear system of 4 variables with a global (integral) scalar condition linking two of my variables (same kind of system as #69) and have an out of memory error.
I then have 4 fields discretized on 256 sites and ask for 40 eigenvalues at each step. Even without saving the eigenvectors or trying to detect fold points, after only 1000 steps, I have 32GB of RAM... (doesn't change if I use the saveToFile option or not).
For example in this part, when I run first "br" on a small region of the parameter (only to have the first bifurcation point) then "br1", the memory of my computer is rapidly full.
Do you have a particular advice for this ?
The text was updated successfully, but these errors were encountered: