Updating code from 2016, my current error is in assigning variables for interpolation

1

I am new to FEniCS. I found a code in a paper I am attempting to update. I have worked on it for a while. I have spent time in the tutorials. I am running Python 3.7.3 with FEniCS installed through Conda.

The code is getting stuck on the lines beneath #initialize solution according to exact solution.

I would also welcome direction to specific resources to help update this code.

Thank you and I apologize for what I am sure is a very simple problem.

from dolfin import *
import numpy
from ufl import nabla_grad
from ufl import nabla_div

def gravity(u):
    '''
    right hand side of Eq. 7
    '''
    Ra = 1.0
    val = as_vector([0.0, Ra*u])
    return val


nx = 10
order = 1
t = 0
dt = 0.01
T = 0.1

mesh = UnitSquareMesh(nx, nx)

DGO = FiniteElement("DG", mesh.ufl_cell(), order)
BDM = FiniteElement("BDM", mesh.ufl_cell(), order+1)
DG1  = FiniteElement("DG", mesh.ufl_cell(), order+1)

W = FunctionSpace(mesh, MixedElement([DGO, BDM, DG1]))

#source terms
f_1 = Expression('(sin(x[0]) + cos(x[1]))*cos(t)', t=t, element = DGO)
f_2 = Expression(('0.0','(-sin(x[0]) - sin(x[1]))*cos(t)'), t=t, element = BDM)
f_3 = Expression('-(sin(x[0]) +sin(x[1]))*sin(t) + (-cos(x[0])*cos(x[0]) + sin(x[1])*cos(x[1]))*cos(t)*cos(t) + (sin(x[0])+sin(x[1]))*cos(t)', t=t, element = DG1)

# exact solutions
p_exact = Expression('(sin(x[0]) + cos(x[1]))*cos(t)', t=t, element = DGO) # pressure
uu_exact = Expression(('-cos(x[0])*cos(t)','sin(x[1])*cos(t)'), t=t, element = BDM) # velocity
u_exact = Expression('(sin(x[0]) + sin(x[1]))*cos(t)', t=t, element = DG1) # tempreture

#(p, uu, u) = TrialFunctions(W)  # store the solution for current step
#(p0, uu0, u0) = TrialFunctions(W) # store the solution for previous step

w = Function(W)
w0 =Function(W)

#initialize solution according to exact solution
assign(w.sub(0), interpolate(p_exact, DGO))
assign(w.sub(1), interpolate(uu_exact, BDM))
assign(w.sub(2), interpolate(u_exact, DG1))
assign(w0.sub(0), interpolate(p_exact, DGO))
assign(w0.sub(1), interpolate(uu_exact, BDM))
assign(w0.sub(2), interpolate(u_exact, DG1))

    # split w into p (pressure), uu (velocity) and u (temperature)
(p, uu, u) = w.split()
(p0, uu0, u0) = w0.split()

    # define test functions
(q, vv, v) = TestFunctions(W)

n = FacetNormal(mesh)

    # penalty weights
alpha = Constant(5000.0)
gamma = Constant(10000.0)

    # cell size
h = CellDiameter(mesh)

    # Eq. 20, flow
F_flo = nabla_div(uu)*q*dx - f_1*q*dx
    #F_flo = inner(nabla_div(uu),q)*dx - inner(f_o91,q)*dx

    # Eq. 21, velocity
F_vel = (dot(uu,vv) - div(vv)*p - inner(gravity(u),vv) - inner(f_2,vv))*dx + dot(n,vv)*p_exact*ds

    # un = un on the outflow facet, otherwise 0
un = (dot(uu,n) +abs(dot(uu,n))) / 2.0

    # Eq. 26 , uu is not divergence-free here
a_a = dot(grad(v), -uu*u)*dx - u*v*nabla_div(uu)*dx - f_3*v*dx + dot(jump(v), un('+')*u('+') - un('-')*u('-'))*dS + dot(v,un*u)*ds

    # Eq. 28
