Errors in use of conditional expressions in variational form

Hello

First of all, I installed FEniCS on a Windows 10 subsystem for Linux with an Ubuntu terminal and I am working with version 2019.1.0. I am attempting to test the error of a fully discretised mixed finite element scheme for a transformed and linearised form of Richards Equation. The exact solution against which the error is measured is:

u_{\text{ex}}(x, y, t) = \begin{cases}-\frac{2(e^{s}-1)}{e^{s}+1},\qquad\text{for } s\geq0,\\ -s,\qquad\ \qquad\text{for } s<0,\end{cases}

with s=x-y-t.

The variational form used in the scheme is as follows, find \left(u^{n,i}_h, \vec{q^{n,i}_h}\right)\in W_h\times V_h such that:

k\left(u^{n,i}_h, w_h \right)_{L^2} + \tau\left(\nabla\cdot\vec{q^{n,i}_h}, w_h\right)_{L^2}=\left(b\left(u^{n-1}_h\right)-b\left(u^{n, i-1}_h\right)+ ku^{n, i-1}_h, w_h\right)_{L^2},\quad \forall w_h\in{W_h},

\left(\vec{q^{n, i}_h}, \vec{v_h}\right)_{L^2} - \left(u^{n,i}_h,\nabla\cdot \vec{v_h} \right)_{L^2} = 0, \qquad \forall \vec{v_h}\in V_h.

Here W_h and V_h are the test spaces which are finite subspaces of L^2\left(\Omega\right) and H\left(\text{div}, \Omega\right) respectively. The superscript n denotes the time step which is being considered and the superscript i indicates the stage within the iterative Newton linearisation procedure performed at each time step n. The value k is a constant involved in the Newton linearisation. Also (\cdot,\cdot)_{L^2} is the L^2 inner product and \Omega=[0,1]^2.

The issue I am having relates to the fact that the function b is a conditional function of u, specifically:

b\left(u\right)=\begin{cases}\frac{\pi^2}{2}-\frac{u^2}{2}, \qquad\text{where }u\leq 0,\\\frac{\pi^2}{2},\qquad\qquad \ \text{where } u>0.\end{cases}

I have attempted to express this function using the conditional command but I either generate nan solutions at each time step or I receive the error message:
NotImplementedError: Cannot take length of non-vector expression.
I provide the working examples for each case below.

Generation of nan solutions:

Running this will show that the numerical solution u has value nan at each vertex on the mesh.

from fenics import *
from dolfin import *
from mshr import *
import numpy as np

#*****************range of mesh refinements and step totals********************

dofs = [4, 8, 16, 32]

steps = [25, 50, 100, 200]

#********solving numerically at each time step and mesh size level************

for j in range(len(dofs)):


#******************* Setting up unit square mesh ***************************

    mesh = UnitSquareMesh(dofs[j], dofs[j])

#***************Defining variables and corresponding function spaces**********

    W = FiniteElement('DG', mesh.ufl_cell(), 0)                                    #Trial/test space of Discontinuous Lagrange elements  
    V = FiniteElement('RT', mesh.ufl_cell(), 1)                                    #Trial/test space of Raviart-Thomas elements
    element = MixedElement([W,V])
    M = FunctionSpace(mesh, element)                                               #Creating mixed finite element space
    W_0 = M.sub(0).collapse()                                                      #Rebuilding trial/test space of constant elements 
    V_0 = M.sub(1).collapse()


#*******************Exact solution and boundary condition**********************


    u_Ex = Expression('x[0]-x[1]-t >= 0 ? -2*(exp(x[0]-x[1]-t)-1)/(exp(x[0]-x[1]-t)-1): -(x[0]-x[1]-t)',
                      degree = 3, t = 0)

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(M.sub(0), u_Ex, boundary)

#*******************Initial Condition******************************************

    u_ = project(u_Ex, W_0)

#*****************Setting values of certain parameters************************

    k = 2.0                                                # Linearisation constant

 
#************** Fixing duration of evolution and timestep size****************

    T = 1                                                   #Final time
    N = steps[j]                                            #Number of steps
    tau = T/N                                               #Time step size
    I = 8                                                   #Number of interim steps in linearization method


