Solving a complex second order PDEs coupling with ODEs

General Usage
,
1

image
Hi, I try to solve a system including two 2nd order PDEs (as shown in the figure) coupling with six 1st ODEs.
When I set the k’’ to zero, the PDEs transform to the 1st PDEs, and i can solve it using the ROCK4 in GPU. However, once the k’’ is not zero, i try to convert the 2nd order PDEs to 1st order by du[i, 7] = u[i, 9], du[i, 8] = u[i, 10], but the solver becomes unstable. the k’’ is set to 1e-20 which is small but is the real value. The ts is a scale parameter which is set to 1 at moment. The core codes of two situation are below:

function no_gvd_kernel!(du, u, p, N, t)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i < 1 || i > N
        return 
    end

    f0, mu, hbar, T_ul, T_lg, T_ug, T_2, D,  tao, center_omega, c, n, epsilon, L, r, alpha, beta, l, J, dx, tr, wn, ntot, ts = p

    @inbounds if i >= 1 && i <= N
        if i == 1
            u[i, 7] = -r * u[i, 8]
            du[i, 8] = ts * (u[i+1, 8] - u[i, 8]) / dx
            du[i, 8] -= im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 6] * ts
            du[i, 8] += im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 8] * ts
            du[i, 8] -= (l / 2.0) * u[i, 8] * ts
            du[i, 8] *= c / n
        elseif i >= 2 && i <= N - 1
            du[i, 7] = -ts * (u[i, 7] - u[i-1, 7]) / dx
            du[i, 7] -= im * (1.0 + im * alpha) * (ntot * tao * mu *  2 * pi * f0) / (n * epsilon * c * L) * u[i, 5] * ts
            du[i, 7] += im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 7] * ts
            du[i, 7] -= (l / 2.0) * u[i, 7] * ts
            du[i, 7] *= c / n
            du[i, 8] = ts * (u[i+1, 8] - u[i, 8]) / dx
            du[i, 8] -= im * (1.0 + im * alpha) * (ntot * tao * mu *  2 * pi * f0) / (n * epsilon * c * L) * u[i, 6] * ts
            du[i, 8] += im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 8]  * ts
            du[i, 8] -= (l / 2.0) * u[i, 8] * ts
            du[i, 8] *= c / n
        elseif i == N
            u[i, 8] = -r * u[i, 7]          
            du[i, 7] = -ts * (u[i, 7] - u[i-1, 7]) / dx
            du[i, 7] -= im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 5] * ts
            du[i, 7] += im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 7] * ts
            du[i, 7] -= (l / 2.0) * u[i, 7] * ts
            du[i, 7] *= c / n
        end
        du[i, 1] = u[i, 3] / T_ul * ts - u[i, 1] / T_lg * ts - (mu / hbar) * imag(u[i, 7] * conj(u[i, 5]) + u[i, 8] * conj(u[i, 6])) * ts
        du[i, 2] = u[i, 4] / T_ul * ts - u[i, 2] / T_lg * ts + (im * mu / (2 * hbar)) * (u[i, 7] * conj(u[i, 6]) - conj(u[i, 8]) * u[i, 5]) * ts - 4.0 * wn^2 * D * u[i, 2] * ts
        du[i, 3] = ts * J - ts * (1.0 / T_ul + 1.0 / T_ug) * u[i, 3] + (mu / hbar) * imag(u[i, 7] * conj(u[i, 5]) + u[i, 8] * conj(u[i, 6])) * ts
        du[i, 4] = -ts * (1.0 / T_ul + 1.0 / T_ug) * u[i, 4] - (im * mu / (2 * hbar)) * (u[i, 7] * conj(u[i, 6]) - conj(u[i, 8]) * u[i, 5]) * ts - 4.0 * wn^2 * D * u[i, 4] * ts
        du[i, 5] = -ts * u[i, 5] / T_2 + (im * mu / (2 * hbar)) * (u[i, 7] * (u[i, 3] - u[i, 1]) + u[i, 8] * (u[i, 4] - u[i, 2])) * ts
        du[i, 6] = -ts * u[i, 6] / T_2 + (im * mu / (2 * hbar)) * (u[i, 8] * (u[i, 3] - u[i, 1]) + u[i, 7] * (conj(u[i, 4]) - conj(u[i, 2]))) * ts
    end
    return
