Time dependent source function in reaction diffusion system

1

How can I make the source term time dependent inside a certain region?
I want to define region in which the concentration in changed as a function of time. At a certain time point I want to stop the time dependent source term and let the system evolve on its own.

I have example code ready to show you

T =400           # final time
num_steps = 1600# number of time steps
dt = T / num_steps # time step size


mesh = RectangleMesh( Point(0, 0), Point(10, 10),
                     70, 70, "crossed" )
# Boundary COndition
pbc = PeriodicBoundary(2)
# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)#,constrained_domain=pbc)

# Define test functions
vRa, vMa, vGa, = TestFunctions(V)

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

GT = 0.18
MT = 1.24
RT = 0.4
#initial conditions
# R active/passive
Ra_0 = Constant(0)

# M active/passive
Ma_0 = Constant(0)

#G active/passive
Ga_0 = Constant(0)


#Create initial InitialConditions
Ra_n= interpolate(Ra_0, V.sub(0).collapse())
Ma_n = interpolate(Ma_0, V.sub(1).collapse())
Ga_n = interpolate(Ga_0, V.sub(2).collapse())


assign(u_n, [Ra_n, Ma_n, Ga_n, ])


Ra, Ma, Ga = split(u)
Ra_n, Ma_n, Ga_n = split(u_n)


tol = 1e-14
f = Expression('sqrt(pow((x[0]-5),2)+pow((x[1]-5),2)) <= 1 + tol ? 10*t : 0', degree=0,t=0,tol=tol)

fGa = project(f, V.sub(2).collapse())


# Define expressions used in variational forms k2 = 2.04*0.0745, k6 = 0.00509*0.784
k1,k2,k3,k4,k5,k6 = Constant(0),Constant(0),Constant(0),Constant(0),Constant(0),Constant(0)#Constant(3.88*2.42),Constant(2.04*0.0745),Constant(1.19),Constant(3.98),Constant(0.4715*0.014),Constant(0.00509*0.784) #Constant(0),Constant(0),Constant(0),Constant(0),Constant(0),Constant(0)
Km1,Km2,Km5,Km6 = Constant(2.42),Constant(0.0745),Constant(0.014),Constant(0.784)
DRa,DMa,DGa = Constant(0.28),Constant(0.03),Constant(0.028)
k = Constant(dt)

# Create VTK files for visualization output
vtkfile_Ra = File('reaction_system/Ra.pvd')

vtkfile_Ma = File('reaction_system/Ma.pvd')

vtkfile_Ga = File('reaction_system/Ga.pvd')


# define NonlinearProblem
class RDSystem(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)


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

    t += dt


    if t > 1 and t < 100:
        f.t =  t
        fGa = project(f, V.sub(2).collapse())

    else:
        f.t =  0
        fGa = project(f, V.sub(2).collapse())



    # Solve variational problem for time step
    # Define variational problem
    FRa = ((Ra-Ra_n)/k)*vRa*dx  - k1*Ga*(RT-Ra)/(Km1+RT-Ra)*vRa*dx  + k2*Ra/(Km2+Ra)*vRa*dx + DRa*dot(grad(Ra),grad(vRa))*dx
    FMa = ((Ma-Ma_n)/k)*vMa*dx - k5*Ra*(MT-Ma)/(Km5+(MT-Ma))*vMa*dx + k6*Ma/(Km6+Ma)*vMa*dx +DMa*dot(grad(Ma),grad(vMa))*dx
    FGa = ((Ga-Ga_n)/k)*vGa*dx - k3*(GT-Ga)*Ra*vGa*dx + k4*Ga*(MT-Ma)*vGa*dx + DGa*dot(grad(Ga),grad(vGa))*dx -fGa*vGa*dx

    L = FRa  + FMa + FGa
    a = derivative(L, u, du)
    #
    problem = RDSystem(a, L)
    solver = NewtonSolver()
    solver.parameters["linear_solver"] = "lu"
    solver.parameters["convergence_criterion"] = "incremental"
    solver.parameters["relative_tolerance"] = 1e-6

    solver.solve(problem, u.vector())

    # Save solution to file (VTK)
    _Ra, _Ma, _Ga = u.split()

    vtkfile_Ra << (_Ra, t)

    vtkfile_Ma << (_Ma, t)

    vtkfile_Ga << (_Ga, t)

    # Update previous solution
    u_n.assign(u)
    # Update current time

Best,
Dennis

2

Please format your code properly, and remove all code that is not needed, instead of commenting it out. Please try to reduce the code to as small of an example as possible. This makes it easier for people to give you suggestions.