1d transport assertion error

1

I am trying to compute the steady state of a substance that is transported along a 1d velocity field, with prescribed rate of inflow:

0 = div(velocity*rho)
# boundary condition:
velocity(0)*rho(0) = fixed

I translated this to the following code:

from fenics import *

class InFlowDomain(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary
        # return on_boundary

mesh = UnitIntervalMesh(5)
V = FunctionSpace(mesh, "DG", 0)
velocity = Expression(("1+x[0]",), degree=2, domain=mesh)
v   = TestFunction(V)
rho = Function(V)
n = FacetNormal(mesh)
F = inner(jump(velocity*rho, n), v)*ds + inner(div(velocity * rho), v) * dx
inflow_rate = 2.0 / velocity(0)
bc = DirichletBC(V, inflow_rate, InFlowDomain())
u = Function(V)
solve(F==0, u)

However this throws and error:


/usr/lib/python3/dist-packages/ffc/uflacs/analysis/balancing.py in balance_modified_terminal(expr)
     52     layers = [expr]
     53     while not expr._ufl_is_terminal_:
---> 54         assert expr._ufl_is_terminal_modifier_
     55         expr = expr.ufl_operands[0]
     56         layers.append(expr)

AssertionError:

Suggestions whats wrong or different approach to solve the transport problem would be appreciated.

2

Some issues I see right away are:

  1. The measure ds integrates over exterior facets (on the boundary of the mesh), so the interpretation of jump is not clear, since there’s no “other side”.

  2. Assuming the desired integral is over interior facets, the appropriate measure is dS (with a capital S), but then, because v is discontinuous at interior facets, the given integrand is ill-defined.

  3. The unknown function u in the nonlinear solve is defined after and does not appear in the residual F. (Should it perhaps be solve(F==0, rho)? Note that the usage of the nonlinear solver is different from that of the linear solver, where the residual involves a TrialFunction and then the result of the solve goes into a new Function.)

1 Like
3

Thanks for pointing these out. I fixed them and now I get a solution, which is rho = 0. Okay I guess I the problem is that rho is allowed to be discontinuous at the boundary, so my boundary condition is irrelevant. How can I fix this?

from fenics import *

class InFlowDomain(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary
        # return on_boundary

mesh = UnitIntervalMesh(5)
V = FunctionSpace(mesh, "DG", 0)
velocity = Expression(("1+x[0]",), degree=2, domain=mesh)
velocity = Expression(("1",), degree=2, domain=mesh)
v   = TestFunction(V)
rho = Function(V)
n = FacetNormal(mesh)
F = inner(jump(velocity*rho, n), avg(v))*dS + inner(div(velocity * rho), v) * dx# + inner(jump(rho), avg(v)) * ds
inflow_rate = 2.0 / velocity(0)
bc = DirichletBC(V, inflow_rate, InFlowDomain())
u = Function(V)
solve(F==0, rho)
plot(rho)
4

So, there are a few things going on. If you want to use a strong BC, you need to consider:

  1. bc is not being passed to solve().
  2. For DG methods, you need to pass the optional argument "pointwise" to the DirichletBC constructor.
  3. For a degree-0 DG space, the degree of freedom to which you want to apply the BC is in the center of the first element, so the InflowDomain won’t include it.

However, it is, in some sense, more consistent with the DG idea to enforce the BC on the inflow domain weakly, in the same way that you enforce inter-element continuity, but where the solution on the “missing” side opposite the boundary is replaced with the desired boundary data, and v is evaluated on the inflow boundary.

The treatment of v in the weakly-enforced inflow BC also brings me to another issue, which is that the formulation given is not stable. The dS terms should be formulated in an “upwind” way, so that v is being evaluated on the side of the interior facet for which velocity is flowing in. (Basically, you want the degree-0 DG formulation for pure advection to give you the same algebraic equations as the classical upwind finite difference scheme.)

2 Likes
5

Thanks for the great help! So I think on paper I think I see how to do the “upwind” formulation.
I need to replace

inner(jump(velocity*rho, n), avg(v))*dS

by

inner(jump(velocity*rho, n), right_hand_side_cell_value(v))*dS

Also for the weak boundary condition, I need to add

inner(rho - inflow_rate, v) * dx(first cell only)

How to implement these using fenics? Here is what I tried, sadly the solution is still all zero.

from fenics import *

nknots = 5
class IndicatorFirstCell(UserExpression):
    def eval(self, value, pt):
        if pt[0] <= (1.0/nknots):
            value[0] = 1.0
        else:
            value[0] = 0.0
        
    def value_shape(self):
        return ()

mesh = UnitIntervalMesh(nknots)
V = FunctionSpace(mesh, "DG", 0)
velocity = Expression(("1",), degree=2, domain=mesh)
v   = TestFunction(V)
rho = Function(V)
n = FacetNormal(mesh)

inflow_rate = 2.0 / velocity(0)

def plus(v):
    return jump(v) + 2*avg(v)

F = inner(jump(velocity*rho, n), plus(v))*dS 
+ inner(rho * IndicatorFirstCell(), v) * dx - inner(inflow_rate * IndicatorFirstCell(), v) * dx

solve(F==0, rho)
plot(rho)
1 Like