#*****************Defining variational form************************************

    u_i_ = u_                                                                                     #Taking initial value within linearisation iteration as value from previous timestep

    pi = np.pi
    
    u_i, q_i = TrialFunctions(M)                                                                  #Defining trial and test functions of the mixed formulation

    w, v = TestFunctions(M)

    bn_ = conditional(le(0, u_), pow(pi,2)/2 - pow(u_,2)/2, pow(pi,2)/2)                          #The function b evaluated at the value of u from the previous time step n-1
    
    bni_ = conditional(le(0, u_i_), pow(pi,2)/2 - pow(u_i_,2)/2, pow(pi,2)/2)                     #The function b evaluated at the value of u from the previous step of the nested iteration (n, i-1)
    
    a1 = k*(u_i*w*dx) + tau*((div(q_i))*w*dx)

    L1 = (bn_- bni_ + k*u_i_)*w*dx 
    
    a2 = dot(q_i, v)*dx - u_i*(div(v))*dx

    L2 = 0

    a = a1 + a2

    L = L1 + L2     

#****************************Advancing scheme in time*************************

    t = 0                                                   #Initial time


    uq = Function(M)                                        #Solution

    for n in range(N):
    
        t += tau                                            #Updating current time
    
        u_Ex.t = t                                          #Updating boundary condition and exact solution
               
        u_Ex_nodes = u_Ex.compute_vertex_values(mesh)       #Extracting vertex values of the exact solution
        
        print(u_Ex_nodes)
        
        for i in range(I):                                  #Linearisation iteration
        
            print(i)
        
            solve(a==L, uq, bc)
        
            u_i, q_i = uq.split(True)                       
        
            u_i_.assign(u_i)                               #assigning the current solution to be the previous solution in the next run of the loop

        u_.assign(u_i_)                                    #once iteration is finished the value of u at this time is recorded and entered as the previous value in the next run of the time loop.
        
        u_nodes = u_.compute_vertex_values(mesh)           #storing vertex values of numerical solution at this time step 
        
        print(u_nodes)

NotImplementedError: Cannot take length of non-vector expression:

In this case the error message appears upon definition of the variational form.

from fenics import *
from dolfin import *
from mshr import *
import numpy as np

#*****************range of mesh refinements and step totals********************

dofs = [4, 8, 16, 32]

steps = [25, 50, 100, 200]

#********solving numerically at each time step and mesh size level************

for j in range(len(dofs)):


#******************* Setting up unit square mesh ***************************

    mesh = UnitSquareMesh(dofs[j], dofs[j])


#***************Defining variables and corresponding function spaces**********

    W = FiniteElement('DG', mesh.ufl_cell(), 0)                                    #Trial/test space of Discontinuous Lagrange elements  
    V = FiniteElement('RT', mesh.ufl_cell(), 1)                                    #Trial/test space of Raviart-Thomas elements
    element = MixedElement([W,V])
    M = FunctionSpace(mesh, element)                                               #Creating mixed finite element space
    W_0 = M.sub(0).collapse()                                                      #Rebuilding trial/test space of constant elements 
    V_0 = M.sub(1).collapse()


#*******************Exact solution and boundary condition**********************


    u_Ex = Expression('x[0]-x[1]-t >= 0 ? -2*(exp(x[0]-x[1]-t)-1)/(exp(x[0]-x[1]-t)-1): -(x[0]-x[1]-t)',
                      degree = 3, t = 0)

    print(u_Ex)
    
    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(M.sub(0), u_Ex, boundary)

#*******************Initial Condition******************************************

    u_ = project(u_Ex, W_0)

#*****************Setting values of certain parameters************************

    k = 2.0                                                # Linearisation constant

 
#************** Fixing duration of evolution and timestep size****************

    T = 1                                                   #Final time
    N = steps[j]                                            #Number of steps
    tau = T/N                                               #Time step size
    I = 8                                                   #Number of interim steps in linearization method



#*****************Defining variational form************************************

    u_i_ = u_                                                                                     #Taking initial value within linearisation iteration as value from previous timestep

    pi = np.pi
    
    u_i, q_i = TrialFunctions(M)                                                                  #Defining trial and test functions of the mixed formulation

    w, v = TestFunctions(M)

    bn_ = conditional(le(0, u_), pow(pi,2)/2 - pow(u_,2)/2, pow(pi,2)/2)                          #The function b evaluated at the value of u from the previous time step n-1
    
    bni_ = conditional(le(0, u_i_), pow(pi,2)/2 - pow(u_i_,2)/2, pow(pi,2)/2)                     #The function b evaluated at the value of u from the previous step of the nested iteration (n, i-1)
    
    a1 = k*(u_i*w*dx) + tau*((div(q_i))*w*dx)

    L1 = (bn_(u_)- bni_(u_i_) + k*u_i_)*w*dx 
    
    a2 = dot(q_i, v)*dx - u_i*(div(v))*dx

    L2 = 0

    a = a1 + a2

    L = L1 + L2     

