Evaluating time dependent Neumann boundary condition using tensors

I/O
1

Hi everyone,
I am writing a time dependent coupled temperature, flow, and pressure solver based on the Cahn-Hilliard FEniCS demo. I have boundary conditions and source terms that depend on time as well, and I have been using classes to define these. For example,

class uExactDef(UserExpression):
    def __init__(self, t, **kwargs):
        super().__init__(kwargs)
        self._t = t

    def eval_cell(self, value, x, ufc_cell):
        value[0] = -np.exp(-c*self._t(x))*np.sin(2*np.pi*(1-x[1]))
        value[1] = np.exp(-c*self._t(x))*np.sin(2*np.pi*x[0])

    def value_shape(self):
        return(2,)

to define the exact solution I am comparing against. Then,

t = Constant(0.0)
uexact = uexactDef(t)

to evaluate it the first time. I use this class for Dirichlet boundary conditions successfully, for example

bcu_wall = DirichletBC(W.sub(0), uexact, wall)

which works well. The timestepping is then done in a while loop, as in the Cahn-Hilliard demo, i.e.

t2 = +dt
t.assign(t2)

where I use t2 as a helper, (because t.assign(t2)) makes sense to me). This timesteps the boundary condition as well (and other time dependent parts). However, I am also imposing the normal stress at one of the boundaries, for which I need uexact and pexact . I want to do

def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
                          [v[1].dx(0), v[1].dx(1)]]))



# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-p*Identity(2)

b1 = - dot(dot(sigma(u_exact, p_exact), v), n) * ds(1)

where ds(1) is the boundary where I want to impose the condition, and b1 is the boundary term added to the variational form. However, this does not work. I get the error

 error("Cannot determine geometric dimension from expression.")

I am currently using the alternative of defining a separate expression

u_Neumann = Expression(('-exp(-c*t)*sin(2*pi*(1-x[1]))', 'exp(-c*t)*sin(2*pi*x[0])'), c=c, t=0.0, degree=2, domain=mesh)
p_Neumann = Expression('-exp(-c*t)*(cos(2*pi*x[0])*cos(2*pi*(1-x[1])))', c=c, t=0.0, degree=2, domain=mesh)

and then

u_Neumann.t = t2
p_Neumann.t = t2

which then works as expected, namely

b1 = - dot(dot(sigma(u_Neumann, p_Neumann), v), n) * ds(1)

gives no errors and the correct answer but with significant slow down. My questions are:

  • Is there a way to use the class uExactDef to impose the normal stress condition I need without errors?
  • Is there another way of imposing it instead of Expression ? This makes my computations incredibly slow.
  • If I want everything to be consistent, would you suggest leaving everything as Expressions or classes?

NB I realise that I didn’t supply the MWE, but it is because the issue I am having is only relevant for the time dependent coupled case (as the time derivative is in the temperature problem). Therefore my MWE is actually quite long, and I think that my question is more related to the correct form of uexact so that it can be used in sigma(uexact, pexact) . If you require the MWE, I will copy it as well (as simplified as I can).

2

I would guess that any unexpected slowness is actually coming from the UserExpression, not the Expression. The UserExpression executes the Python code in eval_cell every time it is evaluated, whereas the C++ string passed to the Expression constructor is compiled once into a kernel that can be executed more efficiently.

1 Like
3

I see! This is incredibly helpful. I have around 15 Expressions that are time dependent, so my approach would be to do, inside the while time-stepping loop,

TimeDepExpr.t = t2
TimeDepExpr2.t = t2
...
TimeDepExpr15.t = t2

in order to advance each one in time. Is there a way to do it all in one go? Or is this the best way?

4

It might be a little cleaner to put all of the Expressions in a Python list instead of defining them as separate variables, then use a for loop over that list to update their time values.

1 Like
5

And maybe even a little cleaner to wrap t in a Constant, pass it as an argument to each Expression, and assign it at each step of the loop, e.g.:

from fenics import *

mesh = UnitSquareMesh(1, 1)
t = Constant(0)
e1 = Expression('t', t=t, degree=0)
e2 = Expression('sin(t)', t=t, degree=0)
dx = Measure('dx', mesh)

for tt in range(7):
    t.assign(tt)
    print(assemble(e1*dx), assemble(e2*dx))

producing

0.0 0.0
1.0 0.8414709848078965
2.0 0.9092974268256817
3.0 0.1411200080598672
4.0 -0.7568024953079282
5.0 -0.9589242746631385
6.0 -0.27941549819892586
3 Likes
6

Thank you so much! This makes everything way cleaner (and tt is so much better than t2…)

1 Like