# -*- coding: utf-8 -*-
"""
Limiters for high order DG methods
"""
import numpy as nm
from sfepy.discrete.dg.poly_spaces import iter_by_order
from sfepy.discrete.dg.fields import get_raveler, get_unraveler
from sfepy.base.base import output
MACHINE_EPS = 1e-30
[docs]
def minmod(a, b, c):
"""Minmod function of three variables, returns:
0 , where sign(a) != sign(b) != sign(c)
min(a,b,c) , elsewhere
Parameters
----------
a : array_like
c : array_like
b : array_like
Returns
-------
out : ndarray
"""
seq = (nm.sign(a) == nm.sign(b)) & (nm.sign(b) == nm.sign(c))
res = nm.zeros(nm.shape(a))
res[seq] = nm.sign(a[seq]) * nm.minimum.reduce([nm.abs(b[seq]),
nm.abs(a[seq]),
nm.abs(c[seq])])
return res
[docs]
def minmod_seq(abc):
"""Minmod function of n variables, returns:
0 , where sign(a_1) != sign(a_2) != ... != sign(a_n)
min(a_1, a_2, a_3, ... , a_n) , elsewhere
Parameters
----------
abc : sequence of array_like
Returns
-------
out : ndarray
"""
seq = nm.hstack([nm.sign(x) for x in abc])
seq = seq[:, 0, None] == seq
seq = seq.prod(axis=1).astype(bool)
res = nm.zeros(nm.shape(abc[0]))
res[seq] = nm.sign(abc[0][seq]) * \
nm.minimum.reduce([nm.abs(x[seq]) for x in abc])
return res
[docs]
class DGLimiter:
name = "abstract DG limiter"
def __init__(self, field, verbose=False):
self.field = field
self.extended = field.extended
self.n_el_nod = field.n_el_nod
self.n_cell = field.n_cell
self.ravel = get_raveler(self.n_el_nod, self.n_cell)
self.unravel = get_unraveler(self.n_el_nod, self.n_cell)
self.verbose = verbose
output("Setting up limiter: {} for {} field {}.".format(self.name,
self.field.family_name,
self.field.name))
def __call__(self, u):
raise NotImplementedError("Called abstract limiter")
[docs]
class IdentityLimiter(DGLimiter):
"""Neutral limiter returning unchanged solution."""
name = "identity"
def __call__(self, u):
if self.verbose: output(self.name + " limiter")
return u
[docs]
class MomentLimiter1D(DGLimiter):
""" Moment limiter for 1D based on [1]_
.. [1] Krivodonova (2007):
Limiters for high-order discontinuous Galerkin methods"""
name = "moment_1D_limiter"
def __call__(self, u):
""""
Parameters
----------
u : array_like
raveled solution at time step n in shape (order * n_space_nod, ...)
Returns
-------
u : ndarray
unraveled limited solution
"""
# for convenience do not try to limit FV
if self.n_el_nod == 1:
if self.verbose: output(self.name + " no limiting for FV.")
return u
u = self.unravel(u).swapaxes(0, 1)
idx = nm.arange(nm.shape(u[0, 1:-1])[0])
idx_bc = nm.arange(nm.shape(u[0, :])[0])
nu = nm.copy(u)
tilu = nm.zeros(u.shape[1:])
for ll in range(self.n_el_nod - 1, 0, -1):
tilu[idx] = minmod(nu[ll, 1:-1][idx],
nu[ll - 1, 2:][idx] - nu[ll - 1, 1:-1][idx],
nu[ll - 1, 1:-1][idx] - nu[ll - 1, :-2][idx])
idx = idx[nm.where(abs(tilu[idx] - nu[ll, 1:-1][idx])
> MACHINE_EPS)[0]]
if self.verbose:
output("{} limiting in {} cells out of {} field {}:".
format(self.name, len(idx_bc), self.n_cell, self.field.name))
output(idx_bc)
if len(idx_bc) == 0:
break
nu[ll, 1:-1][idx] = tilu[idx]
return self.ravel(nu.swapaxes(0, 1))[:, 0]
[docs]
class MomentLimiter2D(DGLimiter):
""" Moment limiter for 2D based on [1]_
.. [1] Krivodonova (2007):
Limiters for high-order discontinuous Galerkin methods"""
name = "moment_limiter_2D"
def __call__(self, u):
"""
Parameters
----------
u : array_like
raveled solution at time step n in shape (order * n_space_nod, ...)
Returns
-------
u : ndarray
unraveled limited solution
"""
if self.n_el_nod == 1:
if self.verbose: output(self.name + " no limiting for FV.")
return u
ex = self.extended
nbrhd_idx = self.field.get_facet_neighbor_idx()
inner_mask = nbrhd_idx[:, :, 0] > 0
idx = nm.where(inner_mask.prod(axis=1))[0]
nbrhd_idx = nbrhd_idx[:, :, 0]
u = self.unravel(u).swapaxes(0, 1)
nu = nm.zeros((self.field.approx_order + 1,) * 2 + u.shape[1:])
tilu = nm.zeros(u.shape[1:])
for ll, (ii, jj) in enumerate(iter_by_order(self.field.approx_order, 2,
extended=ex)):
nu[ii, jj, ...] = u[ll]
for ii, jj in reversed(list(iter_by_order(self.field.approx_order, 2,
extended=ex))):
minmod_args = [nu[ii, jj, idx]]
nbrhs = nbrhd_idx[idx]
if ii - 1 >= 0:
alf = nm.sqrt((2 * ii - 1) / (2 * ii + 1))
# right difference in x axis
dx_r = alf * (nu[ii - 1, jj, nbrhs[:, 1]] - nu[ii - 1, jj, idx])
# left differnce in x axis
dx_l = alf * (nu[ii - 1, jj, idx] - nu[ii - 1, jj, nbrhs[:, 3]])
minmod_args += [dx_r, dx_l]
if jj - 1 >= 0:
alf = nm.sqrt((2 * jj - 1) / (2 * jj + 1))
# right i.e. element "up" difference in y axis
dy_up = alf * (nu[ii, jj - 1, nbrhs[:, 2]] - nu[ii, jj - 1, idx])
# left i.e. element "down" difference in y axis
dy_dn = alf * (nu[ii, jj - 1, idx] - nu[ii, jj - 1, nbrhs[:, 0]])
minmod_args += [dy_up, dy_dn]
tilu[idx] = minmod_seq(minmod_args)
idx = idx[nm.where(abs(tilu[idx] - nu[ii, jj, idx]) > MACHINE_EPS)[0]]
if self.verbose:
output("{} limiting {} in {} cells out of {} field {}:".
format(self.name, (ii, jj), len(idx), self.n_cell, self.field))
output(idx)
if len(idx) == 0:
break
nu[ii, jj, idx] = tilu[idx]
resu = nm.zeros(u.shape)
for ll, (ii, jj) in enumerate(iter_by_order(self.field.approx_order, 2,
extended=ex)):
resu[ll] = nu[ii, jj]
return self.ravel(resu.swapaxes(0, 1))[:, 0]
[docs]
class ComposedLimiter:
def __init__(self, fields, limiters, verbose=False):
self.fields = fields
self.limiters = [limiter(field, verbose=verbose) for field, limiter in zip(fields, limiters)]
self.limiter_slices = []
start = 0
for field in self.fields:
self.limiter_slices.append((start, start + field.n_bubble_dof))
start += field.n_bubble_dof
def __call__(self, u):
lu = nm.zeros(u.shape, dtype=u.dtype)
for lslice, limiter in zip(self.limiter_slices, self.limiters):
lu[lslice[0]:lslice[1]] = limiter(u[lslice[0]:lslice[1]])
return lu