🟡 Temperature model (codim 2)

• 🟡 Temperature model (codim 2)
• • Continuation of Fold points
  • Using GMRES or another linear solver

This is a classical example from the Trilinos library.

This is a simple example in which we aim at solving $\Delta T+\alpha N(T,\beta)=0$ with boundary conditions $T(0) = T(1)=\beta$. This example is coded in examples/chan.jl. We start with some imports:

using Revise, BifurcationKit, LinearAlgebra, Plots
const BK = BifurcationKit

N(x; a = 0.5, b = 0.01) = 1 + (x + a*x^2)/(1 + b*x^2)

We then write our functional:

function F_chan(x, p)
	(;α, β) = p
	f = similar(x)
	n = length(x)
	f[1] = x[1] - β
	f[n] = x[n] - β
	for i=2:n-1
		f[i] = (x[i-1] - 2 * x[i] + x[i+1]) * (n-1)^2 + α * N(x[i], b = β)
	end
	return f
end

We want to call a Newton solver. We first need an initial guess:

n = 101
sol0 = [(i-1)*(n-i)/n^2+0.1 for i=1:n]

# set of parameters
par = (α = 3.3, β = 0.01)

Finally, we need to provide some parameters for the Newton iterations. This is done by calling

optnewton = NewtonPar(tol = 1e-11, verbose = true)

We call the Newton solver:

prob = BifurcationProblem(F_chan, sol0, par, (@optic _.α),
	# function to plot the solution
	plot_solution = (x, p; k...) -> plot!(x; ylabel="solution", label="", k...))
sol = @time newton( prob, optnewton)

┌─────────────────────────────────────────────────────┐
│ Newton step         residual      linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│       0     │       2.3440e+01     │        0       │
│       1     │       1.3774e+00     │        1       │
│       2     │       1.6267e-02     │        1       │
│       3     │       2.4336e-06     │        1       │
│       4     │       6.7452e-12     │        1       │
└─────────────┴──────────────────────┴────────────────┘
  0.000859 seconds (360 allocations: 2.302 MiB)

Note that, in this case, we did not give the Jacobian. It was computed internally using Automatic Differentiation.

We can perform numerical continuation w.r.t. the parameter $\alpha$. This time, we need to provide additional parameters, but now for the continuation method:

optcont = ContinuationPar(dsmin = 0.01, dsmax = 0.2, ds= 0.1, p_min = 0., p_max = 4.2,
	newton_options = NewtonPar(max_iterations = 10, tol = 1e-9))

Next, we call the continuation routine as follows.

br = continuation(prob, PALC(), optcont; plot = true)
nothing #hide

The parameter axis lens = @optic _.α is used to extract the component of par corresponding to α. Internally, it is used as get(par, lens) which returns 3.3.

You should see

scene = title!("") #hide

The left figure is the norm of the solution as function of the parameter $p=\alpha$, the y-axis can be changed by passing a different recordFromSolution to BifurcationProblem. The top right figure is the value of $\alpha$ as function of the iteration number. The bottom right is the solution for the current value of the parameter. This last plot can be modified by changing the argument plotSolution to BifurcationProblem.

Continuation of Fold points

We get a summary of the branch by doing

br

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 58
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter α starts at 3.3, ends at 4.2
 ├─ Algo: PALC
 └─ Special points:

- #  1,       bp at α ≈ +4.04103799 ∈ (+4.04103799, +4.04115545), |δp|=1e-04, [converged], δ = ( 1,  0), step =   5
- #  2,       bp at α ≈ +3.15565319 ∈ (+3.15565221, +3.15565319), |δp|=1e-06, [converged], δ = (-1,  0), step =  23
- #  3, endpoint at α ≈ +4.20000000,                                                                     step =  57

We can take the first Fold point, which has been guessed during the previous continuation run and locate it precisely. However, this only works well when the jacobian is computed analytically. We use automatic differentiation for that

# index of the Fold bifurcation point in br.specialpoint
indfold = 2

outfold = newton(
	#index of the fold point
	br, indfold)
