conditional having no effect

Question

I am trying to write a UFL form to represent a flux in a DG problem. I want there to be no flux across an edge if the value of a rho on the side the flow would be into exceeds a maximum value. I thought I could implement this with conditional. The relevant code is:

import fenics as fe
...

compute a flow velocity vector self.v. In what follows, self.rho is a Function on the DG space. Then

        self.flux = self.v * self.rho
        self.vn = fe.max_value(fe.dot(self.v, self.n), 0)
        self.facet_flux1 = (
            self.vn('+')*self.rho('+') - self.vn('-')*self.rho('-')
        )
        self.facet_flux2 = fe.conditional(
            fe.ge(self.rho('+'), self.rho_max),
                fe.min_value(0, self.facet_flux1),
                self.facet_flux1
        )
        self.facet_flux = fe.conditional(
            fe.ge(self.rho('-'), self.rho_max),
            fe.max_value(0, self.facet_flux2),
            self.facet_flux2
        )
        self.rho_flux_jump = -self.facet_flux*fe.jump(self.wrho)*fe.dS

self.rho_flux_jump then becomes one term in a more complicated form. This works in the sense that I get a weak form that I can assemble, and I can even solve the DE. However, self.rho_max has no effect, even if I set it absurdly low, so that fe.ge(self.rho('+-'), self.rho_max) is expected to trigger almost everywhere. The results are indistinguishable from what i get with just:

        self.facet_flux = (
            self.vn('+')*self.rho('+') - self.vn('-')*self.rho('-')
        )

Is this a valid usage of conditional? If not, is there a way to accomplish what I'm after here?

Thanks.

lavery · Answer 1 · Jul 7, 2017

I found a way of doing it. I don't know why my previous attempt didn't work, but this does:

        self.rho_blocked = fe.conditional(
            fe.ge(self.rho, self.rho_max),
            fe.conditional(
                fe.gt(fe.dot(fe.grad(self.rho), self.v), 0),
                0,
                self.rho
            ),
            self.rho
        )
        self.flux = self.v * self.rho_blocked
        self.vn = fe.max_value(fe.dot(self.v, self.n), 0)
        self.facet_flux = (
            self.vn('+')*self.rho_blocked('+') - self.vn('-')*self.rho_blocked('-')
        )
        self.rho_flux_jump = -self.facet_flux*fe.jump(self.wrho)*fe.dS
        ...

The effect of the conditional is to set rho_blocked to 0 wherever rho >= rho_max and the flow velocity points up the gradient of rho. Thus there is no flux into regions where rho exceeds the max.

conditional having no effect

1 Answer