This produces the error:

Traceback (most recent call last):
    File "r_k_m_c_conditional_mwe_notimplerr.py", line 84, in <module>                                                        
        L1 = (bn_(u_)- bni_(u_i_) + k*u_i_)*w*dx
    File "/usr/lib/python3/dist-packages/ufl/exproperators.py", line 341, in _call                                            
        return _eval(self, arg, mapping, component)                                                                           
    File "/usr/lib/python3/dist-packages/ufl/exproperators.py", line 331, in _eval                                            
        return f.evaluate(coord, mapping, component, index_values)                                                            
    File "/usr/lib/python3/dist-packages/ufl/conditional.py", line 250, in evaluate                                           
        c = self.ufl_operands[0].evaluate(x, mapping, component, index_values)                                                
    File "/usr/lib/python3/dist-packages/ufl/conditional.py", line 134, in evaluate                                           
        b = self.ufl_operands[1].evaluate(x, mapping, component, index_values)                                                
    File "/usr/lib/python3/dist-packages/ufl/core/terminal.py", line 71, in evaluate                                          
        return self.ufl_evaluate(x, component, derivatives)                                                                   
    File "/usr/lib/python3/dist-packages/dolfin/function/function.py", line 270, in ufl_evaluate                              
        return self(*x)                                                                                                       
    File "/usr/lib/python3/dist-packages/ufl/core/expr.py", line 406, in __len__                                              
        raise NotImplementedError("Cannot take length of non-vector expression.")                                           NotImplementedError: Cannot take length of non-vector expression.

I accept this is a large post but I would really appreciate any help at all so that I can finally incorporate this conditional expression into my variational form.

Thank you

Your problem is nonlinear. You shouldn’t be evaluating conditional expressions of your basis (TrialFunctions). You should reformulate your problem as a nonlinear residual.

Thank you for replying so quickly. The function b as I have defined is nonlinear and so is the original
continuous problem, but in my variational form I have already applied a linearisation procedure and as a result, am only applying b to the initial condition u_{ex}, to the solutions at previous timesteps u^{n-1}_h, or to previous iterations of the linearisation procedure u^{n,i-1}_h, all of which should be known. From what I can see, I do not attempt at any point in the script, to include b(u^{n,i}_h) (a conditional expression of a basis function) in the variational form. Sorry if I am wrong but I just wanted to make sure before I try something else.

This short script below reproduces the same message:
NotImplementedError: Cannot take length of non-vector expression
perhaps that will better illustrate what it is I am unsure of.

from fenics import *
from dolfin import *
from mshr import *
import numpy as np



p = Expression('-2*(exp(x[0]-x[1]-t)-1)/(exp(x[0]-x[1]-t)-1)',
                  degree = 3, t = 0)

bn = conditional(le(0, p), pow(pi,2)/2 - pow(p,2)/2, pow(pi,2)/2)


print(bn(p))

Thanks a lot for your help.

Must admit I skim read your code.

Why are you treating bn_ as a function and invoking bn_(u_)? this doesn’t make sense. It’s a ufl form.

That’s fair enough, there is a lot.

In my initial attempt I defined bn_ and bni_ as functions then defined the corresponding component of the linear part of the variational form as:
L1 = (bn_- bni_ + k*u_i_)*w*dx.

However when I ran the code and checked the values the numerical solution took across every node at each time step, I saw nan at every node. Running the example titled Generation of nan solutions should show what I mean.

After running into this problem I wondered if I needed to invoke bn_(u_) instead of bn_, so I just tried it.

Hello again.
I have condensed the code slightly by taking out the iterative procedure and time-stepping. I have also altered certain sections while referring to https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/ I understand it’s an old post but there don’t seem to be many more relevant ones out there.

The function p(x,y) is now given by:
p(x,y)=\begin{cases}-\frac{2\left(e^{(x-y)}-1\right)}{e^{(x-y)}+1},\qquad\text{for } x-y\geq 0,\\ y-x,\qquad\qquad\ \ \text{for } x-y< 0.\end{cases}

The function b(p) is:

b(p)=\begin{cases}\frac{\pi^2}{2}-\frac{u^2}{2},\qquad\text{for } p\leq 0,\\ \frac{\pi^2}{2},\qquad\qquad\ \text{for } p> 0.\end{cases}

