Transpose needed ? within backsolve for Crank Nicolson when integrating advection/diffusion equation

[aquillen]

I have been trying to simulate in skfem the advection diffusion equation

$$ \partial_t u +( {\bf a} \cdot \nabla ) u = D_u \Delta u $$

on the unit disk with rotation $${\bf a} = \omega (x \hat {\bf y} - y \hat {\bf x}) $$ and $\omega$ an angular rotation rate.

I have been following examples 25,28, 42, 50 for advection (but which use solve) and examples 19, 39 which use the Crank Nicolson method and pythons sparse matrix splu/backsolve for time dependent problems.

I could not get advection to take place unless there is a .T on the argument within the backsolve call. The examples 19, 39 which work, didn't have a .T within the backsolve call. Perhaps the laplacian operator is symmetric and the .T is not needed in those cases? I thought I would point this out just in case somebody else is confused by this issue. I am not 100% sure on where the inconsistency is. Is it with the splu call, the definition of the bilinear form for advection or with the backsolve call?

Below is the code I have been running:

from scipy.sparse.linalg import splu  # for sparse matrices, 
# returns something like invA with LU decomp which can be applied with solve to solve Ax = b for x

from skfem.models.poisson import laplace, mass

mcirc= skfem.MeshTri.init_circle(5,smoothed=False)  # make a triangular mesh of the unit circle

element = skfem.ElementTriP1()  #  a triangular element
basis = skfem.Basis(mcirc, element)

dt = .005   # timestep 
D_u = 1e-2  # diffusion coefficent 
om_circ = 1.0;    # angular rotation rate for advection 

@skfem.BilinearForm   # for advection 
def advection(u, v, w):
    from skfem.helpers import grad
    return v*(w.x[1]*grad(u)[0] - w.x[0]*grad(u)[1])  

L_advec = skfem.asm(advection, basis) 
L_lap = skfem.asm(laplace, basis) 
M = skfem.asm(mass, basis)  # mass matrix (sizes of triangles) 

L = D_u * L_lap + om_circ * L_advec
theta = 0.5                     # Crank–Nicolson algorithm which is meant for diffusive equations not advection 
A = M + theta * L * dt
B = M - (1 - theta) * L * dt  
backsolve = splu(A.T).solve  # .T as splu prefers CSC which is a sparse matrix format (compressed sparse column)

# do a single timestep, return u field and time 
# requires globals:  dt, backsolve, B
def one_step(t,u):
    t += dt 
    new_u = backsolve(B.T @ u )  # the advection does not work unless .T is on B 
    # -- It took me a while to figure this out!
    # the use of .T here conflicts with example ex19.py
    return t,new_u

# display single field, calls skfem's plot routine 
def disp(mesh,u,basis,t):
    fig,ax = plt.subplots(1,1,figsize=(3,3),sharex=True,sharey=True)
    ax.set_aspect('equal');
    ax_j0 = plot(mesh, u[basis.nodal_dofs.flatten()], shading='gouraud', ax = ax)
    field0 = ax.get_children()[0]  # vertex-based temperature-colour
    fig.colorbar(field0,shrink=0.7)
    title = ax.set_title(f'$t$ = {t:.3f}')

# initial conditions 
x_p = basis.doflocs[0]; y_p = basis.doflocs[1];  t = 0 
x_0 =0.5; y_0= 0.0; sig = 0.1
r2 = (x_p - x_0)**2 + (y_p - y_0)**2 ;
u_init = np.exp(-0.5*r2/sig**2)   # make u have an offset gaussian peak 
disp(mcirc,u_init,basis,t)

# run a short integration and show results
u = u_init;  t=0; nj = int(1/dt)+1
disp(mcirc,u,basis,t)
for j in range(4):
    for i in range(nj):
        t, u= one_step(t,u)  
    disp(mcirc,u,basis,t)
disp(mcirc,u-u_init,basis,t)