0
$\begingroup$

I am trying to solve the following linear PDEs $$\frac{\partial u_x}{\partial x}=-[i\omega b_{||}+\nabla_\perp u_\perp],$$ $$\frac{\partial b_{||}}{\partial x}=-\frac{i}{\omega}\mathcal{L}u_x,$$ $$\mathcal{L}u_\perp = i\omega\nabla_\perp b_{||},$$ for $u_x=\texttt{ux}$, $b_{||}=\texttt{b_par}$ and $u_\perp=\texttt{u_perp}$ which are functions of $x$ and $z$. Where $$\nabla_\perp = ik_\perp - \sin\alpha\frac{\partial }{\partial z},$$ $$\mathcal{L} = \cos^2\alpha\frac{\partial^2}{\partial z^2}+\frac{\omega^2}{v_A^2(x,z)}.$$

To solve this I am trying to use the method of lines. I discretize the $z$-derivatives which gives set of ODEs. The code is for the case where $v_A(x,z)=v_{A0}$. The attached graph shows that the solutions go to infinity. Does this mean the scheme is numerically unstable? Note that if I change the line of code which reads

    dux_dx = -(1j * omega * b_par + 1j * k_perp * u_perp - \
                np.sin(alpha) / (2 * dz) * (np.roll(u_perp,-1) - np.roll(u_perp,1)))

to

    dux_dx = -(1j * omega * b_par_exact(x,z) + nabla_perp0 * u_perp_exact(x,z))

then the code works perfectly. This suggests the matrix inversion is causing the solution to go to infinity. Do you know how to fix this problem?

Code:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from scipy.integrate import solve_ivp
from scipy.linalg import solve_banded

def ux_exact(x,z):
    return ux0 * np.exp(1j * (kx * x + kz * z))

def b_par_exact(x,z):   # Note that b_par denotes b_{||}
    return b_par0 * np.exp(1j * (kx * x + kz * z))

def u_perp_exact(x,z): # Note that perp is short for perpendicular
    return u_perp0 * np.exp(1j * (kx * x + kz * z))

def calc_u_perp(x, b_par):
    # Calculate u_perp by inverting the matrix ab,
    # Note that A is nearly tridiagonal and so a computationally efficent algorithm can be used
    # See https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)
    [ab, u, v] = calc_A_u_v(x)
    d = vector_d(b_par)

    y = solve_banded((1,1), ab, d)
    q = solve_banded((1,1), ab, u)

    vTy = v[0] * y[0] + v[-1] * y[-1]
    vTq = v[0] * q[0] + v[-1] * q[-1]
    u_perp = y - vTy / (1 + vTq) * q

    return u_perp

def calc_A_u_v(x):
    a = np.full(nz, np.cos(alpha) ** 2 / dz ** 2, dtype=complex) # Off diagonal elements
    b = -2 * a + omega ** 2 / vA0 ** 2                           # Diagonal elements
    b[0]  += b[1]
    b[-1] += a[0] ** 2 / b[1]

    ab = np.zeros((3, nz), dtype=complex)
    ab[0,:] = a
    ab[1,:] = b
    ab[2,:] = a

    u = np.zeros(nz, dtype=complex)
    v = np.zeros(nz, dtype=complex)
    u[0] = -b[1]
    u[-1] = a[0]
    v[0] = 1
    v[-1] = -a[0] / b[1]

    return [ab, u, v]

def vector_d(b_par):
    return -omega * (k_perp * b_par + \
            1j * np.sin(alpha) / (2 * dz) * (np.roll(b_par,-1) - np.roll(b_par,1)))

def dUdx(x, U):
    ux = U[0:nz]
    b_par = U[nz:]
    u_perp = calc_u_perp(x,b_par)
    # dux_dx = -(1j * omega * b_par_exact(x,z) + nabla_perp0 * u_perp_exact(x,z))
    dux_dx = -(1j * omega * b_par + 1j * k_perp * u_perp - \
                np.sin(alpha) / (2 * dz) * (np.roll(u_perp,-1) - np.roll(u_perp,1)))
    db_par_dx = -1j / omega * ( \
                np.cos(alpha) ** 2 / dz ** 2 * (np.roll(ux,-1) - 2 * ux  + np.roll(ux,1))  + \
                omega ** 2 / vA0 ** 2 * ux)
    return np.concatenate((dux_dx, db_par_dx))

Lz = 1
kx = 1
kz = np.pi / Lz

alpha = 0.25 * np.pi
k_perp = 0.5 * np.pi
vA0 = 1
omega = vA0 * np.sqrt(kx ** 2 + k_perp ** 2 + kz ** 2 - 2 * kz * k_perp * np.sin(alpha) + 0j)

nx = 256
x_min = 0
x_max = 6
dx = (x_max - x_min) / (nx - 1)
x = np.linspace(x_min, x_max, nx)

nz = 128
z_min = 0
z_max =  2 * Lz
dz = (z_max - z_min) / nz
z = np.linspace(z_min + dz / 2, z_max - dz / 2, nz)

X, Z = np.meshgrid(x, z)

nabla_perp0 = 1j * (k_perp - kz * np.sin(alpha))
nabla_par0  = 1j * kz * np.cos(alpha)
L0          = nabla_par0 ** 2 + omega ** 2 / vA0 ** 2

ux0     = 1
b_par0  = -L0 / (omega * kx) * ux0
u_perp0 = -1j * kx * nabla_perp0 / (L0 + nabla_perp0 ** 2) * ux0

ux_x_min    = ux_exact(x_min, z)
b_par_x_min = b_par_exact(x_min, z)
U_x_min = np.concatenate((ux_x_min, b_par_x_min))
sol = solve_ivp(dUdx, [x_min, x_max], U_x_min, t_eval=x)

ux = sol.y[0:nz]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Z, np.real(ux))
ax.set_xlabel('x')
ax.set_ylabel('z')
ax.set_title(r'Re($u_x$)')
plt.savefig('Figures/real_part_ux.png')

Ouputted figure: enter image description here

$\endgroup$
2
  • $\begingroup$ The third equation for u_perp does not involve x derivative. How are you able to use ODE approach ? The set of equations is second order and would be a boundary value problem. What are your boundary conditions, how are you enforcing them ? $\endgroup$
    – cfdlab
    Jun 24, 2020 at 3:26
  • $\begingroup$ Yes, the third equation I think is the problem. I discretize in z to turn it into a matrix inversion problem. The BCs are periodic in z and a driver at x=x_min. I enforce them by using np.roll when taking z-derivatives and I treat the system as in IVP in x to impose the driver. $\endgroup$
    – Peanutlex
    Jun 24, 2020 at 8:54

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.