Sparse Direct Solvers cannot work for HDG

[zhyxue0]

Hello, developers. I tested the example in https://github.com/firedrakeproject/firedrake/blob/master/tests/slate/test_hdg_poisson.py, but I changed the solver to a direct solver(mumps) and got the wrong results.

Here is my code：

from firedrake import *
def Norm_func(u, p, quad_deg):
    return assemble(abs(u)**p*dx(degree=quad_deg))**(1/p)

def Norm_Vec_func(u1, u2, p, quad_deg):
    return assemble(  abs(u1)**p*dx(degree=quad_deg)
                    + abs(u2)**p*dx(degree=quad_deg) )**(1/p)

def run_LDG_H_problem(r, degree, quads=False):
    """
    Solves the Dirichlet problem for the elliptic equation:
    -div(grad(u)) = f in [0, 1]^2, u = g on the domain boundary.
    The source function f and g are chosen such that the analytic
    solution is:
    u(x, y) = sin(x*pi)*sin(y*pi).
    This test uses an HDG discretization (specifically LDG-H by Cockburn)
    and Slate to perform the static condensation and local recovery.
    """
    # Set up problem domain
    # mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=quads)
    Nx = 2**r
    Ny = Nx
    Lx = 2*pi
    Ly = Lx
    Hx = Lx / Nx
    print( "degree = ", degree, ",r = ", r, ",Nx = ",Nx)

    mesh = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=Qk_flag, reorder=None, diagonal="left")
    x, y = SpatialCoordinate(mesh)
    n = FacetNormal(mesh)

    # Set up function spaces
    U = VectorFunctionSpace(mesh, "DG", degree)
    V = FunctionSpace(mesh, "DG", degree)
    T = FunctionSpace(mesh, "HDiv Trace", degree)

    # Mixed space and test/trial functions
    W = U * V * T
    s = Function(W).assign(0.0)
    # q, u, uhat = split(s)
    q, u, uhat = TrialFunctions(W)
    v, w, mu = TestFunctions(W)

    # Analytical solutions for u and q
    V_a = FunctionSpace(mesh, "DG", degree + 3)
    U_a = VectorFunctionSpace(mesh, "DG", degree + 3)

    u_a = Function(V_a)
    # a_scalar = sin(pi*x[0])*sin(pi*x[1])
    # q_exact_0 = pi*cos(pi*x[0])*sin(pi*x[1])
    # q_exact_1 = pi*sin(pi*x[0])*cos(pi*x[1])
    a_scalar = sin(x)*sin(y)
    q_exact_0 = -cos(x)*sin(y)
    q_exact_1 = -sin(x)*cos(y)
    u_a.interpolate(a_scalar)

    q_a = Function(U_a)
    a_flux = -grad(a_scalar)
    q_a.project(a_flux)

    Vh = FunctionSpace(mesh, "DG", degree + 3)
    f = Function(Vh).interpolate(-div(grad(a_scalar)))

    # Stability parameter
    tau = Constant(1.0)

    # Numerical flux
    qhat = q + tau*(u - uhat)*n

    # Formulate the LDG-H method in UFL
    a = ((inner(q, v) - inner(u, div(v)))*dx
         + inner(uhat('+'), jump(v, n=n))*dS
         + inner(uhat, dot(v, n))*ds
         - inner(grad(w), q)*dx
         + inner(jump(qhat, n=n), w('+'))*dS
         + inner(dot(qhat, n), w)*ds
         # Transmission condition
         + inner(jump(qhat, n=n), mu('+'))*dS
         # BC condtion
         + inner(uhat, mu)*ds
         )

    L = inner(f, w)*dx

    solver_parameters_fd = {'ksp_type': 'preonly',
              'mat_type': 'matfree',
              'pmat_type': 'matfree',
              'ksp_type': 'preonly',
              'pc_type': 'python',
              'pc_python_type': 'firedrake.SCPC',
              'pc_sc_eliminate_fields': '0, 1',
              'condensed_field': {'ksp_type': 'preonly',
                                  'pc_type': 'redundant',
                                  "redundant_pc_type": "lu",
                                  "redundant_pc_factor_mat_solver_type": "mumps",
                                  "redundant_mat_mumps_icntl_14": 200},
                'ksp_converged_reason': None,
                'ksp_monitor_true_residual': None,}

    solver_parameters_lu ={
                "ksp_type": "preonly",
                "pc_type": "lu",
                "pc_factor_mat_solver_type": "mumps",
                # "pc_factor_mat_solver_type": "umfpack",
                # "pc_factor_mat_solver_type": "superlu_dist",
                'ksp_converged_reason': None,
                'ksp_monitor_true_residual': None,
                # 'ksp_view': None
                }

    problem = LinearVariationalProblem(a, L, s)
    # solver = LinearVariationalSolver(problem, solver_parameters=solver_parameters_fd)
    solver = LinearVariationalSolver(problem, solver_parameters=solver_parameters_lu)
    solver.solve()

    # Computed flux, scalar, and trace
    q_h, u_h, uhat_h = s.split()

    scalar_error = Norm_func(a_scalar - u_h, quad_deg, 2)
    flux_error = Norm_Vec_func(q_exact_0 - q_h.sub(0),
                               q_exact_1 - q_h.sub(1),
                               quad_deg, 2)
    return scalar_error, flux_error

