How to implement a slope limiter?

General
1

Hi,

I am implementing a slope limiter, for example, presented in “Python framework for hp-adaptive discontinuous Galerkin methods for two-phase flow in porous media”.

Basically, how to implement the limiter is that firstly we access each element and get solution values at quadrature point in each element, secondly calculate the slope factor according to the theory of slope limiter and define a linear function, and lastly compute L^2 projection to get DOF values and substitute the new DOF values to the original solution vector.

I have figured out how to do the first and second, but not the last. Could you please teach me how to compute L^2 projection of such a define linear function in each element and replace the new DOF values?

In my code below, sh is the DG solution so contains DOF values. deg and domain are given.

points, _ = basix.quadrature.make_quadrature(
    CellType.triangle, deg*8, QuadratureType.xiao_gimbutas
)

num_cells = domain.topology.index_map(domain.topology.dim).size_local

cell_indices = np.arange(num_cells, dtype=np.int32)
cell_values  = np.arange(1, num_cells+1, dtype=np.int32)

cell_tag = mesh.meshtags(domain, domain.topology.dim, cell_indices, cell_values)

dx = ufl.Measure(
    "dx",
    domain=domain,
    subdomain_data=cell_tag,
    metadata={
        "quadrature_degree": deg*2,
        "quadrature_rule": "xiao_gimbutas",
    },
)

one = fem.Constant(domain, 1.0)

for cell_idx in range(num_cells):
    cell_dofs = sh.function_space.dofmap.cell_dofs(cell_idx)

    vol = fem.assemble_scalar(fem.form(one * dx(cell_idx+1)))
    avg = fem.assemble_scalar(fem.form(sh * dx(cell_idx+1))) / vol

    cell_coords = domain.geometry.x[
        domain.topology.connectivity(domain.topology.dim, 0).links(cell_idx)
    ]

    qp_phys = domain.geometry.cmap.push_forward(points, cell_coords)

    s_vals = np.zeros((points.shape[0], 1))
    sh.eval(qp_phys,
            np.full(points.shape[0], cell_idx, dtype=np.int32),
            s_vals)

    chi = []
    for sa in s_vals[:, 0]:
        denom = avg - sa
        if abs(denom) < 1e-12:
            chi.append(1.0)
        else:
            a = abs((avg - 0.0) / denom)
            b = abs((1.0 - avg) / denom)
            chi.append(min(1.0, a, b))

    theta = min(chi)
# --- Apply limiter to DOFs ---
        #need to define a linear function like \bar{s}(x) = theta*(s(x)-avg) + avg, x is quadrature point
#need to compute L^2 projection of  \bar{s}(x) in each element 
#need to substitute the new DOF values to the right places in sh

Thank you