Expression and class(UserExpression) don't yield the same results

variational formulation
1

Dear community,

due to code robustness I would like to use the same class for the initial values and my Manufactured Solution.

For this example case, I used the ft_03_heat.py example of the tutorial. The Manufactured Solution is implemented as following:

u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)

To my surprise, the following implementation as a class doesn’t yield the same results

class U_D(UserExpression):
    
    def __init__(self, t, alpha, beta, **kwargs):
        super().__init__(**kwargs)
        self.t = t
        self.alpha = alpha
        self.beta = beta
        
    def update(self, t):
        self.t = t
    
    def eval(self, value, x):
        
        u_e = 1.0 + x[0]*x[0] + self.alpha*x[1]*x[1] + self.beta*self.t
    
        value = u_e

like the execution of the hidden MWCE shows.

MWCE 
"""
FEniCS tutorial demo program: Heat equation with Dirichlet conditions.
Test problem is chosen to give an exact solution at all nodes of the mesh.
  u'= Laplace(u) + f  in the unit square
  u = u_D             on the boundary
  u = u_0             at t = 0
  u = 1 + x^2 + alpha*y^2 + \beta*t
  f = beta - 2 - 2*alpha
"""

from fenics import *
import numpy as np

t = 0
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 boundary condition

class U_D(UserExpression):
    
    def __init__(self, t, alpha, beta, **kwargs):
        super().__init__(**kwargs)
        self.t = t
        self.alpha = alpha
        self.beta = beta
        
    def update(self, t):
        self.t = t
    
    def eval(self, value, x):
        
        u_e = 1.0 + x[0]*x[0] + self.alpha*x[1]*x[1] + self.beta*self.t
    
        value = u_e


u_D = U_D(t, alpha, beta, degree = 2)
"""   
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)
"""
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)

for n in range(num_steps):

    # Update current time
    t += dt
    u_D.update(t)
    #u_D.t = t

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

    # Plot solution
    plot(u)

    # Compute error at vertices
    u_e = interpolate(u_D, V)
    L2_error = errornorm(u_e, u, 'L2', degree_rise=3)
    vertex_values_u = u.compute_vertex_values(mesh)
    u_avg = np.average(vertex_values_u)
    print('t = %.2f: u_avg = %.5f, error = %.3g' % (t, u_avg, L2_error))

    # Update previous solution
    u_n.assign(u)

How do I have to adjust the class to yield the same results as the Expression?

Thanks!

2

I haven’t run your code, but I’m pretty sure the above needs the following correction:

    def eval(self, value, x):
        
        u_e = 1.0 + x[0]*x[0] + self.alpha*x[1]*x[1] + self.beta*self.t
    
        value[0] = u_e

In your original implementation you overwrite the variable value whereas you should assign quantities to its members.

2 Likes