Problem defining initial conditions : Expression problem?

I/O
1

Hi everybody !

I am using Fenics 2018.1.0, installed with Docker, running in a Jupyter Notebook on a Windows 10 machine.

I am trying to solve simple reaction-diffusion problem and I face many problems with my code. I was used to older versions of Fenics but I just came back to the new version and it seems that things changed a bit.

I am trying to put Initial Conditions to my problem and it gives me this as Input/Output :

from dolfin import *
from numpy.random import random
import numpy
import random

class TuringPattern(NonlinearProblem):

    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self,b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

# Load mesh from file
mesh = UnitSquareMesh(48,48)

# Define function spaces (P2-P1)
U = FiniteElement('P', triangle, 1)
Wi = U * U
W = FunctionSpace(mesh,Wi)

# Define trial and test functions
du = TrialFunction(W)
(q, p) = TestFunctions(W)

# Define functions
w = Function(W)
w0 = Function(W)

# Split mixed functions
(dact, dhib) = split(du)
(act, hib) = split(w)
(act0, hib0) = split(w0)

# Set parameter values
dt = 0.05
T = 20.0

# Initial conditions
class IC(Expression):

    def eval(self,values,x):

        values[0] = Expression(fenics.Expression('1.0*random() + 0.25'))
        values[1] = Expression(fenics.Expression('1.0*random() + 0.25'))

    def value_shape(self):
        return(2,)
    
w_init = IC(element = W.ufl_element())
w.interpolate(w_init)
w0.interpolate(w_init)

L0 = act*q*dx - act0*q*dx \
+ dt*0.0005*inner(grad(act), grad(q))*dx \
- dt*inner(act*act*hib,q)*dx \
+ 1.0*dt*inner(act,q)*dx

L1 = hib*p*dx - hib0*p*dx \
+ dt*0.1*inner(grad(hib), grad(p))*dx \
+ dt*inner(act*act*hib,p)*dx \
- dt*inner(Constant(1.0),p)*dx

L = L0 + L1

a = derivative(L, w, du)

# Create files for storing solution

ufile = File("results/pattern.pvd")

# Create nonlinear problem and Newton solver
problem = TuringPattern(a, L)
solver = NewtonSolver()

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:
    print ("t =", t)
w0.vector()[:] = w.vector()
solver.solve(problem, w.vector())

# Save to file
ufile << w.split()[0]

# Move to next time step
t += dt

And I get in return :
RuntimeError Traceback (most recent call last)
in ()
51 return(2,)
52
—> 53 w_init = IC(element = W.ufl_element())
54 w.interpolate(w_init)
55 w0.interpolate(w_init)

/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py in init(self, cpp_code, *args, **kwargs)
367
368 if not isinstance(cpp_code, (str, tuple, list)):
–> 369 raise RuntimeError(“Must supply C++ code to Expression. You may want to use UserExpression”)
370 else:
371 params = kwargs

RuntimeError: Must supply C++ code to Expression. You may want to use UserExpression

I think the issue is with the definition of my InitialCondition class… But I can’t see what I am doing wrong…

Thank you for your help !

As I am new to this forum, I also want to apologize if my Syntax is not right or if my code is too long… Please tell me and I will rewrite my post.

2

Custom expressions defined in Python are no longer subclasses of Expression, but, instead, of a separate class, UserExpression. However, I will note that what it looks like you’re trying to do can be accomplished with the standard Expression constructor, by passing a tuple of C++ snippets:

w_init = Expression(2*('1.0*random() + 0.25',),
                    element = W.ufl_element())

I will also add that I’m not sure the C++ random() function is doing what you think it is; this will return very large numerical values. To get random floating point numbers between, say, 0 and 1, you would want something like

random()/((double)RAND_MAX)

because random() on its own returns an integer between 0 and RAND_MAX (which is about 2 billion).

EDIT: Alternatively, if you are trying to use the NumPy random() function (based on the import statements at the top) you could do something like

from numpy.random import random
class IC(UserExpression):
    def eval(self,values,x):
        values[0] = 1.0*random() + 0.25
        values[1] = 1.0*random() + 0.25
    def value_shape(self):
        return(2,)
w_init = IC(element = W.ufl_element())
5 Likes
3

Thanks kamenski ! Your answer was more than complete and solved my problem! Thank you very much :slight_smile: