Initial condition for nonlinear problem

1

Im trying to implement a nonlinear problem with two species like the example problem of reaction-difussion-advection but for Schnakenberg model. I have to set a random initial condition in a specific range for both variables, nevertheless i don’t know how to implement this.

from __future__ import print_function
from fenics import *

T = 10.0            # final time
num_steps = 2000    # number of time steps
dt = T / num_steps # time step size
d = 8.6676        # diffusion coefficient
gamma = 230.82 
a = 0.1
b = 0.9

# Create mesh and define function space
nx = ny = 50
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_1, v_2 = TestFunctions(V)

# Define functions for concentrations
u = Function(V)
u_n = Function(V)

# Split system functions to access components
u_1, u_2= split(u)
u_n1, u_n2 = split(u_n)

# Create intial conditions and interpolate
u_init = InitialConditions()
u.interpolate(u_init)
un.interpolate(u_init)

# Define expressions used in variational forms
k = Constant(dt)
d1 = Constant(d)
d2 = Constant(d)
a = Constant(a)
b = Constant(b)

# Define variational problem
F = ((u_1 - u_n1) / k)*v_1*dx + d1*dot(grad(u_1), grad(v_1))*dx\
   -(gamma*(u_1*u_1*u_2-u_1+a))*v_1*dx \
  + ((u_2 - u_n2) / k)*v_2*dx + d2*dot(grad(u_2), grad(v_2))*dx\
  -(gamma*(-u_1*u_1*u_2+b))*v_2*dx 
  
# Create time series for reading velocity data

# Create VTK files for visualization output
vtkfile_u_1 = File('reaction_system/u_1.pvd')
vtkfile_u_2 = File('reaction_system/u_2.pvd')

# Create progress bar
progress = Progress('Time-stepping')
#conda set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Solve variational problem for time step
    solve(F == 0, u)

    # Save solution to file (VTK)
    _u_1, _u_2 = u.split()
    vtkfile_u_1 << (_u_1, t)
    vtkfile_u_2 << (_u_2, t)

    # Update previous solution
    u_n.assign(u)

    # Update progress bar
    #progress.update(t / T)

# Hold plot

the initial value for u_1 and u_2 must be random in a range, do you know how to implement that? Please

2

Take a look at the Cahn-Hilliard demo.

1 Like
3

Something like this should work

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        np.random.seed(seed=0)
        super().__init__(**kwargs)
    def eval(self, values, x):
        values[0] = 2.e0*np.random.rand() - 1.e0
        values[1] = 2.e0*np.random.rand() - 1.e0
    def value_shape(self):
        return (2,)

u, u_n, v = fn.Function(W), fn.Function(W), fn.TestFunction(W)
u0 = InitialConditions(degree=1)
u.interpolate(u0)
u_n.interpolate(u0)
1 Like
4

Thanks, I am very grateful