Dear all,
I need a bound-preserving limiter to work for scalar DG. (I am posting the same question again… )
The problem is that I have a DG-Pk grid function, whose cell averages lie in the bound [r1, r2], but the variations inside each element violate this bound.
So, for each element, I need to access the grid function at certain quad points, and find the min/max values at these quad points. Then, I apply a scaling limiter so that the min/max values at these quad points stay within bound.
This is cheap to perform, if implemented correctly, and is crucial for the robustness of DG compressible flow simulations…
Here I have a plain python implementation of the limiter, but way tooooo slow…
ir = IntegrationRule(TRIG, 3)
nd = (order+1)*(order+2)/2
DG space: all-dofs-together...
def limit(gfrho):
for i, e in enumerate(fes_rho.Elements()):
trafo = e.GetTrafo()
id0 = int(nd*i)
rho0 = gfrho.vec[id0]
val = []
for p in ir :
val.append(gfrho(trafo(p)))
rhoS = []
rhoSum = 0.
for k in range(3*(order+1)):
tmp = val[k] # density at quad pt
rhoS.append(tmp)
rhoSum += ir.weights[k]*tmp
if order > 1:
rho3 = (rho0-rhoSum)/(1-2*mw)
rhoS.append(rho3)
rhomin = min(rhoS)
rhomax = max(rhoS)
theta = 1
if (rhomin < r1-eps) :
theta = (rho0-r1)/(rho0-rhomin)
if (rhomax > r2+eps):
theta1 = (rho0-r2)/(rho0-rhomax)
theta = min(theta, theta1)
if theta < 1 :
# this is too aggressive, FIXME later
for k in range(nd-1):
gfrho.vec[id0+k+1] = theta*gfrho.vec[id0+k+1] # the update
I realized that there is no way for me to put this piece of code to C++ all by myself, so I am turning to the experts for help… (for sure I will buy you drink if we meet at some imaginary conference later :>)
Best,
Guosheng