Euler's method for matrices

Dear community,

In order to calculate a concentration at the borders of the domain, named c_s, by the Euler’s method, I wrote the following code

# dc_s / dt = integral^{Gamma} (u - c_s) dA
def func_c_s(t, c_s, u_prev ):    
    return ((u_prev-c_s)*ds)
   
def euler_c_s( t0, c_s, dt, t, u_prev ):
    temp = -0

    # Iterating till the point at which we 
    # need approximation 
    while t0 < t: 
        temp = c_s
        c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
        t0 = t0 + dt
    return c_s
       
def value_shape(self):
    return ()

, which is run by

while t < T:
    t = t+dt
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)

,whereas c_s_0 is the intial value for t=0 of c_s and described by

class My_c_s_t0(UserExpression):
    def eval(self, value, x):
                
        # boundaries
        if x[0] == 0:
            value[0] = 0.2
        
        elif x[0] == 1:
            value[0] = 0.2

        elif x[1] == 0:
            value[0] = 0.2
        
        elif x[1] == 1:
            value[0] = 0.2
                
        # domain without boundaries
        else:
            value[0] = 0

def value_shape(self):
    return ()
c_s_t0 = My_c_s_t0()
c_s_0 = interpolate(c_s_t0, V)

and u_prev, which is usually determined by solve(a == L, u), but here, in order to reduce the complexity of the MWCE given by

class My_u_t0(UserExpression):
    def eval(self, value, x):
                
        if x[0] >= 0 and x[0] <= 0.4:
            value[0] = u_avg_domain_t0

        elif x[0] >= 0.6 and x[0] <= 1.0:
            value[0] = u_avg_domain_t0
                
        else:
            value[0] = 0
        
def value_shape(self):
    return ()
u_t0 = My_u_t0()
u_prev = interpolate(u_t0, V)

When running the MWCE, I get the following error code

  File "20201217_c_s.py", line 150, in <module>
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
  File "20201217_c_s.py", line 43, in euler_c_s
    c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
TypeError: unsupported operand type(s) for +: 'Function' and 'Form'

I´m unsure, if the error message is a result of a wrong implementation of the Euler´s method (Before, I only used Euler´s method for problems with scalar initial values) or func_c_s can’t handle the matrix of c_s in its current form?

Thanks a lot for your time,
Cou

Full MWCE
from fenics import *
import numpy as np


############### define concentration c_s (on boundaries) for t=0 ##############

class My_c_s_t0(UserExpression):
    def eval(self, value, x):
                
        # boundaries
        if x[0] == 0:
            value[0] = 0.2
        
        elif x[0] == 1:
            value[0] = 0.2

        elif x[1] == 0:
            value[0] = 0.2
        
        elif x[1] == 1:
            value[0] = 0.2
                
        # domain without boundaries
        else:
            value[0] = 0

def value_shape(self):
    return ()


############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
def func_c_s(t, c_s, u_prev ):    
    return ((u_prev-c_s)*ds)
   
def euler_c_s( t0, c_s, dt, t, u_prev ):
    temp = -0

    # Iterating till the point at which we 
    # need approximation 
    while t0 < t: 
        temp = c_s
        c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
        t0 = t0 + dt
    return c_s
       
def value_shape(self):
    return ()  


######################## define u for t=0 #####################################

class My_u_t0(UserExpression):
    def eval(self, value, x):
                
        if x[0] >= 0 and x[0] <= 0.4:
            value[0] = u_avg_domain_t0

        elif x[0] >= 0.6 and x[0] <= 1.0:
            value[0] = u_avg_domain_t0
                
        else:
            value[0] = 0
        
def value_shape(self):
    return ()


######################## define mesh and function space #######################

nx = ny = 16
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'CG', 1)


################################ geometry #####################################

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

class Obstacle_left(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0.0, 0.4)))
    
class Obstacle_right(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0.6, 1.0)))
    
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle_left = Obstacle_left()
obstacle_right = Obstacle_right()

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, space_dim-1)
domains.set_all(0)
obstacle_left.mark(domains, 1)
obstacle_right.mark(domains, 2)

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, space_dim-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)


##################### define constants #########################################
t = 0
t0 = 0
T = 2.0                     
num_steps = 10              
dt = T / num_steps          

u_avg_domain_t0 = 1   
u_t0 = My_u_t0()
c_s_t0 = My_c_s_t0()

u_prev = interpolate(u_t0, V)
c_s_0 = interpolate(c_s_t0, V)

c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))

while t < T:
    t = t+dt
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
    print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))

As you have not supplied a minimal working code example (one that can be executed and reproduce your error by copy-pasting it into a text editor and running it), it slightly tricky to debug.
My best guess is that if you remove the outer brackets on the following term:

