Mixed Function Space Problem

I try to solve the heat equation

$https://latex.codecogs.com/svg.image?\left\{\begin{matrix} u_t - \Delta u = 0, \quad \Omega = [0,2\pi]\times[0,2\pi] \\u(x,y,0) = \sin(x)\sin(y) \ \end{matrix}\right.$

with periodic boundary conditions . The exact solution is $https://latex.codecogs.com/svg.image?u(x,y,t) = e^{-2t} \sin(x)\sin(y)$

Local Discontinuous Galerkin(LDG) method is used for spatial discretization, we first introduce an auxiliary vector-valued variable $https://latex.codecogs.com/svg.image?\boldsymbol{p}$

then find $https://latex.codecogs.com/svg.image?u \in V, \boldsymbol{p} \in \Sigma, \forall v \in V, \boldsymbol{w} \in \Sigma$ , s.t.

$https://latex.codecogs.com/svg.image?\begin{cases}\int_{\Omega}u_t v dx &= -\int_{\Omega}\boldsymbol{p} \cdot \nabla v dx + \int_{\Gamma_{int}} \hat{\boldsymbol{p}} \cdot [[v]]_{\boldsymbol{n}} dS,\\ \int_{\Omega} {\boldsymbol{p}} \cdot\boldsymbol{w} dx &= -\int_{\Omega} u \nabla \cdot \boldsymbol{w} dx +\int_{\Gamma_{int}} \hat{u} [[\boldsymbol{w}]]_{\boldsymbol{n}} dS , \end{cases}$

where the flux is

$https://latex.codecogs.com/svg.image?\hat{\boldsymbol{p}} = \boldsymbol{p}^{+}, \hat{u} = u^{-}.$

Two-stage, second-order DIRK method is used for time discretization, to discretise $u_t = L(u),$ we set

$https://latex.codecogs.com/svg.image?\begin{align*} u^{(1)} &= u^{n} + \Delta t a_{11} L(u^{(1)}), \nonumber \\u^{(2)} &= u^{n} + \Delta t \left( a_{21} L(u^{(1)}) + a_{22} L(u^{(2)}) \right) , \\u^{n+1} &= u^{n} + \Delta t \left( b_1 L(u^{(1)}) + b_{2} L(u^{(2)}) \nonumber \right),\end{align*}$

where $https://latex.codecogs.com/svg.image?a_{11} = a_{22} = \alpha, a_{12} = 1- \alpha, b_{1} = 1-\alpha, b_{2} = \alpha, \alpha = 1-\frac{\sqrt{2}}{2}$

I follow the example of Mixed formulation for the Poisson , but there are some problems. Here is my code:

from firedrake import *
import math

Nx = 10; Ny = 10; 
Pk = 1
CFL = 0.1

Lx = 2*math.pi; Ly = 2*math.pi
Hx = Lx/Nx
mesh = PeriodicRectangleMesh(Nx, Ny, Lx, Ly, direction='both',quadrilateral=False, reorder=None, diagonal="right")
x, y = SpatialCoordinate(mesh)
n = FacetNormal(mesh)

V = FunctionSpace(mesh, "DG", Pk)
Simga = VectorFunctionSpace(mesh, "DG", Pk)
W = V * Simga

T = 0.001
dt = CFL * Hx
dtc = Constant(dt)

u_init = Function(V).interpolate(sin(x) * sin(y))
u = Function(W.sub(0))
u.assign(u_init)

# coefficients of SDIRK 
alpha = 1.0 - math.sqrt(2.0)/2.0
RK_a = Constant(alpha); RK_a21 = Constant(1.0 - alpha)
RK_b1 = Constant(1.0 - alpha); RK_b2 = Constant(alpha)

u_trial, p_trial = TrialFunctions(W)
v, w = TestFunctions(W)

LHS = u_trial * v *dx \
    - dtc*RK_a * ( -dot(p_trial,grad(v))*dx +  dot(p_trial('+'),jump( v,n))*dS ) \
    + dot(p_trial,w)*dx \
    - (-u_trial*div(w)*dx + u_trial('-')*jump(w,n)*dS)  
A = assemble(LHS)                                                                           
M = assemble( u_trial*v *dx)

p1 = Function(W.sub(1))
p2 = Function(W.sub(1))

L1 = u*v *dx
L2 =  u*v *dx + dtc * RK_a21 *(-dot(p1,grad(v))*dx +  dot( p1('+'),jump(v,n) )*dS)          
L = dtc * (RK_b1 * (-dot(p1,grad(v))*dx +  dot( p1('+'), jump(v,n) )*dS) \
        + RK_b2 * (-dot(p2,grad(v))*dx +  dot( p2('+'), jump(v,n) )*dS) )                  

w_solution = Function(W)

b = assemble(L1)
solve(A, w_solution, b)
p1.assign(w_solution.sub(1))

b = assemble(L2)
solve(A, w_solution, b)
p2.assign(w_solution.sub(1))

b = assemble(L)
solve(M, w_solution, b)
u.assign(u + w_solution.sub(0))

u_exact = pow(e, -2*T) * sin(x) * sin(y)
L2_error = norms.errornorm( u_exact, u, norm_type = "L2")

print("Nx = ",Nx ,", L2_error = ", L2_error)

Error:

Traceback (most recent call last):
  File "/home/mycode/Code/firedrake_code/Heat_LDG.py", line 71, in <module>
    solve(M, w_solution, b)
  File "PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/mycode/firedrake/src/firedrake/firedrake/adjoint/solving.py", line 53, in wrapper
    output = solve(*args, **kwargs)
  File "/home/mycode/firedrake/src/firedrake/firedrake/solving.py", line 133, in solve
    return _la_solve(*args, **kwargs)
  File "/home/mycode/firedrake/src/firedrake/firedrake/solving.py", line 244, in _la_solve
    solver.solve(x, b)
  File "PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/mycode/firedrake/src/firedrake/firedrake/linear_solver.py", line 162, in solve
    raise ConvergenceError("LinearSolver failed to converge after %d iterations with reason: %s", self.ksp.getIterationNumber(), solving_utils.KSPReasons[r])
firedrake.exceptions.ConvergenceError: ('LinearSolver failed to converge after %d iterations with reason: %s', 0, 'DIVERGED_PCSETUP_FAILED')

firedrakeproject / firedrake

Mixed Function Space Problem #2389