need help to implement the bound limiter (again...)

Guosheng Fu
Topic Author
Offline
Elite Member

6 years 2 months ago #761 by Guosheng Fu

need help to implement the bound limiter (again...) was created by Guosheng Fu

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....

Code:

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

christopher
Offline
Administrator

6 years 2 months ago - 6 years 2 months ago #766 by christopher

Replied by christopher on topic need help to implement the bound limiter (again...)

Hi Guosheng,
not sure if it helps, but we implemented a new feature last month: You can create numpy arrays of mapped integrationpoints and evaluate a CoefficientFunction/GridFunction in a tight C++ loop on all the points. You could then use numpy functionality to get your desired values efficiently.
I'm thinking of something like this:

Code:

# create 3 numpy arrays with x,y and z corrdinates of integration points, i.e. m points per element
xpts, ypts, zpts = np.array(...), np.array(...), np.array(...)
# create ngsolve integration points from them
mips = mesh(xpts, ypts, zpts)
# then in each step:
vals = gf(mips)
# then reshape to mxnelements numpy array and use numpy min/max on the rows

Does that help you? This would remove the need for a C++ dependency in your code...
The evaluation of cf/gf makes use of multithreading when inside a TaskManager environment.

Best
Christopher

Last edit: 6 years 2 months ago by christopher.

Guosheng Fu
Topic Author
Offline
Elite Member

6 years 2 months ago #767 by Guosheng Fu

Replied by Guosheng Fu on topic need help to implement the bound limiter (again...)

Hi Christopher,
It works much faster now!
Here is the code, maybe you can help me optimizing it

...... (the for loop over elements makes me nervous)

I am using a rectangular mesh with DG-Q1 finite element.
I need to evaluate the function at four vertices of each square element.
First, I create the integration pionts (inside the reference square, but very close to the four vertices) I used the interior points on the reference square to avoid conflicts for evaluation of L2 functions on the mesh boundaries:

Code:

eps = 1e-8
ir = IntegrationRule(points=[(0+eps,0+eps),(0+eps,1-eps), (1-eps,0+eps),(1-eps,1-eps)],                                      weights=[1/4,1/4,1/4, 1/4])      

Then, I obtain the mapped ir points of the whole mesh:

Code:

xpts = []
ypts = []
zpts = []
 for e in fes.Elements():
            trafo = e.GetTrafo()
            for p in ir :
               xpts.append(trafo(p).point[0])
               ypts.append(trafo(p).point[1])
               zpts.append(0)
mips = mesh(xpts, ypts, zpts)

Then, I define the following scaling limiter:

Code:

nd = 4   # dofs per element
npts = 4 # number of quad pts
nelems = mesh.ne
def limit(gfrho):
    rho0 = gfrho.vec[0:nelems] # cell averages
    vals = gfrho(mips).reshape(nelems, npts) # gf at quad points
    rhomin = vals.min(axis=1)
    rhomax = vals.max(axis=1)
    for i in range(nelems):
        nn = nelems +(nd-1)*i
        theta = 1
        if (rhomin[i] < r1-eps) :
            theta = (rho0[i]-r1)/(rho0[i]-rhomin[i])
        if (rhomax[i] > r2+eps):
            theta2 = (rho0[i]-r2)/(rho0[i]-rhomax[i])
            theta = min(theta, theta2)
        if theta < 1 :
            for k in range(nd-1):
                gfrho.vec[nn+k] = theta*gfrho.vec[nn+k]

Best,
Guosheng

christopher
Offline
Administrator

6 years 2 months ago - 6 years 2 months ago #768 by christopher

Replied by christopher on topic need help to implement the bound limiter (again...)

You can use numpy boolean array index assignment for that:
Something like this should work:

Code:

rho0 = gfrho.vec[0:nelems].NumPy()
theta = np.ones(nelems,dtype=float)
mask = rhomin < r1-eps
theta[mask] = (rho[mask]-r2)/(rho0[mask]-rhomin[mask])
mask = rhomax > r2+eps
theta2 = (rho0[mask] - r2)/(rho0[mask]-rhomax[mask])
theta[mask] = np.minimum(theta[mask],theta2)

The scaling is a bit tricky, but can be done using reshape and np.newaxis (so that you multiply 4 elements of the left hand side with 1 element of the array on the right hand side)

Code:

gfu.vec[nn:].NumPy().reshape(nelems,4) *= theta[:,np.newaxis]

Note that neither .NumPy() nor reshape copy the vector, but only provide views on it, so you manipulate the original GridFunction.

Best
Christopher

Last edit: 6 years 2 months ago by christopher.

Guosheng Fu
Topic Author
Offline
Elite Member

6 years 2 months ago #1171 by Guosheng Fu

Replied by Guosheng Fu on topic need help to implement the bound limiter (again...)

Dear Christopher,

Thanks for the help.
I have a fast working version now.

I noticed a bug when calling the command

Code:

mips = mesh(xpts, ypts, zpts)

It gives me negative element number, maybe there is something wrong with the mesh...

So, I used the following more direct (and more accurate) one:
(there is one slight issue needs clarification: the forth argument of mip is 'meshptr', whose number is something like 3.38634175e-316, which changes for every run...
what is this number used for? seems to be a pointer...but I am not sure)

Code:

mpts = mesh([0],[0],[0])
mpt = mpts[0][3]  # this is  the number for meshptr
ir =  [(0,0),(1,0),(1,1),(0,1)] # reference integration nodes
npts = len(ir)                                                               
mips = np.zeros((npts*nelems,), dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8'), 
             ('meshptr', '<f8'), ('VorB', '<i4'), ('nr', '<i4')])

for i in range(nelems):
       for j in range(len(ir)):                                                 
             mips[4*i+j] = (ir[j][0], ir[j][1], 0, mpt, 0, i)       

For the scaling, the reshape(...) *= do not work, so I just enlarged the theta vector to be the same dimension as the h.o. dofs:

Code:

gfrho.vec.FV().NumPy()[nelems:] *= (theta.reshape(nelems,1)*np.ones(nd-1)).ravel()    

I will buy you a drink, say at the next NGSolve user meeting... (hopefully I can come to join the party this time

)

Best,
Guosheng

christopher
Offline
Administrator

6 years 2 months ago #1175 by christopher

Replied by christopher on topic need help to implement the bound limiter (again...)

Hm, it returns -1 as element number when the point is not in the mesh (not a well documented "feature" but if you want to evaluate on a grid on a complex geometry its handy - you can just put in the grid and filter out all points which are not in the mesh by removing all with el = -1).
You are right, the 4th entry is a pointer to the mesh encoded as a double. The reason is that to be able to be used in numpy, the underlying c++ class must have a very simple structure. Amongst other things it must fullfil the requirements of StandardLayoutType . That's why we created the MeshPoint class. The BaseMappedIntegrationPoint which is then created and that is taken into the CF evaluation has a too complicated layout (it stores the trafo as a matrix, has virtual functions,...). But using the meshpointer we can construct a BaseMappedIntegrationPoint implicitly on the fly when passing to the CF.
The only pitfall there is, that the MeshPoint doesn't keep the mesh alive, so be aware to keep the Python mesh object alive as long as you use points on them.

Ah in my last line was a mistake, you can't assign to a function call... This should work and should be a little cheaper

Code:

view = gfu.vec[nn:].NumPy().reshape(nelems,4)
view *= theta[:,np.newaxis]

See you then in Vienna next year

Best
Christopher

Time to create page: 0.107 seconds