end

function gvd_kernel!(du, u, p, N, t)
......
       if i == 1
            u[i, 7] = r * u[i, 8]
            du[i, 7] = u[i, 9] * ts
            du[i, 8] = u[i, 10] * ts
            du[i, 9] =  (u[i+1, 7] - u[i, 7]) / dx
            du[i, 9] += (n / c) * u[i, 9]
            du[i, 9] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 5]
            du[i, 9] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 7] 
            du[i, 9] += (l / 2.0) * u[i, 7]
            du[i, 9] *= 2.0 / (im * k) * ts^2
            du[i, 10] = - (-11 * u[i,8] + 18 * u[i+1,8] - 9 * u[i+2,8] + 2 * u[i+3,8]) / (6 * dx)
            #du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx
            du[i, 10] += (n / c) * u[i, 10]
            du[i, 10] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 6]
            du[i, 10] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 8]
            du[i, 10] += (l / 2.0) * u[i, 8]
            du[i, 10] *= 2.0 / (im * k) * ts^2
        elseif i == 2 || i == 3
            du[i, 7] = u[i, 9] * ts
            du[i, 8] = u[i, 10] * ts
            du[i, 9] = (u[i+1, 7] - u[i, 7]) / dx 
            du[i, 9] += (n / c) * u[i, 9]
            du[i, 9] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 5]
            du[i, 9] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 7]
            du[i, 9] += (l / 2.0) * u[i, 7]
            du[i, 9] *= 2.0 / (im * k) * ts^2
            du[i, 10] = - (-11 * u[i,8] + 18 * u[i+1,8] - 9 * u[i+2,8] + 2 * u[i+3,8]) / (6 * dx)
            #du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx
            du[i, 10] += (n / c) * u[i, 10]
            du[i, 10] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 6]
            du[i, 10] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 8]
            du[i, 10] += (l / 2.0) * u[i, 8]
            du[i, 10] *= 2.0 / (im * k) * ts^2
        elseif i >= 4 && i <= N - 3
            du[i, 7] = u[i, 9] * ts
            du[i, 8] = u[i, 10] * ts
            du[i, 9] = (11 * u[i, 7] - 18 * u[i-1, 7] + 9 * u[i-2, 7] - 2 * u[i-3, 7]) / (6 * dx)
            #du[i, 9] = (u[i, 7] - u[i-1, 7]) / dx
            du[i, 9] += (n / c) * u[i, 9]
            du[i, 9] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 5]
            du[i, 9] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 7]
            du[i, 9] += (l / 2.0) * u[i, 7]
            du[i, 9] *= 2.0 / (im * k) * ts^2 
            du[i, 10] = -(-11 * u[i,8] + 18 * u[i+1,8] - 9 * u[i+2,8] + 2 * u[i+3,8]) / (6 * dx)
            #du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx
            du[i, 10] += (n / c) * u[i, 10]
            du[i, 10] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 6]
            du[i, 10] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 8]
            du[i, 10] += (l / 2.0) * u[i, 8]
            du[i, 10] *= 2.0 / (im * k) * ts^2
        elseif i == N-1 || i == N-2
            du[i, 7] = u[i, 9] * ts
            du[i, 8] = u[i, 10] * ts
            du[i, 9] =  (11 * u[i, 7] - 18 * u[i-1, 7] + 9 * u[i-2, 7] - 2 * u[i-3, 7]) / (6 * dx)
            #du[i, 9] = (u[i, 7] - u[i-1, 7]) / dx 
            du[i, 9] += (n / c) * u[i, 9]
            du[i, 9] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 5]
            du[i, 9] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 7]
            du[i, 9] += (l / 2.0) * u[i, 7]
            du[i, 9] *= 2.0 / (im * k) * ts^2
            du[i, 10] = -(u[i, 8] - u[i-1, 8]) / dx
            du[i, 10] += (n / c) * u[i, 10]
            du[i, 10] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 6]
            du[i, 10] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 8]
            du[i, 10] += (l / 2.0) * u[i, 8]
            du[i, 10] *= 2.0 / (im * k) * ts^2
        elseif i == N
            u[i, 8] = r * u[i, 7]
            du[i, 7] = u[i, 9] * ts
            du[i, 8] = u[i, 10] * ts
            du[i, 9] = (11 * u[i, 7] - 18 * u[i-1, 7] + 9 * u[i-2, 7] - 2 * u[i-3, 7]) / (6 * dx)
            #du[i, 9] = (u[i, 7] - u[i-1, 7]) / dx 
            du[i, 9] += (n / c) * u[i, 9]
            du[i, 9] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 5]
            du[i, 9] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 7]
            du[i, 9] += (l / 2.0) * u[i, 7]
            du[i, 9] *= 2.0 / (im * k) * ts^2
            du[i, 10] = -(u[i, 8] - u[i-1, 8]) / dx
            du[i, 10] += (n / c) * u[i, 10]
            du[i, 10] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 6]
            du[i, 10] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) * u[i, 8]
            du[i, 10] += (l / 2.0) * u[i, 8]
            du[i, 10] *= 2.0 / (im * k) * ts^2
        end
 du[i, 1]=.... (same as the no_gvd_kernel!)
2

ROCK methods do not have good stability for large complex eigenvalues. Did you try something else?

3

THx, I have also try to use the AB4, but also not stable (even without k’'). In all articles I have read, the ROCK4 is used to solve the PDEs [Optics Express 24 (20), 23232(2016)]

I have done some debugging. When I increase the k’’ to 1e-10, the solver is stable but the results is wrong (the u[7] and u[8] just grow up).

I think the unstable reason of gvd_kernel! function is the set of boundary (such as i==N) u[i, 8] = r * u[i, 7]; du[i, 7] = u[i, 9]; du[i, 8] = u[i, 10]; du[i, 9] = ....; du[i, 10]=...., the value of u[8] is update by the u[i, 7] and u[i, 10], so it may make unstable. In the i == N of no_gvd_kernel function, u[i, 8] = r * u[i, 7]; du[i, 7] =... the u[8] is only update by the u[7]. So i also delete the du[i, 8] = u[i, 10] but is still unstable

4

Yes, that will lead to instability. The docs are very explicit about this. The boundary needs to not be solved for if it’s not an unknown

5

Well, the u[7] and u[8] is unknow in boundary. They are just satisfied with the reflection boundary. I have a problem about the du[i, 8] = u[i, 10] which may be conflicted with the u[i, 8] = r * u[i, 7]. Reasonable expression should be like u[i,10] = du[i, 8] in the boundary, because the update relationship is “u[9] → u[7] → u[8] → u[10]” in the really flow process.

6

They are not unknowns though. Your degrees of freedom are off. You’re solving for more variables than you have differential equations for. You either need to treat this as a DAE, or keep the boundary points as pseudo equations (i.e. the ghost point method) where you explicitly solve for the boundary values and then substitute them into the other equations and only solve the differential equations of the interior.

7

Really thanks for your reply. It’s helpful.

I’m not so familiar with the numerical solver. What i understand is that there are 10 variables corresponding to 10 differential equations. In the boundary point, one equation is replaced by the boundary condition u[i, 8] = r * u[i, 7], therefore only 9 differential equations exit in the boundary which makes the solver unstable. Is my understanding correct?

Firstly, I am curious about solving the no_gvd_kernel function process which has 8 variables (the two PDEs are 1st order). I also replace one equation by the u[i, 8] = r * u[i, 7], the ROCK4 solver is working well and the results are reasonability (but I can’t judge whether the results are correct）

The other question is about the boundary condition u[i, 8] = r * u[i, 7] which is strict. So may I differentiate the boundary conditions twice then replace it by du[i, 10] = r * du[i, 9] which can keep 10 differential equations in the boundary point?

# in the interior
du[i, 7] = u[i, 9] 
du[i, 8] = u[i, 10] 
du[i, 9] = f(u, x)
du[i,10] = g(u, x)
# in the boundary i==N
du[i, 7] = u[i, 9]
du[i, 8] = u[i, 10] 
du[i, 9] = f(u, x)
du[i, 10] = r * du[i, 9]
8

9 differential equations and one algebraic equation.

That is one way to do it. It can have numerical drift though.

The most exact way is just to eliminate u[i, 8] from the equations entirely. Any place you would use it, just use r * u[i, 7]. You don’t need to solve for it if it’s known.

9

Well, in the last two days, I have tried to do something, but looks like no useful and the solver needs larger maxiters, which results from some errors in the model i don’t find (considering the PDEs of k’'=0 being solved well).
I should explain the detail background about the two PDEs, where the E+ is propagating to the forward direction (from i=1 to i = N) and the E- is propagating to the backward direction. So the item of dE/dz is derivate to (u[i, 7] - u[i-1, 7]) / dx # u[:,7] is the E+ and (u[i+1, 8] - u[i, 8]) / dx # u[:,8] is the E-
These are what I have revised in the codes:

# in the boundary i==N,  u[i, 8] = r * u[i, 7]
   du[i, 1]=.... 
   ....
   du[i, 7] = u[i, 9]
   du[i, 8] = u[i, 10]
   du[i, 9] = f(u, x)
   du[i, 10] = g(u, x)  
# all u[i, 8] is replaced by the r * u[i, 7]

When I check the solutions, I found the u[i, 8] = r * u[i, 7] in the boundary is not achieved. I guess it that the propagating equations of du[i, 9] = f(u, x) and du[i, 10] = g(u, x) cannot result in the u[i, 8] = r * u[i, 7]. Then I do these things:

# in the boundary i==N,  u[i, 8] = r * u[i, 7]
   du[i, 1]=.... 
   ....
   du[i, 7] = u[i, 9]
   du[i, 8] = u[i, 10]
   du[i, 9] = f(u, x)
   du[i, 10] = r * du[i, 9]
# all u[i, 8] is replaced by the r * u[i, 7]

And the solver needs larger maxiters. I also consider to solve the PDEs in the interior points but I need the value of the boundary to solve the PDEs of the interior ( i == N-1, the du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx ...., the item of u[i+1, 8] is the value at boundary point). The next things I try to do is replace all u[i, 8] and u[i, 10] items by the r * u[i, 7] and r * u[i, 9] in the boundary i == N. To be honestly, I doubt it but i don’t have any ideas.

Still thanks for your reply and I am looking forward your reply. if u need the full workflows, I can send it (a little long）.

10

u[i, 8] shouldn’t be in the array at all?

11

u[i, 8]is replaced by ther * u[i, 7]at the boundary point (i = N), and for other spatial point, the points where needs the value ofu[N,8]are also employ the value of the r * u[N, 7]. Consideringdu[i, 8] = u[i, 10], I use theu[i, 10]to update theu[i, 8] at anyplace. Should I ignore the value ofu[N, 8]andu[N,10], and user * u[i, 7]explicitly to update theu[i, 8]for the boundary point after the solving PDEs?

12

it just shouldn’t be in the equations at all?

13

Sorry, I don’t know what “at all” means. Does it mean the boundary point or all spatial points? In the revised code, there is not du[i, 8], u[i, 8], and u[i, 10] in the equations at the boundary point. Now I calculate u[i, 8] after the solver having done by the u[i, 8] = r * u[i, 7] at the boundary point.

....
elseif i == N # right boundary
            du[i, 7] = u[i, 9] 
            du[i, 9] = (u[i, 7] -  u[i-1, 7]) / dx
            du[i, 9] += (n / c) * u[i, 9]
            du[i, 9] += im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u[i, 5]
            du[i, 9] -= im * beta * (abs2.(u[i, 7]) + abs2.(u[i, 8])) .* u[i, 7]
            du[i, 9] += (l / 2.0) * u[i, 7]
            du[i, 9] *= 2.0 / (im * k) 
         
            du[i, 1] = u[i, 3] / T_ul  - u[i, 1] / T_lg - (mu / hbar) * imag(u[i, 7] * conj(u[i, 5]) + r * u[i, 7] * conj(u[i, 6])) * ts
            du[i, 2] = u[i, 4] / T_ul  - u[i, 2] / T_lg  + (im * mu / (2 * hbar)) * (u[i, 7] * conj(u[i, 6]) - conj(r * u[i, 7]) * u[i, 5]) - 4.0 * wn^2 * D * u[i, 2]
            du[i, 3] = J - (1.0 / T_ul + 1.0 / T_ug) * u[i, 3] + (mu / hbar) * imag(u[i, 7] * conj(u[i, 5]) + r * u[i, 7] * conj(u[i, 6]))
            du[i, 4] = -(1.0 / T_ul + 1.0 / T_ug) * u[i, 4] - (im * mu / (2 * hbar)) * (u[i, 7] * conj(u[i, 6]) - conj(r * u[i, 7]) * u[i, 5]) - 4.0 * wn^2 * D * u[i, 4] 
            du[i, 5] =  -u[i, 5] / T_2 + (im * mu / (2 * hbar)) * (u[i, 7] * (u[i, 3] - u[i, 1]) + r * u[i, 7] * (u[i, 4] - u[i, 2]))
            du[i, 6] = -u[i, 6] / T_2 + (im * mu / (2 * hbar)) * (r * u[i, 7] * (u[i, 3] - u[i, 1]) + u[i, 7] * (conj(u[i, 4]) - conj(u[i, 2])))
14

Let’s simplify. Take a simple differential algebraic equation. I’ll just write this in ModelingToolkit.jl notation for simplicity:

D(x) ~ x + y
y ~ x + 1

Notice that PDE semi-discretizations give you differential-algebraic equations, where the interior are differential equations and the boundary conditions are algebraic equations. There are a few different ways to solve this. One is that you solve the DAE directly. This is shown in the documentation Differential Algebraic Equations · DifferentialEquations.jl. For example, here you would treat it as a mass matrix DAE:

Mu' = f(u)

where M is singular, via:

using DifferentialEquations
function massmatrixdae!(du, u, p, t)
    x,y = u
    du[1] = x + y
    du[2] = x + 1 - y
    nothing
end
M = [1.0   0.0
         0.0  0.0]
f = ODEFunction(massmatrixdae!, mass_matrix = M)
prob_mm = ODEProblem(f, [1.0, 2.0], (0.0, 1.0))
sol = solve(prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8)

That directly encodes the algebraic conditions and boom you’re off to the races. Or you could use DAEProblem, which is

f(u',u,p,t) = 0

and can be written as:

using Sundials
function f2(resid, du, u, p, t)
    dx, dy = du
    x, y = u
    resid[1] = x + y - dx
    resid[2] = x + 1 - y
end
u₀ = [1.0, 2.0]
du₀ = [-0.04, 0.04] # guess for the initial derivatives, corrected during initialization
tspan = (0.0, 1.0)
differential_vars = [true, false]
prob = DAEProblem(f2, du₀, u₀, tspan, differential_vars = differential_vars)
sol = solve(prob, IDA())

Both of those will work, but the problem with the DAE approach is that it requires the use of an implicit (or semi-implicit) method to handle the DAE constraint equations. For many systems, for example stiff equations, this is totally fine! If your equation is already stiff, yeah whatever you already are use an implicit method, this is a simple change. If your equation is non-stiff, then this can add computational overhead.

But you can then also just factor equations out. For example, if

D(x) ~ x + y
y ~ x + 1

then

D(x) ~ 2x + 1

you can just solve 1 ODE, and if you need to compute y then you just take the solution and add 1. Note that this is the kind of thing that ModelingToolkit.jl will do automatically (so I would recommend MethodOfLines.jl for this kind of thing because it’ll do all of this automatically, but wait a month or so for a few more optimizations for large PDE cases that are scheduled to arrive by the end of summer).

But you can do it by hand: that’s what I was recommending by “at all”: just completely get rid of u[i, 8] from the array, it’s not even in the equations, solve without it, and then post-process to build your extended equations.

This will be faster than the DAE approach if you want to use an explicit method. If you’re using an implicit method, it’s not worth the trouble.

There’s also three more approaches for DAE handling. One is that you can assume D(y) ~ 0 in steps and then project the algebraic conditions after steps. This is done via a callback Manifold Projection · DiffEqCallbacks.jl and can be a reasonably easy solution. It actually keeps the same order for the solution, easy proof via the triangle inequality. However, it has a downside that the continuous (dense) interpolation is ruined by the projection of the callback, so you only can trust the step values not the intermediate saveat values (so you’d need to tstops the save points). But, this can be a reasonably quick way to handle some hard DAEs.

Another approach is to psuedo-remove the algebraic equations. Instead of giving the solver:

D(x) ~ x + y
y ~ x + 1

you would instead give the solver

D(x) ~ x + y

where you inside of f you have a subroutine that solves for the algebraic equations y ~ x + 1. For simple equations this is equivalent to the elimination approach I described, but it has good scaling to other cases. It can be tedious to do by hand though, and can have issues with local minima and Newton guesses, so it generally isn’t something I’d do by hand.

Finally, you can do a singular perturbation approach. Instead of:

D(x) ~ x + y
0 ~ x + 1 - y

you’d relax it to:

D(x) ~ x + y
eps*D(y) ~ x + 1 - y

where as eps -> 0 it limits to your DAE. You flip this with mu = -1/eps

D(x) ~ x + y
D(y) ~ -mu*(x + 1 - y)

where the sign flip then ensures that the algebraic condition is the fixed point of the system, so more error pushes you back towards x + 1 - y, and then you just make mu “big enough”. This approach introduces a numerical error but can be a good way to relax some algebraic constraints in a way that makes them simpler to solve. But if mu is too large you just get a stiff ODE, in which case you should just solve the DAE directly.

So that’s the 5 ways to handle it. I’m just suggesting the elimination route as the fastest, but you can also just use the DAE solver to make it easier. Or use ModelingToolkit.jl as a tool that helps with some of the simplifications and other things. Or singular perturbation or manifold projection. Take your pick. But whatever it is, the answer doesn’t have u[i] = ... in the f.

Hopefully that full description helps.

6 Likes
15

I appreciate your reply. It’s helpful to revise my workflows. In my last reply, there is a mistake where u[i, 8] is not replaced by the r*u[i,7], and I have fixed it. du[i, 9] -= im * beta * (abs2.(u[i, 7]) + abs2.(r *u[i, 7])) .* u[i, 7]. Bad new is the solver still unstable.

Therefore, I try to use the MethodOfLines.jl. But the discretization needs too long time (In fact, I don’t run the code successfully no matter how many spatial grids is set to). And the real value of tspan is 10000*tr which is far larger than the test value. It is the reason why I avoid to use implict solver.

using DifferentialEquations, OrdinaryDiffEq, Random, CairoMakie, DelimitedFiles, Statistics, Base.Threads, HDF5, ModelingToolkit, MethodOfLines, DomainSets
const c = 3e8               # Speed of light, m/s
const e = 1.6e-19             # elementary charge
const mu = 3.8e-9 * e        # Dipole moment, m·C
const hbar = 1.05e-34  # Reduced Planck constant, unit:J·s
const T_ul = 50e-12           # Upper-level lifetime
const T_lg = 10e-12           # Lower-level lifetime
const T_ug = 100e-12          # Ground-state lifetime
const T_2 = 2e-13           # Coherence time
const D = 46e-4               # Diffusion coefficient
const f0 = 4.2e12             # Central frequency
const n = 3.6                 # Refractive index
const lamda = c / f0 / n     # Wavelength
const tao = 0.25              # Overlap parameter
const center_omega = 2 * pi * f0  # Optical frequency
const wn = n * center_omega / c      # Wavenumber
const epsilon = 8.8e-12   # Vacuum permittivity, unit: F·m-1
const Lc = 6e-3               # Cavity length
const L = 130e-9              # Period length
const r = 0.32                # Reflectivity
const k = 1e-24               # GVD fs^2/mm 
const alpha = 0.7             # Linewidth enhancement factor
const beta = 0      # Nonlinear Kerr
const l = 400.0               # Linear loss
#const ntot = 11.9e22     
const ntot = 1                             
const T_1 = T_ul * T_ug / (T_ug + T_ul)
const tr = 2 * Lc * n / c
Jth = hbar * l * epsilon * n * L *c/ (T_2 * tao * mu^2  * T_1 * pi * f0) 
pp = 1.4
J = pp * Jth
N = Int(ceil(16 * Lc / lamda))
x_grid = range(0, stop = Lc, length = N)
dx = step(x_grid)
dt =  n * dx / c 
fmax = 5e12;
dt_save =  1 / (2 * fmax)
@parameters t x
@variables u1(..) u2(..) u3(..) u4(..) u5(..) u6(..) u7(..) u8(..)
Dt = Differential(t)
Dtt = Differential(t)^2
Dx = Differential(x)
func = [Dt(u1(x, t)) ~ u3(x, t) / T_ul - u1(x, t) / T_lg - (mu / hbar) * imag(u7(x, t) * conj(u5(x, t)) + u8(x, t) * conj(u6(x, t))),
    Dt(u2(x, t)) ~ u4(x, t) / T_ul - u2(x, t) / T_lg + (im * mu / (2 * hbar)) * (u7(x, t) * conj(u6(x, t)) - conj(u8(x, t)) * u5(x, t)) - 4.0 * wn^2 * D * u2(x, t),
    Dt(u3(x, t)) ~ J - (1.0 / T_ul + 1.0 / T_ug) * u3(x, t) + (mu / hbar) * imag(u7(x, t) * conj(u5(x, t)) + u8(x, t) * conj(u6(x, t))),
    Dt(u4(x, t)) ~ -(1.0 / T_ul + 1.0 / T_ug) * u4(x, t) - (im * mu / (2 * hbar)) * (u7(x, t) * conj(u6(x, t)) - conj(u8(x, t)) * u5(x, t)) - 4.0 * wn^2 * D * u4(x, t),
    Dt(u5(x, t)) ~ -u5(x, t) / T_2 + (im * mu / (2 * hbar)) * (u7(x, t) * (u3(x, t) - u1(x, t)) + u8(x, t) * (u4(x, t) - u2(x, t))),
    Dt(u6(x, t)) ~ -u6(x, t) / T_2 + (im * mu / (2 * hbar)) * (u8(x, t) * (u3(x, t) - u1(x, t)) + u7(x, t) * (conj(u4(x, t)) - conj(u2(x, t)))),
    Dtt(u7(x, t)) ~ (n / c * Dt(u7(x, t)) + Dx(u7(x, t)) + im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u5(x, t) - im * beta * (abs2.(u7(x, t)) + abs2.(u8(x, t))) * u7(x, t) + (l / 2.0) * u7(x, t)) / (im * k),
    Dtt(u8(x, t)) ~ (n / c * Dt(u8(x, t)) - Dx(u8(x, t)) + im * (1.0 + im * alpha) * (ntot * tao * mu * 2 * pi * f0) / (n * epsilon * c * L) * u6(x, t) - im * beta * (abs2.(u7(x, t)) + abs2.(u8(x, t))) * u8(x, t) + (l / 2.0) * u8(x, t)) / (im * k)]
Random.seed!(42)
bcs = [u1(x, 0) ~ abs.(ComplexF64.(randn(), randn())),
    u2(x, 0) ~ ComplexF64.(abs.(randn()), randn()),
    u3(x, 0) ~ abs.(ComplexF64.(randn(), randn())),
    u4(x, 0) ~ ComplexF64.(abs.(randn()), randn()),
    u5(x, 0) ~ abs.(ComplexF64.(randn(), randn())),
    u6(x, 0) ~ ComplexF64.(abs.(randn()), randn()),
    u7(x, 0) ~ ComplexF64.(abs.(randn()), randn()),
    u8(x, 0) ~ ComplexF64.(abs.(randn()), randn()),
    0 ~ r * u8(0, t) - u7(0, t), 
    0 ~ r * u7(Lc, t) - u8(Lc, t)]
M = 2
T_patch = 1
domains = [x ∈ Interval(0, Lc),
         t ∈ Interval(0, M * T_patch * tr)]
@named pdesys = PDESystem(func, bcs, domains, [x, t], [u1(x, t), u2(x, t), u3(x, t), u4(x, t), u5(x, t), u6(x, t), u7(x, t), u8(x, t)])
discretization = MOLFiniteDifference([x => dx], t)
prob = discretize(pdesys, discretization)
sol = solve(prob, ROCK4(), dt=dt, saveat=dt_save, dense=false, save_everystep=false, abstol=1e-4, reltol=1e-4)
16

What does the instability look like? And using what solver?

Yeah, MOL still needs a codegen fix. I think I mentioned above, I’ll get to that this summer.

17

I used the ROCK4 solver. It works well for the 1st order PDEs (k’’ =0) but needs large maxiters in 2nd order PDEs (k’’ =! 0). I think it results from the unstable. For the long tspan, I divide the tspan into M patches and save the results independently. The last results of M-1 patches will be as the initial value of the next patches. When the solver runs after some patches (just two or three patches), there is a warning “The solver needs large maxiters”. I have checked the saved results, and compared it with the right calculated results (k’’ =0)

In the case of solver unstable, the results are significantly larger than the right results and keep growing up (I’m not sure if it will keep increasing because the solver is no longer working). The true value should be oscillated between positive and negative.

18

ROCK4 has issues with large complex eigenvalues. Did you verify the stability with a simple method like Tsit5 first?

This suggests an issue with the boundary condition implementation. How did you do the upwinding?

19

I have used Tsit5, which has same trouble with ROCK4. I checked the results and found it that the u7 or u8 would grow up to the level of 1e7, then the solver would be unstable and needs large maxiters. The right value of u7 or u8 should be the level of 1e5.

Here is the manual discretization, I used the toward/backward difference scheme to achieve the spatial discretization. Thx for your help

@inbounds if i >= 1 && i <= N
        if i == 1        #u[i, 7] = r * u[i, 8]
           du[i, 8] = u[i, 10]
           du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx
           du[i, 10] +=.....
           du[i, 1] +=.....
           du[i, 2] +=.....
           ..... # no u[i, 7] neither u[i,9]
        elseif i == 2
           du[i, 7] = u[i, 9]
           du[i, 8] = u[i, 10]
           du[i, 9] =  (u[i, 7] - r * u[i-1, 8]) / dx
           du[i, 9] +=.....
           du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx
           du[i, 10] +=.....
           .....
        elseif i >= 3 && i <= N - 2
           du[i, 7] = u[i, 9]
           du[i, 8] = u[i, 10]
           du[i, 9] =  (u[i, 7] - u[i-1, 7]) / dx
           du[i, 9] +=.....
           du[i, 10] = -(u[i+1, 8] - u[i, 8]) / dx
           du[i, 10] +=.....
           .....
        elseif i == N-1
           du[i, 7] = u[i, 9]
           du[i, 8] = u[i, 10]
           du[i, 9] =  (u[i, 7] - u[i-1, 7]) / dx
           du[i, 9] +=.....
           du[i, 10] = -(r * u[i+1, 7] - u[i, 8]) / dx
           du[i, 10] +=.....
           .....
        elseif i == N     #u[i, 8] = r * u[i, 7]
           du[i, 7] = u[i, 9]
           du[i, 9] =  (u[i, 7] - u[i-1, 7]) / dx
           du[i, 9] +=.....
           ..... # no u[i, 8] neither u[i,10]
       end
end
20

I’m not sure that’s going to be stable. You probably need to just wind the dE/dz term.