Adapting Tutorial code ft03 to 2019 syntax: attribute and type errors

Errors
1

I am trying to merge a few ideas from the tutorials described in the 2017 book and simultaneously make the code a bit more Python-style so I can debug with Spyder.

Also, because I know it answers a few questions: I am running Ubuntu 20.04.2 (64-bit) on a mini-server (Meerkat5) with an i7 processor, and I installed FEniCS 2019 via Anaconda/conda-forge because I am an engineer and I like to use the Spyder IDE to debug my code. I do have mshr and matplotlib verified as properly installed and connected from running the previous tutorial codes.

When trying to apply the FEniCS command assign from ft03, I started with the updates and then added the Python-style of importing things (e.g. import fenics as fe instead of from fenics import *). The errors I get are associated with the syntax of last line of the file…

If I use u_n.fe.assign(u) I get:

AttributeError: 'Function' object has no attribute 'fe'

If I use u_n=fe.assign(u) I get:

TypeError: assign(): incompatible function arguments. The following argument types are supported:
    1. (arg0: object, arg1: object) -> None

Invoked with: Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 669), FiniteElement('Lagrange', triangle, 1)), 680)

I understand from the book that I can’t use u_n = u because that would set them to be the same variable; and we need to keep them as two variables.

the pydoc fenics.assign command in a terminal window gives me this:

Help on built-in function assign in fenics:
fenics.assign = assign(...) method of builtins.PyCapsule instance
    assign(arg0: object, arg1: object) -> None
(END)

which I don’t find to be very helpful (possibly because I am still new at Python!)

When I simply use the original u_n.assign(u)
it appears that every single time step gives the same graph
image

Since the original ft03 didn’t include an image of the plot nor a way to visualize it, in neither the book nor online (to my knowledge) - I am unsure of whether the code is actually finding the changes in heat at each step, and I’m unsure how to add a sense of scale to the plot.

Since I am dealing with tutorial code, there’s no minimum working example…but here’s the code that doesn’t have the errors but doesn’t appear to update at each timestep:

"""
FEniCS tutorial demo program: Time-dependent 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
"""
# Updated from ft03_heat.py file
# Last changed: 2021-06-09 by Kristin Schaefer
#  Revised for syntax updates in Python 3.8 & FEniCS 2019.1.0
#  and to include ft03 fixes at https://github.com/hplgit/fenics-tutorial/pull/74/commits/9c67e6d0c844a17f02b46a47d1bcdf2029c6fb25
#  ~> 2017 Tutorial book Ch. 3.1 (test problem 1 in 3.1.3)

import fenics as fe
import numpy as np
import matplotlib.pyplot as plt

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

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

# Define boundary condition
u_D = fe.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 = fe.DirichletBC(V, u_D, boundary)

# Define initial value ~ two choices for variational problem (use # to suppress one u_n)
u_n = fe.interpolate(u_D, V)   # recovers exact solution to machine precision (allows verification of code)
# u_n = project(u_D, V)   # computes approximate values at nodes (u_0 as L^2 projection into the FE space)

# Define variational problem: define F, ask for a & L
u = fe.TrialFunction(V)
v = fe.TestFunction(V)
f = fe.Constant(beta - 2 - 2*alpha)

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

# Create VTK format file for saving solutions as time iterates
vtkfile = fe.File("ft03/heat_solution.pvd")

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

    # Update expression objects with current time
    t += dt
    u_D.t = t

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

    # Save to file and Plot solution to screen
    vtkfile << (u, t)   # appends heat in body and time to file
    fe.plot(u, title="ft03: Solution u")
    plt.show()

    # Compute error at vertices by 
    u_e = fe.interpolate(u_D, V)
    error = np.abs(u_e.vector() - u.vector()).max()
    print('t = %.2f: error = %.3g' % (t, error))

    # Update previous solution, keeping u and u_n as separate variables
    u_n.assign(u)
1 Like
2

Using u_n.assign(u) is the correct thing to do. I think what confuses you is that you do not have a colormap in your plot (i.e. the scale changes per plot).
Adding:

    plt.figure()
    p = fe.plot(u, title=f"ft03: Solution u t:{t:.2e}")
    plt.colorbar(p ,cmap='jet')
    # Optional save to png
    plt.savefig(f"sol_{t}.png")

Should show you the different ranges for each plot.

3 Likes
3

Thanks so much Jorgen! Now the first one looks like this:
image

and the last is:
image

which shows the progression that I couldn’t “see” before.

1 Like
4

@dokken I just noticed a the following warning:

MatplotlibDeprecationWarning: The 'cmap' parameter to Colorbar has no effect 
because it is overridden by the mappable; it is deprecated since 3.3 and will 
be removed two minor releases later.

about the line plt.colorbar(p, cmap='jet')
I verified that if I changed it to plt.colorbar(p, cmap='coolwarm') there was no effect, indeed the coloring looks to be the default “viridis” rather than “jet”.

After digging around the matplotlib help and a small, interesting stumble onto why jet may not be good in the case of greyscale and viridis is fine.

I found a few examples (here’s one) where it seemed like the line should actually be two:

plt.imshow(p, cmap='coolwarm')
plt.colorbar()

But this gives a different TypeError: Image data of dtype object cannot be converted to float

The addition of the scale bar is the important thing, and that works wonderfully, but I was wondering how much more time is worth spending on trying to get the figure colors to actually respond to the request… especially since lots of visualization is now done in ParaView and this plot is really just a check for functional code…

Should I just use plt.colorbar(p) and call it a day?

5

I personally like viridis more than other color maps (I use it in paraview as well). I wouldn’t spend too much time tweaking matplotlib visualization, as it is only useful for smaller 2D problems.

2 Likes
6

This is the usage that I employ. I think plt.colorbar takes colormap information from the object p. (This makes sense if you think about it; if p is plotting some field with the colormap viridis, you wouldn’t want your colorbar to be displayed using jet…)

1 Like
8

You should not create a new Expression at every timestep because the boundary condition will not be able to see that you have updated the expression. Likewise, you should not interpolate the initial condition at every timestep, because the variational form will not see the updated initial condition. Instead you should use

for i in range(num_steps):
    
    # Update current time
    t = (i+1)*dt
   
    u_D.t = t
    
    # Compute solution
    solve(a == L, u, bc)
    # Compute error at vertices
    E1 = errornorm(u_D,u,'H1')
    E = E + E1*E1*dt
    E2 = sqrt(E) 
    print(E2)
    
    # Update previous solution
    u_n.assign(u)
1 Like