The variational form of the problem has been simplified to finding \left(u_h, \vec{q}_h\right)\in W_h \times V_h such that

k\left(u_h,w_h\right)_{L^2}+\tau\left(\nabla\cdot\vec{q}_h, w_h\right)_{L^2} = \left(b(p)+kp, w_h\right)_{L^2}\qquad\forall w_h\in W_h

\left(\vec{q}, \vec{v}_h\right)_{L^2}-\left(u, \nabla\cdot\vec{v}_h\right)_{L^2} = 0\qquad\forall \vec{v}_h\in V_h.

The test spaces W_h and V_h are the same spaces that were given in the original post and k is again a constant.

The code I now use gives
Error: Unable to define linear variational problem a(u, v) == L(v) for all v. *** Reason: Expecting the left-hand side to be a bilinear form (not rank 1).

Here is the code

from fenics import *
from dolfin import *
from mshr import *
import numpy as np
import ufl                                                                     #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/


ufl.algorithms.apply_derivatives.CONDITIONAL_WORKAROUND = True                 #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/    

mesh = UnitSquareMesh(16, 16)

#***************Defining variables and corresponding function spaces**********

W = FiniteElement('DG', mesh.ufl_cell(), 0)                                    #Trial/test space of Discontinuous Lagrange elements  
V = FiniteElement('RT', mesh.ufl_cell(), 1)                                    #Trial/test space of Raviart-Thomas elements
element = MixedElement([W,V])
M = FunctionSpace(mesh, element)                                               #Creating mixed finite element space
W_0 = M.sub(0).collapse()                                                      #Retrieving trial/test space of constant elements 
V_0 = M.sub(1).collapse()                                                      #Retrieving trial/test space of vector elements
    
    
#**************************boundary condition*********************************     

u_Ex = Expression('x[0]-x[1] >= 0 ? -2*(exp(x[0]-x[1])-1)/(exp(x[0]-x[1])-1): -(x[0]-x[1])',
                   degree = 3)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(M.sub(0), u_Ex, boundary)

#*******************Initial Condition******************************************

p = project(u_Ex, W_0)                                                         #Projection of u_Ex onto scalar component of mixed space.

#*****************Setting values of certain parameters************************

k = 2.0                                                                        # Linearisation constant

T = 0.14                                                                       #Final time in corresponding parabolic problem
N = 10                                                                         #Number of steps in parabolic problem
tau = T/N                                                                      #Time step size in parabolic problem


#*****************Defining variational form************************************

uq = Function(M)                                                               #Defining trial functions of the mixed formulation

u, q = uq.split(True)

w, v = TestFunctions(M)                                                        #Defining test functions of the mixed formulation

bn_ = conditional(le(0, p), pow(pi,2)/2 - pow(p,2)/2, pow(pi,2)/2)             #The function b acting on the projection u_p of u_Ex 

a1 = k*(u*w*dx) + tau*((div(q))*w*dx)                                          #bilinear form in equation for u 

L1 = (bn_ + k*p)*w*dx                                                          #linear form in equation for u 
    
a2 = dot(q, v)*dx - u*(div(v))*dx                                              #bilinear form in equation for q

L2 = 0                                                                         #linear form in equation for q  

a = a1 + a2

L = L1 + L2     

#*****************Solving mixed problem****************************************

uq = Function(M)                                                               #Solution

solve(a==L, uq, bc)                                                            #Solving mixed variational form 

As before I am using Fenics 2019.10 on a Windows 10 subsystem for linux.

Any help on why I am receiving the error I am is greatly appreciated,

Thanks

Try the following substitutions:

#uq = Function(M)
uq = TrialFunction(M)
#u, q = uq.split(True)
u,q = split(uq)

The difference between a Function and a TrialFunction is that a Function is not considered a free argument in variational forms, whereas TrialFunctions and TestFunctions are. The “rank” is referring to the number of free arguments to the form, and a bilinear form must have two, by definition. The difference between split(uq) and uq.split() is, I think, generally agreed to be a confusing aspect of the API, as discussed on the issue tracker here, although, for constructing variational forms for mixed formulations, you want to use the former.

Thanks for your help. Should I split in the same way after fenics has solved the problem?

I am currently doing:

#*****************Solving mixed problem****************************************

uq = Function(M)                                                               #Solution

solve(a==L, uq, bc)                                                            #Solving mixed variational form 
        
