Converting a set of functions to a class

Hello,

I tried to convert the following set of functions

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

    # Iterating till the point at which we 
    # need approximation 
    while t0 < t: 
        temp = c_s

        c_s = c_s + assemble(dt * func_c_s(t0, c_s, u)) 
        t0 = t0 + dt
    return c_s

def func_c_s(t, c_s, u):    
    return (u-c_s)*ds
       
def value_shape(self):
    return ()

to a class, derived from an UserExpression, by the following code

############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
class My_c_s(UserExpression):
    
    def __init__(self, t, u, **kwargs):
        super().__init__(**kwargs)
        self.t = t
        self.u = u   
    
    def euler_c_s(self, t0, c_s, dt, t, u):
        temp = -0

        # Iterating till the point at which we 
        # need approximation 
        while t0 < self.t: 
            temp = c_s

            c_s = c_s + assemble(dt * func_c_s(t0, c_s, self.u)) 
            t0 = t0 + dt
        return c_s

    def func_c_s(self, t, c_s, u):    
        return (self.u-c_s)*ds
       
    def value_shape(self):
        return () 

Replacing the set of functions from above with the class My_c_s(UserExpression) in the running code, which also serves as a MWCE,

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


############################## Euler method ###################################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
class My_c_s(UserExpression):
    
    def __init__(self, t, u, **kwargs):
        super().__init__(**kwargs)
        self.t = t
        self.u = u   
    
    def euler_c_s(self, t0, c_s, dt, t, u):
        temp = -0

        # Iterating till the point at which we 
        # need approximation 
        while t0 < self.t: 
            temp = c_s

            c_s = c_s + assemble(dt * func_c_s(t0, c_s, self.u)) 
            t0 = t0 + dt
        return c_s

    def func_c_s(self, t, c_s, u):    
        return (self.u-c_s)*ds
       
    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)
R = FiniteElement("Real", mesh.ufl_cell(), 0)

################################ 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()
u = interpolate(u_t0, V)

c_s_0 = 0
c_sat = 0.9
#c_s = euler_c_s(t0, c_s_0, dt, t, u)
c_s = My_c_s(t, u, degree=1)  

while t < T:
    t = t+dt
    #c_s = euler_c_s(t0, c_s_0, dt, t, u)
    c_s.t = t
    c_s.u = u
    
    if c_s < c_sat:
        print('t = %.1f: c_s = %.9f' % (t, c_s))
    else:
        print('Solution saturated.')

results in the following error message

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
UFL conditions cannot be evaluated as bool in a Python context.
Traceback (most recent call last):
  File "20210412_c_s_class.py", line 136, in <module>
    if c_s < c_sat:
  File "/usr/local/lib/python3.6/dist-packages/ufl/conditional.py", line 46, in __bool__
    error("UFL conditions cannot be evaluated as bool in a Python context.")
  File "/usr/local/lib/python3.6/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.

I tried to evaluate c_s by integrating it over a computational domain by replacing

c_s = My_c_s(t, u, degree=1) 

with

c_s_help = My_c_s(t, u, degree=1)
c_s = assemble(c_s_help*dx(domain=mesh))

but this only resulted in an additional error message

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
Traceback (most recent call last):
  File "20210412_c_s_class.py", line 130, in <module>
    c_s = assemble(c_s_help*dx(domain=mesh))
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 213, in assemble
    assembler.assemble(tensor, dolfin_form)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate expression.
*** Reason:  Missing eval() function (must be overloaded).
*** Where:   This error was encountered inside Expression.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Can somebody tell me how the class My_c_s(UserExpression) should be implemented the right way?

Thanks!

Do you need any further information?

We likely need less information. Most of us don’t have time nor willpower to read long posts to track down the problem. So if you can produce a minimal working example which reproduces your error, you’re more likely to get help. See here.

1 Like

You´re right, @nate!

Basically the question is, how to convert the calculated concentration of the solution c_s to a type format, that can be evaluated in a boolean way by if c_s < c_sat?

Here is my revised MWCE

from fenics import *

################ Euler's method to calculate the differential #################
################## of the concentration of the solution c_s ###################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
   