a_d = dot(grad(v),grad(u))*dx - dot(avg(grad(v)),jump(u,n))*dS - dot(jump(v,n), avg(grad(u)))*dS\
        - dot(grad(v),n*u)*ds - dot(v*n, grad(u))*ds\
        + (alpha('+')/h('+'))*dot(jump(v,n),jump(u,n))*dS + (gamma/h)*v*u*ds\
        + u_exact*dot(grad(v),n)*ds - (gamma/h)*u_exact*v*ds

a_tim = (u-u0)/dt*v*dx

    #Eq. 29
F_a_d = a_tim +a_a +a_d

    #Eq. 32
F = F_flo + F_vel + F_a_d

    # Time-stepping loop
while t<T:
    t += dt
    print ("nx=%d order=%d t=%f" % (nx, order, t))
        # update time in time-dependent expressions
    f_1.t = t; f_2.t = t; f_3.t = t
    u_exact.t = t; p_exact.t = t; uu_exact.t = t
    # solve the system and store solution in w
    solve(F==0, w)
        # update w0 with w
    w0.vector()[:] = w.vector()
        # save solution in .xml file
    File('Solution.xml') << w

ERROR Output

TypeError                                 Traceback (most recent call last)
<ipython-input-51-abd99c760f99> in <module>()
     44 
     45 #intitialize solution according to exact solution
---> 46 assign(w.sub(0), interpolate(p_exact, DGO))
     47 assign(w.sub(1), interpolate(uu_exact, BDM))
     48 assign(w.sub(2), interpolate(u_exact, DG1))

/usr/lib/python3/dist-packages/dolfin/fem/interpolation.py in interpolate(v, V)
     65         Pv = MultiMeshFunction(V)
     66     else:
---> 67         Pv = Function(V)
     68 
     69     # Compute interpolation

/usr/lib/python3/dist-packages/dolfin/function/function.py in __init__(self, *args, **kwargs)
    232 
    233         else:
--> 234             raise TypeError("Expected a FunctionSpace or a Function as argument 1")
    235 
    236         # Set name as given or automatic

TypeError: Expected a FunctionSpace or a Function as argument 1
2

Hi,
As your error suggest, the last argument has to be a FunctionSpace. You are supplying a FiniteElement. Thus by supplying W.sub(0).collapse() the interpolation will work. However, Since you are working with a MixedSpace, i would suggest using a FunctionAssigner for assigning the solutions to w. See: https://fenicsproject.org/docs/dolfin/1.6.0/python/programmers-reference/cpp/function/FunctionAssigner.html#dolfin.cpp.function.FunctionAssigner

4

I have continued to work on this problem. Any guidance would be very appreciated. Either something in the python 2.7 or in the FENICS from 2015 update are tripping the simulation. Any thoughts are greatly appreciated!

The error output is:

DG0-BDM1-DG1

nx =   4
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-31-e7fbb66509a9> in <module>()
    171         h.append(1.0/nx)
    172         order = orders[case]
--> 173         temp = mms_solver(nx, order,t, dt, T)
    174         E.append(temp[0])
    175         DOF.append(temp[1])

<ipython-input-31-e7fbb66509a9> in mms_solver(nx, order, t, dt, T)
     94 
     95     # cell size
---> 96     h = CellSize(mesh)
     97 
     98     # weak form of flow equation

NameError: name 'CellSize' is not defined

The current code is:

from dolfin import *
from math import log as ln
import time
import numpy as np

def calerrornorm(u_e, u, Ve):
    u_Ve = interpolate(u, Ve)
    u_e_Ve = interpolate(u_e, Ve)
    e_Ve = Function(Ve)
    e_Ve.vector()[:] = u_e_Ve.vector().array() - \
                       u_Ve.vector().array()
    error = e_Ve**2*dx
    return sqrt(assemble(error))

def gravity(u):
    Ra = 1.0
    val = as_vector([0.0, Ra*u])
    print("gravity")
    return val