u, q = uq.split(True)                                                          #Splitting solution into functions from the respective components of the mixed space 

Should I instead be doing the following:

#*****************Solving mixed problem****************************************

uq = Function(M)                                                               #Solution

solve(a==L, uq, bc)                                                            #Solving mixed variational form 
        
#u, q = uq.split(True)                                                          #Splitting solution into functions from the respective components of the mixed space 

u, q = split(uq) 

?

When extracting components of the mixed function for plotting/postprocessing, you do want to use the Function.split() method.

Thank you very much. This seems to have resolved my issue. I also had to create the function p using conditional as opposed to expression.

#p_Ex = Expression('x[0]-x[1] >= 0 ? -2*(exp(x[0]-x[1])-1)/(exp(x[0]-x[1])-1): -(x[0]-x[1])',
#                   degree = 3)

x0, x1 = MeshCoordinates(mesh)

p_Ex = conditional(ge(0, x0-x1), -2*(exp(x0-x1)-1)/(exp(x0-x1)-1), -(x0-x1))

My functioning code in full now reads as follows:

from fenics import *
from dolfin import *
from mshr import *
import numpy as np
#import ufl                                                                     #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/


#ufl.algorithms.apply_derivatives.CONDITIONAL_WORKAROUND = True                 #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/    

mesh = UnitSquareMesh(16, 16)

#***************Defining variables and corresponding function spaces**********

W = FiniteElement('DG', mesh.ufl_cell(), 0)                                    #Trial/test space of Discontinuous Lagrange elements  
V = FiniteElement('RT', mesh.ufl_cell(), 1)                                    #Trial/test space of Raviart-Thomas elements
element = MixedElement([W,V])
M = FunctionSpace(mesh, element)                                               #Creating mixed finite element space
W_0 = M.sub(0).collapse()                                                      #Retrieving trial/test space of constant elements 
V_0 = M.sub(1).collapse()                                                      #Retrieving trial/test space of vector elements
    
    
#**************************boundary condition*********************************     

#p_Ex = Expression('x[0]-x[1] >= 0 ? -2*(exp(x[0]-x[1])-1)/(exp(x[0]-x[1])-1): -(x[0]-x[1])',
#                   degree = 3)

x0, x1 = MeshCoordinates(mesh)

p_Ex = conditional(ge(0, x0-x1), -2*(exp(x0-x1)-1)/(exp(x0-x1)-1), -(x0-x1))

def boundary(x, on_boundary):
    return on_boundary

#bc = DirichletBC(M.sub(0), u_Ex, boundary)

#bc = DirichletBC(W_0, u_Ex, boundary)


#*******************Initial Condition******************************************

p = project(p_Ex, W_0)                                                         #Projection of u_Ex onto scalar component of mixed space.

bc = DirichletBC(M.sub(0), p, boundary)

#*****************Setting values of certain parameters************************

k = 2.0                                                                        # Linearisation constant

T = 0.14                                                                       #Final time in corresponding parabolic problem
N = 10                                                                         #Number of steps in parabolic problem
tau = T/N                                                                      #Time step size in parabolic problem


#*****************Defining variational form************************************

#uq = Function(M)                                                               

uq = TrialFunction(M)                                                          #Defining trial functions of the mixed formulation

#u, q = uq.split(True)

u,q = split(uq)

w,v = TestFunctions(M)                                                         #Defining test functions of the mixed formulation

bn_ = conditional(le(0, p), pow(pi,2)/2 - pow(p,2)/2, pow(pi,2)/2)             #The function b acting on the projection u_p of u_Ex 

a1 = k*(u*w*dx) + tau*((div(q))*w*dx)                                          #bilinear form in equation for u 

L1 = (bn_ + k*p)*w*dx                                                          #linear form in equation for u 
  
a2 = dot(q, v)*dx - u*(div(v))*dx                                              #bilinear form in equation for q

L2 = 0                                                                         #linear form in equation for q  

a = a1 + a2

L = L1 + L2     

#*****************Solving mixed problem****************************************

uq = Function(M)                                                               #Solution

solve(a==L, uq, bc)                                                            #Solving mixed variational form 
        
u, q = uq.split(True)                                                          #Splitting solution into functions from the respective components of the mixed space 

print(u)

unodes = u.compute_vertex_values(mesh)                                         #Storing vertex values of numerical solution

print(unodes)                                                                  #Displaying vertex values of numerical solution 

Thanks again Nate and David :slight_smile: .

1 Like