Generating a Random Potential

variational formulation
1

I’m trying to solve the following equation (-\nabla ^2+ V)u = 1 , where V is a random potential. This random potential takes different random values uniformly, say from 0 to N on different nodes.
I can solve the problem when I have a simple expression for V = x^2 + y^2 as follows

mesh = UnitSquareMesh(8,8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('0', degree=0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)
potential = Expression('x[0]*x[0] + x[1]*x[1]', element = V.ufl_element())

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(1.0)
a = dot(grad(u), grad(v))*dx
a = (inner(grad(u), grad(v)) \
     + potential*u*v)*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

How should I go about creating a subclass for Expression that can implement a random potential on my mesh?

1 Like
2

Hi @abelthayil,
you might have a look at the Cahn-Hilliard demo (https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/cahn-hilliard/python/documentation.html). There the initial condition is implemented by a subclass of Expression with random values.

3

Thank you for your reply. When I try to implement a sub class inspired by the example shared above as follows:

# Class representing the random potential
class RandomPotential(Expression):
    def __init__(self):
        random.seed(2)
    def eval(self, values, x):
        values[0] = random.randint(5)
    def value_shape(self):
        return (0,)

I get the following error when I call potential = RandomPotential()

RecursionError: maximum recursion depth exceeded in comparison

4

Hi @abelthayil,
I can not reproduce your error, but the line
values[0] = random.randint(5)
seems suspicious to me. I think you need to give two arguments to random.randint(a, b)… And shouldn’t the function value_shape return (1,)?

5

Hello @causem,

Sorry, the error wasn’t coming from potential = RandomPotential() but instead when I try to evaluate

        a = (inner(grad(u), grad(v)) \
            + potential*u*v)*dx

Here is the complete error statement:

RecursionError                            Traceback (most recent call last)
<ipython-input-57-5c7e016b201e> in <module>
      8 f = Constant(1.0)
      9 a = (inner(grad(u), grad(v)) \
---> 10      + potential*u*v)*dx
     11 # L = f*v*dx
     12 

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ufl/exproperators.py in _mul(self, o)
    191         return NotImplemented
    192     o = as_ufl(o)
--> 193     return _mult(self, o)
    194 
    195 

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ufl/exproperators.py in _mult(a, b)
    122     # - matrix-matrix (A*B, M*grad(u)) => A . B
    123     # - matrix-vector (A*v) => A . v
--> 124     s1, s2 = a.ufl_shape, b.ufl_shape
    125     r1, r2 = len(s1), len(s2)
    126 

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ufl/coefficient.py in ufl_shape(self)
     71     def ufl_shape(self):
     72         "Return the associated UFL shape."
---> 73         return self._ufl_shape
     74 
     75     def ufl_function_space(self):

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/function/expression.py in __getattr__(self, name)
    430             return self._parameters
    431         else:
--> 432             return self._parameters[name]
    433 
    434     def __setattr__(self, name, value):

... last 1 frames repeated, from the frame below ...

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/function/expression.py in __getattr__(self, name)
    430             return self._parameters
    431         else:
--> 432             return self._parameters[name]
    433 
    434     def __setattr__(self, name, value):

RecursionError: maximum recursion depth exceeded in comparison

To address your concerns,
randint defaults a to 0 if there is only one argument. I’ve tried with shape (1,) and the error persists so I doubt that’s the issue. Thanks again!

6

Hi @abelthayil,
I think there has been a change in the way of defining custom Expressions, the Parent class now is UserExpression (instead of Expression). So the RandomPotential would then be:

class RandomPotential(UserExpression):

    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world) )
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = random.randint(0,5)
    def value_shape(self):
        return ()

Sorry for pointing you to the outdated example!

7

Thank you, this seems to work now!