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]