# Source code for firedrake.slope_limiter.vertex_based_limiter

```from firedrake import dx, assemble, LinearSolver
from firedrake.function import Function
from firedrake.cofunction import Cofunction
from firedrake.functionspace import FunctionSpace
from firedrake.parloops import par_loop, READ, RW, MIN, MAX
from firedrake.ufl_expr import TrialFunction, TestFunction
from firedrake.slope_limiter.limiter import Limiter
from firedrake import utils
from ufl import inner
__all__ = ("VertexBasedLimiter",)

[docs]
class VertexBasedLimiter(Limiter):
"""
A vertex based limiter for P1DG fields.

This limiter implements the vertex-based limiting scheme described in
Dmitri Kuzmin, "A vertex-based hierarchical slope limiter for p-adaptive
discontinuous Galerkin methods". J. Comp. Appl. Maths (2010)
http://dx.doi.org/10.1016/j.cam.2009.05.028
"""

def __init__(self, space):
"""
Initialise limiter

:param space : FunctionSpace instance
"""

if utils.complex_mode:
raise ValueError("We haven't decided what limiting complex valued fields means. Please get in touch if you have need.")

self.P1DG = space
self.P1CG = FunctionSpace(self.P1DG.mesh(), 'CG', 1)  # for min/max limits
self.P0 = FunctionSpace(self.P1DG.mesh(), 'DG', 0)  # for centroids

# Storage containers for cell means, max and mins
self.centroids = Function(self.P0)
self.centroids_rhs = Cofunction(self.P0.dual())
self.max_field = Function(self.P1CG)
self.min_field = Function(self.P1CG)

self.centroid_solver = self._construct_centroid_solver()

# Update min and max loop
domain = "{[i]: 0 <= i < maxq.dofs}"
instructions = """
for i
maxq[i] = fmax(maxq[i], q[0])
minq[i] = fmin(minq[i], q[0])
end
"""
self._min_max_loop = (domain, instructions)

# Perform limiting loop
domain = "{[i, ii]: 0 <= i < q.dofs and 0 <= ii < q.dofs}"
instructions = """
<float64> alpha = 1
<float64> qavg = qbar[0, 0]
for i
<float64> _alpha1 = fmin(alpha, fmin(1, (qmax[i] - qavg)/(q[i] - qavg)))
<float64> _alpha2 = fmin(alpha, fmin(1, (qavg - qmin[i])/(qavg - q[i])))
alpha = _alpha1 if q[i] > qavg else (_alpha2 if q[i] < qavg else  alpha)
end
for ii
q[ii] = qavg + alpha * (q[ii] - qavg)
end
"""
self._limit_kernel = (domain, instructions)

def _construct_centroid_solver(self):
"""
Constructs a linear problem for computing the centroids

:return: LinearSolver instance
"""
u = TrialFunction(self.P0)
v = TestFunction(self.P0)
a = assemble(inner(u, v) * dx)
return LinearSolver(a, solver_parameters={'ksp_type': 'preonly',
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu'})

def _update_centroids(self, field):
"""
Update centroid values
"""
assemble(inner(field, TestFunction(self.P0)) * dx, tensor=self.centroids_rhs)
self.centroid_solver.solve(self.centroids, self.centroids_rhs)

[docs]
def compute_bounds(self, field):
"""
Only computes min and max bounds of neighbouring cells
"""
self._update_centroids(field)
self.max_field.assign(-1.0e10)  # small number
self.min_field.assign(1.0e10)  # big number

par_loop(self._min_max_loop,
dx,
{"maxq": (self.max_field, MAX),
"minq": (self.min_field, MIN),

[docs]
def apply_limiter(self, field):
"""
Only applies limiting loop on the given field
"""
par_loop(self._limit_kernel, dx,
"q": (field, RW),

[docs]
def apply(self, field):
"""
Re-computes centroids and applies limiter to given field
"""
assert field.function_space() == self.P1DG, \
'Given field does not belong to this objects function space'

self.compute_bounds(field)
self.apply_limiter(field)

```