def mms_solver(nx, order, t, dt, T):
    nx = 10
    order = 1
    t = 0
    dt = 0.01
    T = 0.1

    cost0 = time.time()
    # create a unit saure mesh (nx * nx)
    mesh = UnitSquareMesh(nx, nx)

    # define function spaces
    # for velocity
    BDM = FiniteElement("BDM", mesh.ufl_cell(), order+1)
    # for pressure
    DG0 = FiniteElement("DG", mesh.ufl_cell(), order)
    # for temperature
    DG1 = FiniteElement("DG", mesh.ufl_cell(), order+1)
    W = FunctionSpace(mesh, MixedElement([BDM, DG0, DG1]))
    
    # spaces to interpolate exact solutions
    VCG = VectorFunctionSpace(mesh, "CG", order+3)
    CG = FunctionSpace(mesh, "CG", order+3)

    # manufactured source terms
    f_1 = Expression('(sin(x[0]) + cos(x[1]))*cos(t)', t=t, element=DG0)
    f_2 = Expression(('0.0', '(-sin(x[0]) - sin(x[1]))*cos(t)'), t=t, element=BDM)
    f_3 = Expression('-(sin(x[0]) + sin(x[1]))*sin(t) +(-cos(x[0])*cos(x[0]) \
        + sin(x[1])*cos(x[1]))*cos(t)*cos(t) + (sin(x[0]) + sin(x[1]))*cos(t)', t=t, element=DG1)

    # exact solution expressions
    u_exact = Expression('(sin(x[0])+ sin(x[1]))*cos(t)', t=t, element=DG0)
    p_exact = Expression('(sin(x[0]) + cos(x[1]))*cos(t)', t=t, element=BDM)
    uu_exact = Expression(('-cos(x[0])*cos(t)', 'sin(x[1])*cos(t)'), t=t, element=DG1)

    # w store current solution, and w0 stores solution from the previous step
    w = Function(W)
    w0= Function(W)
        
    # initialize solution according to exact solution expressions
    #Original
    assign(w.sub(0), interpolate(uu_exact, W.sub(0).collapse()))
    assign(w.sub(1), interpolate(p_exact, W.sub(1).collapse()))
    assign(w.sub(2), interpolate(u_exact, W.sub(2).collapse()))
    assign(w0.sub(0), interpolate(uu_exact, W.sub(0).collapse()))
    assign(w0.sub(1), interpolate(p_exact, W.sub(1).collapse()))
    assign(w0.sub(2), interpolate(u_exact, W.sub(2).collapse()))
    (uu, p, u) = split(w)
    (uu0, p0, u0) = split(w0)

    # define test functions
    (vv, q, v) = TestFunctions(W)

    # n is unit normal vector to facet
    n = FacetNormal(mesh)

    # penalty terms on interior and boundary faces
    alpha = Constant(500000.0)
    gamma = Constant(1000000.0)

    # cell size
    h = CellSize(mesh)

    # weak form of flow equation
    F_flo = nabla_div(uu)*q*dx - f_1*q*dx
    # weak form of darcy velocity equation
    F_vel = (dot(uu, vv) - div(vv)*p - inner(gravity(u), vv) - inner(f_2, vv) )*dx + dot(n, vv)*p_exact*ds

    # un = un on outflow facet, otherwise 0
    un = (dot(uu, n) + abs(dot(uu, n)))/2.0

    # weak form of advectin term (Note uu is not divergence free)
    a_a = dot(grad(v), -uu*u)*dx - u*v*nabla_div(uu)*dx + dot(jump(v), un('+')*u('+') - un('-')*u('-') )*dS\
          + dot(v, un*u)*ds - f_3*v*dx
    # weak form of diffusion term
    a_d = dot(grad(v), grad(u))*dx\
          - dot(avg(grad(v)), jump(u, n))*dS - dot(jump(v, n), avg(grad(u)))*dS\
          - dot(grad(v), n*u)*ds - dot(v*n, grad(u))*ds\
          + (alpha('+')/h('+'))*dot(jump(v, n), jump(u, n))*dS + (gamma/h)*v*u*ds\
          + u_exact*dot(grad(v), n)*ds - (gamma/h)*u_exact*v*ds
    # weak form of time gradient term
    a_tim = (u - u0)/dt*v*dx
    # weak form of advection-diffusion equation
    F_a_d = a_tim + a_a + a_d 
    # weak form of the whole system
    F = F_flo + F_vel + F_a_d
    # time-stepping loop
    while t < T:
        t +=dt
        print ("nx=%d order=%d t=%f" % (nx, order, t))
    # update time in time-dependent expressions
    f_1.t = t
    f_2.t = t
    f_3.t = t
    u_exact.t = t
    p_exact.t = t
    uu_exact.t = t
    # solve the system and store solution in w
    solve(F==0, w, solver_parameters={"newton_solver": {"linear_solver": "gmres", "preconditioner": "ilu"}})
    # update w0 with w
    w0.vector()[:] = w.vector()
    open("solution.xml", 'w+') << w
    
    # calculate differences between numerical and analytical solutions
    (uu, p, u) = w.split()
    p_error = calerrornorm(p_exact, p, CG)
    uu_error = calerrornorm(uu_exact, uu, VCG)
    u_error = calerrornorm(u_exact, u, CG)
    # number of degrees of freedom
    dof = len(w.vector().array())
    # return errors and number of dof
    cost = time.time() - cost0
    return [[p_error, uu_error, u_error], dof, cost]

