How to define u for t=0 by an UserExpression?

Cou · February 29, 2020, 6:17pm

The minimum working code example is based on the ft03_heat.py tutorial, which describes a time dependent poisson equation and is governed by the following equations:

u′=Δu+f in the unit square
u=u_D on the boundary (Dirichlet boundary conditions)
u=u_0 at t = 0

I implemented u for t=0 by the following UserExpression instead of using an Expression like
u_D = Expression(‘1 + x[0]x[0] + alphax[1]x[1] + betat’, degree=2, alpha=alpha, beta=beta, t=0),
like it’s used in the tutorial:

# define u for t=0
class MyExpression_t0(UserExpression):
    def eval(self, value, x):
        x = SpatialCoordinate(mesh)
        tol = 1E-14
        # u on left boundary
        if x[0] < tol:
            u = Constant(2.5)
        # u on left domain
        elif x[0] > 0 and x[0] <= 0.1:
            u = Constant(0.5)
        # u for the rest of the domain and boundaries
        else:
            u = 0
        return u
    
u_n = MyExpression_t0()

Unfortunately, I get the following error message:

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
UFL conditions cannot be evaluated as bool in a Python context.
Traceback (most recent call last):
  File "ft03_heat_u0.py", line 38, in <module>
    u_n = interpolate(u_D, V)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/interpolation.py", line 71, in interpolate
    Pv.interpolate(v._cpp_object)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 365, in interpolate
    self._cpp_object.interpolate(u)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py", line 53, in wrapped_eval
    self.user_expression.eval(values, x)
  File "ft03_heat_u0.py", line 20, in eval
    if x[0] < tol:
  File "/usr/local/lib/python3.6/dist-packages/ufl/conditional.py", line 46, in __bool__
    error("UFL conditions cannot be evaluated as bool in a Python context.")
  File "/usr/local/lib/python3.6/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.

How am I able to fix this? Of course, the use of an UserExpression isn´t mandatory. It was just what I came up with.

Here is the full MWCE:

from fenics import *

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# define u for t=0
class MyExpression_t0(UserExpression):
    def eval(self, value, x):
        x = SpatialCoordinate(mesh)
        tol = 1E-14
        # u on left boundary
        if x[0] < tol:
            u = Constant(2.5)
        # u on left domain
        elif x[0] > 0 and x[0] <= 0.1:
            u = Constant(0.5)
        # u for the rest of the domain and boundaries
        else:
            u = 0
        return u
    
u_D = MyExpression_t0() 

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Compute solution
    solve(a == L, u, bc)

    # Update previous solution
    u_n.assign(u)

# Hold plot
interactive()

kamensky · March 1, 2020, 5:14pm

The problem here is the line x = SpatialCoordinate(mesh) inside of the eval method. You want to use the x argument that is passed to eval.

Cou · March 4, 2020, 11:06am

Of course, thanks a lot @kamensky!

Unfortunately, u for t=0, that I described by UserExpression, isn’t used by the program.
I tried two cases( a) u=2.5 and u=0.5; b) u=500 and u=100), like seen underneath

# define u for t=0
class MyExpression_t0(UserExpression):
    def eval(self, value, x):
        tol = 1E-14
        # u on left boundary
        if x[0] < tol:
            u = Constant(2.5)
            #u = Constant(500)
        # u on left domain
        elif x[0] > 0 and x[0] <= 0.1:
            u = Constant(0.5)
            #u = Constant(100)
        # u for the rest of the domain and boundaries
        else:
            u = 0
        return u
    
u_D = MyExpression_t0()

and I always get the same solution, which shouldn´t be the case in my opinion.
I suspect, it has something to do with the warning, that is displayed, when running the program:

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.

In my opinion u for t=0 is a tensor of second degree, not a scalar. How do I implement this or is there another problem, I´m not aware of?

Here is my updated MWCE:

from fenics import *

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

vtkfile = File("ft03_heat_u0/solution.pvd")

# define u for t=0
class MyExpression_t0(UserExpression):
    def eval(self, value, x):
        tol = 1E-14
        # u on left boundary
        if x[0] < tol:
            u = Constant(2.5)
            #u = Constant(500)
        # u on left domain
        elif x[0] > 0 and x[0] <= 0.1:
            u = Constant(0.5)
            #u = Constant(100)
        # u for the rest of the domain and boundaries
        else:
            u = 0
        return u
    
u_D = MyExpression_t0() 

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Compute solution
    solve(a == L, u, bc)
    vtkfile << (u, t)

    # Update previous solution
    u_n.assign(u)

Thanks a lot!

kamensky · March 4, 2020, 11:32am

The problem is that the evaluation of the UserExpression should be written to the value array, not returned by the eval method. For a scalar-valued expression, you would just set value[0] = u. (When assigning to a NumPy array, u will be cast to a float, so there is also no need to wrap numerical values as Constants, although it is not harmful in this case.) To get rid of the warning, you can implement the value_shape method of UserExpression, i.e.,

    def value_shape(self):
        return ()

for a scalar-valued expression (although that is assumed by default, as stated in the warning message).

Cou · March 6, 2020, 10:25am

Thanks again for your great support @kamensky!

I’ve got one final question:
I want to save u for t=0 to the .pvd file for further postprocessing in ParaView.
Unfortunately…

u = MyExpression_t0()
t = 0
vtkfile << (u, t)

…results in the following error message:

Traceback (most recent call last):
  File "ft03_heat_u0_Forum_Check.py", line 43, in <module>
    vtkfile << (u, t) 
  File "/usr/local/lib/python3.6/dist-packages/dolfin/io/__init__.py", line 19, in __lshift__
    self.write(u[0], u[1])
RuntimeError: Unable to cast Python instance to C++ type (compile in debug mode for details)

Can you explain me how to fix this issue?
Thanks again!

Updated MWCE:

from fenics import *
import numpy as np

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

vtkfile = File("ft03_heat_u0/solution.pvd")

# define u for t=0
class MyExpression_t0(UserExpression):
    def eval(self, value, x):
        tol = 1E-14
        
        # u on left boundary
        if x[0] < tol:
            u = Constant(500)
            value[0] = u
            
        # u on left domain
        elif x[0] > 0 and x[0] <= 0.1:
            u = Constant(100)
            value[0] = u
            
        # u for the rest of the domain and boundaries
        else:
            u = 0
            value[0] = u
        return value[0]
    def value_shape(self):
        return ()
    
# save u for t=0 to vtkfile    
u = MyExpression_t0()
t = 0
vtkfile << (u, t) 

# Define boundary conditions
def boundary(x, on_boundary):
    return on_boundary

u_D = MyExpression_t0()
bc = DirichletBC(V, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)

for n in range(num_steps):

    # Update current time
    t += dt

    # Compute solution
    solve(a == L, u, bc)
    vtkfile << (u, t)

    # Update previous solution
    u_n.assign(u)

kamensky · March 6, 2020, 3:42pm

If you just want to visualize the initial condition, you could instead write the finite element function u_n to vtkfile after interpolating u_D.

Topic		Replies	Views
Problem in defining User-defined expressions by subclassing variational formulation	7	1250	November 17, 2023
User defined initial condition in heat equation?	1	2202	September 1, 2019
Beginner struggels with object-oriented programming to compute a time-, space- and u-dependent source term variational formulation	5	1514	May 21, 2020
Ufl error when defining variational formulations based on UserExpression variational formulation	3	2111	June 18, 2020
Expression and class(UserExpression) don't yield the same results variational formulation	1	400	July 21, 2022

How to define u for t=0 by an UserExpression?

Related Topics