Closed Oriane-Foussadier closed 3 months ago
rveltz
commented
1 year ago A first remark: why do you use EigArpack(1.01, :LM).
saveToFile does not help here. I would say that Field is not enough optimized, it allocates too much
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
rveltz
commented
1 year ago Did you solve it?
Oriane-Foussadier
commented
1 year ago Hello, yes indeed it was exactly the EigArpack solver. I just put again the default one and I don't have problems anymore. I just wonder why I didn't had much problem before and it only appeared in the last month...
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). Do you have any hint ?
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
Hello, yes indeed it was exactly the EigArpack solver. I just put again the default one and I don't have problems anymore. I just wonder why I didn't had much problem before and it only appeared in the last month...
For BK, I had some issues with Arpack and needed to freeze it to version Arpack = "=0.5.3" (see the project.toml)
rveltz
commented
1 year ago 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 multicontinuation. I will push this new branching procedure, but for now, here is what you can do.
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.
--- ```julia 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) end ``` ---
I 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 getFirstPointsOnBranch) and send it to multicontinuation. However, the deflation in getFirstPointsOnBranch seems to prevent convergence. Either the norm is not good, or your equations are very stiff.
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
Oriane-Foussadier
commented
1 year ago 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= $\partial{xx}$ V + Z * $\partial{x}$ C(2C - 3C^2) ).
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!
rveltz
commented
1 year ago 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= V + Z * C(2C - 3C^2) ).
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)
$$\lambda Bu=Ju $$
where $J$ is the jacobian of your rhs and $B=diag(I,0)$, the zero block being for $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 :
The equations depend on the expression of the eigenvectors. I forgot to tell you that I used eig = EigArpack(sigma = 0.1, which = :LM) and that Arpack needs to be this version:
⌅ [7d9fca2a] Arpack v0.5.3And 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.
res = predictor(bp, Val(:exhaustive), δp0; verbose = true, nbfailures=50, ampfactor = 4e-2, optn = NewtonPar(verbose = true, tol = 1e-11, maxIter = 34)) dos not work?
rveltz
commented
9 months ago 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 ?