your code will work.

1 Like

Thanks @dokken and happy New Year to everyone!

I removed the outer brackets, but unfortunately I´m faced with the same error message.

  File "20210106_c_s.py", line 150, in <module>
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
  File "20210106_c_s.py", line 43, in euler_c_s
    c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
TypeError: unsupported operand type(s) for +: 'Function' and 'Form'

That´s my MWCE:

from fenics import *
import numpy as np


############### define concentration c_s (on boundaries) for t=0 ##############

class My_c_s_t0(UserExpression):
    def eval(self, value, x):
                
        # boundaries
        if x[0] == 0:
            value[0] = 0.2
        
        elif x[0] == 1:
            value[0] = 0.2

        elif x[1] == 0:
            value[0] = 0.2
        
        elif x[1] == 1:
            value[0] = 0.2
                
        # domain without boundaries
        else:
            value[0] = 0

def value_shape(self):
    return ()


############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
def func_c_s(t, c_s, u_prev ):    
    return (u_prev-c_s)*ds
   
def euler_c_s( t0, c_s, dt, t, u_prev ):
    temp = -0

    # Iterating till the point at which we 
    # need approximation 
    while t0 < t: 
        temp = c_s
        c_s = c_s + dt * func_c_s(t0, c_s, u_prev) 
        t0 = t0 + dt
    return c_s
       
def value_shape(self):
    return ()  


######################## define u for t=0 #####################################

class My_u_t0(UserExpression):
    def eval(self, value, x):
                
        if x[0] >= 0 and x[0] <= 0.4:
            value[0] = u_avg_domain_t0

        elif x[0] >= 0.6 and x[0] <= 1.0:
            value[0] = u_avg_domain_t0
                
        else:
            value[0] = 0
        
def value_shape(self):
    return ()


######################## define mesh and function space #######################

nx = ny = 16
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'CG', 1)


################################ geometry #####################################

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

class Obstacle_left(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0.0, 0.4)))
    
class Obstacle_right(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0.6, 1.0)))
    
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle_left = Obstacle_left()
obstacle_right = Obstacle_right()

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, space_dim-1)
domains.set_all(0)
obstacle_left.mark(domains, 1)
obstacle_right.mark(domains, 2)

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, space_dim-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)


##################### define constants #########################################
t = 0
t0 = 0
T = 2.0                     
num_steps = 10              
dt = T / num_steps          

u_avg_domain_t0 = 1   
u_t0 = My_u_t0()
c_s_t0 = My_c_s_t0()

u_prev = interpolate(u_t0, V)
c_s_0 = interpolate(c_s_t0, V)

c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))

while t < T:
    t = t+dt
    c_s = euler_c_s(t0, c_s_0, dt, t, u_prev)
    print('t = %.2f, c_s(0,0) = %.3f' % (t, c_s(0,0)))

The initial c_s is a function, with no integration domain, thus the error.
I am not sure what you are trying to achieve in this loop.

Thanks for pointing that out, because after some thinking I´m unsure, if my approach is logical sound.

Anyway, my goal is to calculate c_s, which is the concentration outside of the domain, described by
\frac{c_s(t)}{dt}=\int_{\Gamma}^{}(u(\mathbf{x},t)-c_s(t))dA
and by implementing the Euler´s method
c_{s}^{n+1}=c_s^{n}+\Delta t\int_{\Gamma}^{}(u(\mathbf{x},t)-c_s(t))dA
like shown here.

Because c_s is a scalar value by definition, but I have to substract c_s from u, I decided to “scale” c_s “up” by using (c_s_t0 = df.zero((nx,ny)) was not working)

class My_c_s_t0(UserExpression):
    def eval(self, value, x):
                
        # boundaries
        if x[0] == 0:
            value[0] = 0
        
        elif x[0] == 1:
            value[0] = 0

        elif x[1] == 0:
            value[0] = 0
        
        elif x[1] == 1:
            value[0] = 0
                
        # domain without boundaries
        else:
            value[0] = 0

def value_shape(self):
    return ()

c_s_t0 = My_c_s_t0()
c_s_0 = interpolate(c_s_t0, V)

in order to obtain the same number of rows and columns as u.
By the way, this is an update to the previous MWCE, because c_s(t=0)=0.

So, besides the eventual lack of logical soundness, how am I able to add an integration domain to c_s?

You would Need to assemble this (so that it becomes a scalar and can be added to a float). However, in your code, c_s is a function (that is spatially varying, as it is defined by a function. I guess you could use the “R” space for “c”: see for instance: Bitbucket