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!