class My_c_s(UserExpression):
    
    def __init__(self, t, u, **kwargs):
        super().__init__(**kwargs)
        self.t = t
        self.u = u   
    
    def euler_c_s(self, t0, c_s, dt, t, u):
        temp = -0

        # Iterating till the point at which we 
        # need approximation 
        while t0 < self.t: 
            temp = c_s

            c_s = c_s + assemble(dt * func_c_s(t0, c_s, self.u)) 
            t0 = t0 + dt
        return c_s

    def func_c_s(self, t, c_s, u):    
        return (self.u-c_s)*ds
       
    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)
R = FiniteElement("Real", mesh.ufl_cell(), 0)


##################### define constants #########################################
t = 0   # time
t0 = 0  # t0 used by Euler's method
T = 2.0 # total time of the simulation                    
num_steps = 10              
dt = T / num_steps
          

u_avg_domain_t0 = 1  # concentration of the corresponding domain used by My_u_t0 
u_t0 = My_u_t0()
u = interpolate(u_t0, V)

c_s_0 = 0   # concentration of c_s for t=0
c_sat = 0.9 # saturation concentration of c_s
c_s = My_c_s(t, u, degree=1)

while t < T:
    t = t+dt
    c_s.t = t

    if c_s < c_sat:
        print('t = %.1f: c_s = %.9f' % (t, c_s))
    else:
        print('Solution saturated.')

and the resulting error message

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
UFL conditions cannot be evaluated as bool in a Python context.
Traceback (most recent call last):
  File "20210426_c_s_class.py", line 80, in <module>
    if c_s < c_sat:
  File "/usr/local/lib/python3.6/dist-packages/ufl/conditional.py", line 46, in __bool__
    error("UFL conditions cannot be evaluated as bool in a Python context.")
  File "/usr/local/lib/python3.6/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.

I´m running FEniCS via Docker with the lastest stable version.

Thanks!

1 Like

c_s is a My_c_s which is an expression which is a field defined everywhere in your domain. c_sat is just a number. It doesn’t make sense to perform the boolean comparison if c_s < c_sat.

1 Like

Thanks for clarifying that, @nate!

The only reason I tried to convert the set of functions, which also calculate c_s using the Euler’s method, to a class, is because I know how to update c_s by c_s.t = t inside the for loop for every new time step.

Can somebody explain me how to update the set of functions for every time step after “calling” it initially by c_s = euler_c_s(t0, c_s_0, dt, t, u)?

That’s the revised MWCE, which only differs in using a set of functions to calculate c_s instead of the class My_c_s,


################ Euler's method to calculate the differential #################
################## of the concentration of the solution c_s ###################
# dc_s / dt = integral^{Gamma} (u - c_s) dA
   
def euler_c_s( t0, c_s, dt, t, u):
    temp = -0

    # Iterating till the point at which we 
    # need approximation 
    while t0 < t: 
        temp = c_s

        c_s = c_s + assemble(dt * func_c_s(t0, c_s, u)) 
        t0 = t0 + dt
    return c_s

def func_c_s(t, c_s, u):    
    return (u-c_s)*ds
       
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)
R = FiniteElement("Real", mesh.ufl_cell(), 0)


##################### define constants #########################################
t = 0   # time
t0 = 0  # t0 used by Euler's method
T = 2.0 # total time of the simulation                    
num_steps = 10              
dt = T / num_steps
          

u_avg_domain_t0 = 1  # concentration of the corresponding domain used by My_u_t0 
u_t0 = My_u_t0()
u = interpolate(u_t0, V)

c_s_0 = 0   # concentration of c_s for t=0
c_sat = 0.9 # saturation concentration of c_s
c_s = euler_c_s(t0, c_s_0, dt, t, u)
print('t = %.1f: c_s = %.9f' % (t, c_s))

while t < T:
    t = t+dt
    c_s.t = t
    if c_s < c_sat:
        print('t = %.1f: c_s = %.9f' % (t, c_s))
    else:
        print('Solution saturated.')

which is creating to following error message

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
t = 0.0: c_s = 0.000000000
Traceback (most recent call last):
  File "20210426_c_s_func.py", line 72, in <module>
    c_s.t = t
AttributeError: 'int' object has no attribute 't'

Thanks!