GitHub

🟠 Detonation engine

This is a model of a detonation engine, a new kind of reactor developed for planes. The model[Koch] quantifies the spatio-temporal evolution of a property analogous to specific internal energy, $u(x,t)$, on a one dimensional (1D) periodic domain:

\[\left\{ \begin{gathered} \frac{\partial u}{\partial t}=\nu_{1} \frac{\partial^{2} u}{\partial x^{2}}-u \frac{\partial u}{\partial x}+k q(1-\lambda) \exp \left(\frac{u-u_{c}}{\alpha}\right)-\epsilon u^{2} \\ \frac{\partial \lambda}{\partial t}=\nu_{2} \frac{\partial^{2} \lambda}{\partial x^{2}}+k(1-\lambda) \exp \left(\frac{u-u_{c}}{\alpha}\right) -\frac{s u_{p} \lambda}{1+\exp \left(r\left(u-u_{p}\right)\right)} \end{gathered}\right.\]

Problem discretization

We start by discretizing the above PDE based on finite differences.

using Revise, Parameters
using DiffEqOperators, ForwardDiff, DifferentialEquations, SparseArrays
using BifurcationKit, LinearAlgebra, Plots, Setfield
const BK = BifurcationKit

ω(u, p) = p.k * exp((u - p.uc) / p.α)
β(u, p) = p.s * p.up /(1 + exp(p.r * (u - p.up)) )

# utilities for plotting solutions
function plotsol!(x; k...)
	n = length(x) ÷ 2
	u = @view x[1:n]
	v = @view x[n+1:2n]
	plot!(u; label="u", k...)
	plot!(v; label="λ", k...)
end
plotsol(x; k...) = (plot();plotsol!(x; k...))

# function to build derivative operators
function DiffOp(N, lx; order = 2)
	h = lx/N

	D2 = CenteredDifference(2, order, h, N)
	D = CenteredDifference(1, order, h, N)

	Q = PeriodicBC(Float64)

	Δ = sparse(D2 * Q)[1]
	D = sparse(D * Q)[1]
	return D, Δ
end

# nonlinearity of the model
function NL!(dest, U, p, t = 0.)
	N = p.N
	u = @view U[1:N]
	λ = @view U[N+1:2N]
	dest[1:N]    .= p.q .* (1 .- λ) .* ω.(u, Ref(p)) .- p.ϵ .* u.^2
	dest[N+1:2N] .=        (1 .- λ) .* ω.(u, Ref(p)) .- λ .* β.(u, Ref(p))
	return dest
end

# function which encodes the right hand side of the PDE
@views function Fdet!(f, u, p, t = 0)
	N = p.N
	NL!(f, u, p) 						# nonlinearity
	mul!(f, p.Δ, u, p.ν1, 1) 			# put Laplacian
	f[1:p.N] .-= (p.D * u[1:N].^2) ./ 2	# add drift a la Burger
	return f
end
NL(U, p, t = 0.) = NL!(similar(U), U, p, t)
JNL(U, p, t = 0.) = ForwardDiff.jacobian(x -> NL(x, p), U)
Fdet(x, p, t = 0) = Fdet!(similar(x), x, p, t)
Jdet(x, p) = sparse(ForwardDiff.jacobian(x -> Fdet(x, p), x))
┌ Warning:
│
│
│     # Breaking changes
│     When dense output from `taylorinteg` is requested by the user,
│     the first row of the dense output polynomials matrix now represents
│     the Taylor expansion of the solution around the initial condition
│     w.r.t. the independent variable in the ODE. Previously it represented
│     the initial condition plus a Taylor1 perturbation.
│
└ @ TaylorIntegration ~/.julia/packages/TaylorIntegration/Yq9wW/src/TaylorIntegration.jl:25

We can now instantiate the model

N = 300
lx = 2pi
X = LinRange(0, lx, N)

D, Δ = DiffOp(N, lx)
_ν1 = 0.0075
# model parameters
par_det = (N = N, q = 0.5, α = 0.3, up = 0.0, uc = 1.1, s = 3.5, k = 1., ϵ = 0.15, r = 5.0, ν1 = _ν1, ν2 = _ν1, Δ = blockdiag(Δ, Δ), D = D, Db = blockdiag(D, D))

# initial conditions
u0 = 0.5ones(N)
λ0 = 0.5ones(N)
U0 = vcat(u0, λ0)
U0cons = vcat(copy(U0), 1.)

Jacobian with sparsity detection

Writing the jacobian explicitly is cumbersome. We rely on automatic differentiation to get the sparse jacobian.

# improved jacobian with sparse coloring
using SparseArrays, SparseDiffTools, Test
const L1 = copy(Jdet(rand(2N), par_det))
const colors = matrix_colors(L1)
function JlgvfColorsAD(J, u, p, colors)
	SparseDiffTools.forwarddiff_color_jacobian!(J, (out, x) -> out .= Fdet(x,p), u, colorvec = colors)
	J
end
JdetAD(x,p) = JlgvfColorsAD(L1, x, p, colors)

We are now ready to compute the bifurcation of the trivial (constant in space) solution:

# bifurcation problem
prob = BifurcationProblem(Fdet, U0, setproperties(par_det; q = 0.5), (@lens _.up); J = JdetAD,
	plot_solution = (x, p; k...) -> plotsol!(x; k...),
	record_from_solution = (x, p) -> (u∞ = norminf(x[1:N]), n2 = norm(x)))

prob = re_make(prob, params = (@set par_det.up = 0.56))

# iterative eigen solver
eig = EigArpack(0.2, :LM, tol = 1e-13, v0 = rand(2N))
eig = EigArnoldiMethod(sigma=0.2, which = BifurcationKit.LM(),x₀ = rand(2N ))

# newton options
optnew = NewtonPar(verbose = true, eigsolver = eig)
solhomo = newton(prob, optnew; normN = norminf)
optcont = ContinuationPar(newton_options = setproperties(optnew, verbose = false),
	detect_bifurcation = 3, nev = 50, n_inversion = 8, max_bisection_steps = 25,
	dsmax = 0.01, ds = 0.01, p_max = 1.4, max_steps = 1000, plot_every_step = 50)

br = continuation(
		re_make(prob, params = (@set par_det.q = 0.5), u0 = solhomo.u),
		PALC(), optcont; plot = true)
Scene = title!("")
Example block output
br
 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 152
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter up starts at 0.0, ends at 1.4
 ├─ Algo: PALC
 └─ Special points:

- #  1,     hopf at up ≈ +0.56970295 ∈ (+0.56970262, +0.56970295), |δp|=3e-07, [converged], δ = ( 2,  2), step =  74
- #  2,     hopf at up ≈ +0.58997693 ∈ (+0.58997627, +0.58997693), |δp|=7e-07, [converged], δ = ( 2,  2), step =  76
- #  3,     hopf at up ≈ +0.67245132 ∈ (+0.67243040, +0.67245132), |δp|=2e-05, [converged], δ = ( 2,  2), step =  84
- #  4,     hopf at up ≈ +0.86069748 ∈ (+0.86069747, +0.86069748), |δp|=1e-08, [converged], δ = (-2, -2), step = 102
- #  5,     hopf at up ≈ +1.05349842 ∈ (+1.05349825, +1.05349842), |δp|=2e-07, [converged], δ = (-2, -2), step = 120
- #  6,     hopf at up ≈ +1.15401807 ∈ (+1.15401772, +1.15401807), |δp|=3e-07, [converged], δ = (-2, -2), step = 129
- #  7, endpoint at up ≈ +1.40000000,                                                                     step = 151

We have detected 6 Hopf bifurcations. We now study the periodic orbits branching from them.

Computing the branches of Travelling waves

The periodic orbits emanating from the Hopf points look like travelling waves. This is intuitive because the equation is mostly advective as the diffusion coefficients $\nu_i$ are small. We will thus seek for travelling waves instead of periodic orbits. The advantage is that the possible Neimark-Sacker bifurcation is transformed into a regular Hopf one which allows the study of modulated travelling waves.

As we will do the same thing 3 times, we bundle the procedure in functions. We first use the regular Hopf normal form to create a guess for the travelling wave:

function getGuess(br, nb; δp = 0.005)
	nf = get_normal_form(br, nb; verbose  = false)
	pred = predictor(nf, δp)
	return pred.p, pred.orbit(0)
end

Using this guess, we can continue the travelling wave as function of a parameter. Note that in the following code, a generalized eigensolver is automatically created during the call to continuation which properly computes the stability of the wave.

function computeBranch(br, nb; δp = 0.005, max_steps = 190)
	_p, sol = getGuess(br, nb)
	# travelling wave problem
	probTW = TWProblem(
		re_make(br.prob, params = setproperties(getparams(br); up = _p)),
		getparams(br).Db,
		copy(sol),
		jacobian = :AutoDiff)
	# newton parameters with iterative eigen solver
	# eig = EigArnoldiMethod(sigma=0.2, which = BifurcationKit.LM(),x₀ = rand(2N ))
	eig = EigArpack(nev = 10, which = :LM, sigma = 0.4)
	optn = NewtonPar(verbose = true, eigsolver = eig)
	# continuation parameters
	opt_cont_br = ContinuationPar(p_min = 0.1, p_max = 1.3, newton_options = optn, ds= -0.001, dsmax = 0.01, plot_every_step = 5, detect_bifurcation = 3, nev = 10, max_steps = max_steps)
	# we build a guess for the travelling wave with speed -0.9
	twguess = vcat(sol, -0.9)
	br_wave = continuation(probTW, twguess, PALC(), opt_cont_br;
		verbosity = 3, plot = true, bothside = true,
		record_from_solution = (x, p) -> (u∞ = maximum(x[1:N]), s = x[end], amp = amplitude(x[1:N])),
		plot_solution = (x, p; k...) -> (plotsol!(x[1:end-1];k...);plot!(br,subplot=1, legend=false)),
		callback_newton = BK.cbMaxNorm(1e2),
		finalise_solution = (z, tau, step, contResult; k...) -> begin
			amplitude(z.u[N+1:2N]) > 0.01
		end,
		)
end

We can try this continuation as follows

amplitude(x) = maximum(x) - minimum(x)
br_wave = computeBranch(br, 1; max_steps = 10)
Scene = title!("")
Example block output

Building the full diagram

branches = [computeBranch(br, i) for i in 1:3]
plot(br, branches..., legend=:topleft, xlims = (0.5, 1.25), ylims=(0.5, 2.3))
Example block output

References

  • Koch

    Koch, James, Mitsuru Kurosaka, Carl Knowlen, and J. Nathan Kutz. “Multi-Scale Physics of Rotating Detonation Engines: Autosolitons and Modulational Instabilities.” ArXiv:2003.06655 [Nlin, Physics:Physics], March 14, 2020. http://arxiv.org/abs/2003.06655.