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

2

As your code is incomplete and doesn’t run, the best I can do is give some advice.

As you’re already looping over all elements, you can also query the dof-locations of a linear DG space. As you have an explicit definition of your linear function, you could simply query the value of that function at the dof location and put the value in the function array.

Regarding performing element-local L2 projections, you could check this out: Tips and Tricks — Numerical tours of continuum mechanics using FEniCS master documentation