# file to store convergence data
output_handle = open('mms_table.txt', 'w+')
# low order element: 0, high order element: 1
orders = [0, 1]
# mesh
nxs = [4, 8, 16, 32]
for case in range(len(orders)):
    output = 'DG%d-BDM%d-DG%d\n' %(orders[case], orders[case]+1, orders[case]+1)
    print (output)
    # h is to store element size
    h = []
    # E is to store numerical errors
    E = []
    # DOF is to store number of degrees of freedom
    DOF = []
    # cost is to store computational cost (in seconds per step)
    cost = []
    # dt is time step length
    dt = 1.0e-5
    # t is simulation start time
    t = 4.0
    # steps to run
    steps = 1.0e2
    # T is simulation stop time
    T = t + dt*steps
    for nx in nxs:
        print ("nx =  ", nx)
        h.append(1.0/nx)
        order = orders[case]
        temp = mms_solver(nx, order,t, dt, T)
        E.append(temp[0])
        DOF.append(temp[1])
        cost.append(temp[2])

    # r is to store convergence rate
    r = np.zeros((len(nxs), 3))
    # calculate convergence rate
    for i in range(3):
        for j in range(len(E)):
            if j>0:
                r[j][i] = ln(E[j][i]/E[j-1][i])/ln(h[j]/h[j-1])
            else:
                r[j][i] = 0.0
    # print out output and store in .txt file.
    output += '-'*97 + '\n'
    output += '%8s | %8s | %8s | %8s | %8s | %8s | %8s | %8s | %8s\n'\
           % ('h', 'error(p)', 'rate(p)', 'error(u)', 'rate(u)', 'error(T)', 'rate(T)', 'DOF', 'cost')
    for j in range(len(E)):
        output += '%8s | %8.2E | %8.2f | %8.2E | %8.2f | %8.2E | %8.2f | %8d | %8.2E\n'\
           % ('1/'+str(nxs[j]), E[j][0], r[j][0], E[j][1], r[j][1], E[j][2], r[j][2], DOF[j], cost[j]/steps)
    output += '-'*96 + '\n'
    print (output)
    output_handle.write(output)
output_handle.close()
print ("DONE")
5

Cellsize is deprecated. Use 2*Circumradius. I would suggest checking the Changelog, as there has been many changes in the last 4 years.

6

Thank you. This should keep me busy for a while. That resolved that issue.

7

Hi,
Did you manage to update the code so it works without any errors in the current version of FEniCS?
If so, could you share your code please, because I tried to fix the errors of the original code for the whole day now, without much success…