Errors in use of conditional expressions in variational form

variational formulation
1

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

2

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

3

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.

4

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.

5

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.

6

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

7

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.

8

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)

?

9

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

10

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: .

2 Likes