Adding a random perturbation to a solution at each time step

Question

My code solves a system of reaction-diffusion equations on a spherical shape. The system is solved at each time step, and I wish to add a small random perturbation to each solution as they will become the initial conditions for the next time step integration. The one caveat is that it could be the case that only a subset of the solution needs to be updated, e.g., the elements where the z > 0.5 on a unit sphere.

I suspect there must be a simple way of doing this, but I haven't been able to figure it out.

The guts of the code roughly look like this:

while t <= T:
b_u1 = assemble(L_u1, tensor = b_u1)
bc_u1.apply(A_u1, b_u1)
solve(A_u1, u1.vector(), b_u1)
t += dt
u1_old.assign(u1)

Answer 1

Hi,

I guess there are several ways of doing this. The Cahn-Hilliard demo contains an example of how to generate initial conditions using random values in an Expression. A simple solution to your problem (neglecting caveat) could be:

from numpy import random
eps = 0.001 # Weight of random number
while t <= T:
    b_u1 = assemble(L_u1, tensor = b_u1)
    bc_u1.apply(A_u1, b_u1)
    solve(A_u1, u1.vector(), b_u1)
    t += dt
    u1_old.vector().set_local(u1.vector().array()+eps*(0.5-random.random(u1.vector().size())))

If you only want to add pertubations to some part of your problem you can, e.g., mark this using a SubDomain or use an Expression that evaluates only in that domain, like

class OuterPart(Expression):
    def eval(self, values, x):
        if sqrt(x[0]**2+x[1]**2+x[2]**2) > 0.5:
            values[0] = 0.01*(0.5 - random.random())

You would have to interpolate this Expression onto the function space of u1 each timestep and then add it, like

pert = interpolate(OuterPart(), V)
u1_old.assign(u1)
u1_old.vector().axpy(1.0, pert.vector())

There may be more efficient ways. Perhaps this is better:

radius = interpolate(Expression("sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])"), V)
num_dofs_perturbed = radius.vector()[radius.vector()>0.5].size()
u1_old.assign(u1)
u1_old.vector()[radius.vector()>0.5] += 0.01*(0.5-random.random(num_dofs_perturbed))

where the first two lines are precomputed and only the last two lines go inside your while loop.

Comments

I altered the top suggestion to

u2_old.vector().set_local(u2.vector().array()+0.01*randn())

and it works great, thank you!