def test_hdg_convergence(degree, quads):
    import numpy as np
    diff = np.array([run_LDG_H_problem(r, degree, quads) for r in range(1, 7)])
    conv = np.log2(diff[:-1] / diff[1:])
    # assert (np.array(conv) > rate).all()
    print("Error:\n", diff)
    print("Order: \n", conv)

quad_deg = 10
Qk_flag = False

test_hdg_convergence(1, True)

The outputs:

degree =  1 ,r =  1 ,Nx =  2
    Residual norms for firedrake_0_ solve.
    0 KSP none resid norm 6.655911960500e+00 true resid norm 3.048491655598e+01 ||r(i)||/||b|| 4.580126170072e+00
    1 KSP none resid norm 3.048491655598e+01 true resid norm 3.048491655598e+01 ||r(i)||/||b|| 4.580126170072e+00
  Linear firedrake_0_ solve converged due to CONVERGED_ITS iterations 1
degree =  1 ,r =  2 ,Nx =  4
    Residual norms for firedrake_1_ solve.
    0 KSP none resid norm 3.644057766782e+00 true resid norm 2.197698792159e+04 ||r(i)||/||b|| 6.030910959186e+03
    1 KSP none resid norm 2.197698792159e+04 true resid norm 2.197698792159e+04 ||r(i)||/||b|| 6.030910959186e+03
  Linear firedrake_1_ solve converged due to CONVERGED_ITS iterations 1
degree =  1 ,r =  3 ,Nx =  8
    Residual norms for firedrake_2_ solve.
    0 KSP none resid norm 1.963770418814e+00 true resid norm 7.553133486707e+06 ||r(i)||/||b|| 3.846240586143e+06
    1 KSP none resid norm 7.553133486707e+06 true resid norm 7.553133486707e+06 ||r(i)||/||b|| 3.846240586143e+06
  Linear firedrake_2_ solve converged due to CONVERGED_ITS iterations 1
degree =  1 ,r =  4 ,Nx =  16
    Residual norms for firedrake_3_ solve.
    0 KSP none resid norm 1.000868840903e+00 true resid norm 4.743938337282e+06 ||r(i)||/||b|| 4.739820187631e+06
    1 KSP none resid norm 4.743938337282e+06 true resid norm 4.743938337282e+06 ||r(i)||/||b|| 4.739820187631e+06
  Linear firedrake_3_ solve converged due to CONVERGED_ITS iterations 1
degree =  1 ,r =  5 ,Nx =  32
    Residual norms for firedrake_4_ solve.
    0 KSP none resid norm 5.028479863453e-01 true resid norm 1.736002421009e+08 ||r(i)||/||b|| 3.452340405350e+08
    1 KSP none resid norm 1.736002421009e+08 true resid norm 1.736002421009e+08 ||r(i)||/||b|| 3.452340405350e+08
  Linear firedrake_4_ solve converged due to CONVERGED_ITS iterations 1
degree =  1 ,r =  6 ,Nx =  64
    Residual norms for firedrake_5_ solve.
    0 KSP none resid norm 2.517269657441e-01 true resid norm 3.308018183271e+09 ||r(i)||/||b|| 1.314129447154e+10
    1 KSP none resid norm 3.308018183271e+09 true resid norm 3.308018183271e+09 ||r(i)||/||b|| 1.314129447154e+10
  Linear firedrake_5_ solve converged due to CONVERGED_ITS iterations 1
Error:
 [[1.32162006e+16 9.83371216e+15]
 [9.84163661e+17 7.45248115e+17]
 [6.32317771e+21 8.78838065e+21]
 [1.99096321e+21 1.52297812e+21]
 [4.63012836e+22 7.16914976e+22]
 [2.71477715e+22 2.82819853e+22]]
Order: 
 [[ -6.21851886  -6.24384088]
 [-12.6494639  -13.52558892]
 [  1.6671832    2.52870215]
 [ -4.53951372  -5.55683491]
 [  0.77021841   1.34191863]]

[wence-]

The hybridised system is a big saddle-point system, so you will probably have to hit mumps quite hard to do the pivoting correctly to invert it (there's a big zero block on the diagonal for the trace variables). I forget which icntl_ things handle that, if it all.