[dune-fem] Problem with flux limiters

[Torben.Prill at dlr.de]

Hi Everyone,



my name is Torben Prill and I work for the German Aerospace Center (DLR) at their site in Stuttgart. Mostly I work on reactive flow/transport problems in porous media.

I have been using the dune-fem python bindings since last year in two of my projects with great success, so thank you for providing such a versatile and ,still, performant library.



Right now, I am using dune-fem/-fem-dg for developing a simulator for granular flow. More specifically, I am using dune-fem and dune-fem-dg at version 2.8.0.dev20210318 from pip.

The problem consists of a convection equation and an equation for the convection velocity, which is computed from an equation using a strain rate dependent stress tensor:

d \rho /dt = \nabla \rho(x,t) \cdot v(x,t)

0 = \sigma ( \epsilon(v), \rho) (x,t)



The spatial discretization of \rho is using first order DG, the mechanical problem is solved with first order FE.



As \rho must be in certain limits to be physical (say \rho \in [0,1]), I was looking at using the scaling limiter. However, the limiter seems to have no effect when applied to a function. I attached a small example where a function undershoots into the negatives, which should be "scaled away", but isn't.



Some help (or an explanation if I am getting something wrong) would be greatly appreciated.



Many thanks and best regards,



Torben



Example:



from dune.grid import onedGrid, cartesianDomain
from dune.fem.space import dgonb as dgSpace
from dune.femdg import createLimiter

import numpy as np
import matplotlib.pyplot as pyplot


numCells = 20
cellSize=1/20
v = [cellSize*i for i in range(int(numCells)+1)]
e = [(i,i+1) for i in range(1,len(v))]
gridView = onedGrid(constructor={"vertices":v, "simplices":e})

spaceRho = dgSpace(gridView, order=1,storage="fem")
rho = spaceRho.interpolate(0.0,name="rho")

a = np.array([ 9.66522837e-01, -3.53578012e-03,  9.55669630e-01, -1.08714750e-02, 9.26611586e-01, -1.83159778e-02,  7.20361930e-01, -1.34572790e-01, 1.73962463e-01, -1.61475760e-01,  5.55255520e-03, -8.51611320e-03, 8.34673602e-04, -4.74279488e-04,  2.63059452e-04, -9.00403039e-05, 1.26471902e-04, -1.98571619e-05,  9.47407424e-05, -1.97633785e-14, 9.47407401e-05,  1.97633341e-14,  1.26471900e-04,  1.98571621e-05, 2.63059451e-04,  9.00403041e-05,  8.34673602e-04,  4.74279489e-04, 5.55255520e-03,  8.51611320e-03,  1.73962463e-01,  1.61475760e-01, 7.20361930e-01,  1.34572790e-01,  9.26611586e-01,  1.83159778e-02, 9.55669630e-01,  1.08714750e-02,  9.66522837e-01,  3.53578012e-03])
for i in range(len(a)):
    rho.dofVector[i] = a[i]

pyplot.plot(gridView.tesselate()[0], rho.pointData(), '-')

tmp = rho.copy()
limiterOp = createLimiter(spaceRho)

limiterOp(rho,tmp)

pyplot.plot(gridView.tesselate()[0], tmp.pointData(), '-')
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune-fem/attachments/20210415/fc59d784/attachment.htm>