BK.converged(outfold) && printstyled(color=:red, "--> We found a Fold Point at α = ", outfold.u.p, ", β = 0.01, from ", br.specialpoint[indfold].param,"\n")

--> We found a Fold Point at α = 3.155650731611167, β = 0.01, from 3.1556531887068466

We can finally continue this fold point in the plane $(α, β)$ by performing a Fold Point continuation. In the present case, we find a Cusp point.

outfoldco = continuation(br, indfold,
	# second parameter axis to use for codim 2 curve
	(@optic _.β),
	# we disable the computation of eigenvalues, it makes little sense here
	ContinuationPar(optcont, detect_bifurcation = 0))
scene = plot(outfoldco, plotfold = true, legend = :bottomright)

Using GMRES or another linear solver

We continue the previous example but now using Matrix Free methods. The user can pass its own solver by implementing a version of LinearSolver. Some linear solvers have been implemented from KrylovKit.jl and IterativeSolvers.jl (see Linear solvers (LS) for more information), we can use them here. Note that we can also use preconditioners as shown below. The same functionality is present for the eigensolver.

# derivative of N
dN(x; a = 0.5, b = 0.01) = (1-b*x^2+2*a*x)/(1+b*x^2)^2

# Matrix Free version of the differential of F_chan
# Very easy to write since we have F_chan.
# We could use Automatic Differentiation as well
function dF_chan(x, dx, p)
	(;α, β) = p
	out = similar(x)
	n = length(x)
	out[1] = dx[1]
	out[n] = dx[n]
	for i=2:n-1
		out[i] = (dx[i-1] - 2 * dx[i] + dx[i+1]) * (n-1)^2 + α * dN(x[i], b = β) * dx[i]
	end
	return out
end

# we create a new linear solver
ls = GMRESKrylovKit(dim = 100)

# and pass it to the newton parameters
optnewton_mf = NewtonPar(verbose = true, linsolver = ls, tol = 1e-10)

# we change the problem with the new jacobian
prob = re_make(prob;
	# we pass the differential a x,
	# which is a linear operator in dx
	J = (x, p) -> (dx -> dF_chan(x, dx, p))
	)

# we can then call the newton solver
out_mf = @time newton(prob, optnewton_mf)

┌─────────────────────────────────────────────────────┐
│ Newton step         residual      linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│       0     │       2.3440e+01     │        0       │
│       1     │       1.3774e+00     │       68       │
│       2     │       1.6267e-02     │       98       │
│       3     │       2.4336e-06     │       74       │
│       4     │       5.6275e-12     │       63       │
└─────────────┴──────────────────────┴────────────────┘
  0.003094 seconds (38.07 k allocations: 3.333 MiB, 0.05% compilation time)

We can improve this computation, i.e. reduce the number of Linear-Iterations, by using a preconditioner

using SparseArrays

# define preconditioner which is basically Δ
P = spdiagm(0 => -2 * (n-1)^2 * ones(n), -1 => (n-1)^2 * ones(n-1), 1 => (n-1)^2 * ones(n-1))
P[1,1:2] .= [1, 0.];P[end,end-1:end] .= [0, 1.]

# define gmres solver with left preconditioner
ls = GMRESIterativeSolvers(reltol = 1e-4, N = length(sol.u), restart = 10, maxiter = 10, Pl = lu(P))
optnewton_mf = NewtonPar(verbose = true, linsolver = ls, tol = 1e-10)
out_mf = @time newton(prob, optnewton_mf)

┌─────────────────────────────────────────────────────┐
│ Newton step         residual      linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│       0     │       2.3440e+01     │        0       │
│       1     │       1.3777e+00     │        3       │
│       2     │       1.6266e-02     │        3       │
│       3     │       2.3699e-05     │        2       │
│       4     │       4.8929e-09     │        3       │
│       5     │       5.6333e-12     │        4       │
└─────────────┴──────────────────────┴────────────────┘
  0.000479 seconds (552 allocations: 1.487 MiB)