Update "UserExpression" at previous time step, expression not depending on time

General
1

hi guys, i was wondering how to update those functions:
\begin{cases} f(u,v)=u-\alpha v +\gamma uv -u^3\\ g(u,v)=u- \beta v \end{cases}
where \alpha,\beta,\gamma are scalar constants.

f(u,v),\thinspace g(u,v) do not directly depend on time but i would like to update them at the previous time step. From what I got defining functions like:

ff = Expression('x[0] - A*x[1] + G*x[0]*x[1] - pow(x[0],3)', A=Constant(6.0), G=Constant(0.9),degree=1)
gg = Expression('x[0] - B*x[1]', B=Constant(4.0),degree=1)

is kind of update x,y that is space, but I would like to update u,v that are in my case the fields that I’m solving for! Therefore I tried to define a class in order to overcome the problem:

# Define UserExpression for ff and gg
class ForcingTerms(UserExpression):
    def __init__(self, u_n1, u_n2, **kwargs):
        super().__init__(kwargs)
        self.u_n1 = u_n1
        self.u_n2 = u_n2

    def eval(self, values, x):
        values[0] = self.u_n1(x[0], x[1]) - alpha * self.u_n2(x[0], x[1]) + gamma * self.u_n1(x[0], x[1]) * self.u_n2(x[0], x[1]) - self.u_n1(x[0], x[1]) ** 3
        values[1] = self.u_n1(x[0], x[1]) - beta * self.u_n2(x[0], x[1])

    def value_shape(self):
        return (2,)

and then instantiate the class using:

ff_gg = ForcingTerms(u_n1, u_n2)

while in the time-loop I update using:

while (t<T):
# Update current time
    t += dt
    #u_n.assign(u)
    u_n.assign(u)
    solve(F==0, u,J=J)

But is giving me an error related to the expression ff_gg

here the complete code:

from dolfin import *
import numpy as np
from numpy.random import random


# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
        return bool((near(x[0], 0) or near(x[1], 0)) and \
                (not ((near(x[0], 0) and near(x[1], 1)) or \
                        (near(x[0], 1) and near(x[1], 0)))) and on_boundary)





    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        if near(x[0], 1) and near(x[1], 1):
            y[0] = x[0] - 1.
            y[1] = x[1] - 1.
        elif near(x[0], 1):
            y[0] = x[0] - 1.
            y[1] = x[1]
        else:   # near(x[1], 1)
            y[0] = x[0]
            y[1] = x[1] - 1.

#Initial condtions (Random)
class IC(UserExpression):
    def eval(self,values,x):
        values[0] = 0.1*random() -0.1*random()
        values[1] = 0.1*random() -0.1*random()
    def value_shape(self):
        return(2,)
    

# Define UserExpression for ff and gg
class ForcingTerms(UserExpression):
    def __init__(self, u_n1, u_n2, **kwargs):
        super().__init__(kwargs)
        self.u_n1 = u_n1
        self.u_n2 = u_n2

    def eval(self, values, x):
        values[0] = self.u_n1(x[0], x[1]) - alpha * self.u_n2(x[0], x[1]) + gamma * self.u_n1(x[0], x[1]) * self.u_n2(x[0], x[1]) - self.u_n1(x[0], x[1]) ** 3
        values[1] = self.u_n1(x[0], x[1]) - beta * self.u_n2(x[0], x[1])

    def value_shape(self):
        return (2,)



# Class for interfacing with the Newton solver
class TuringInstability(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)



#Parameters
gamma = 0.9
alpha = 6
beta = 4
d = 20
dt = 5.0e-06

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
#parameters["form_compiler"]["representation"] = "quadrature"

# Create mesh and finite element
mesh = UnitSquareMesh(60, 60)
P1 = FiniteElement("P", triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element, constrained_domain=PeriodicBoundary())
plot(mesh)

# Condensed definition and encapsulation of functions
du       =TrialFunction(V)
v_1, v_2 = TestFunctions(V)
u = Function(V)
u_n = Function(V)

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

#Initializzatioin
u1init = IC(element = u.ufl_element())
u.interpolate(u1init)
u_n.interpolate(u1init)



# Create instance of the ForcingTerms class
ff_gg = ForcingTerms(u_n1, u_n2)
# Forcing terms
"""ff = Expression('x[0] - A*x[1] + G*x[0]*x[1] - pow(x[0],3)', A=Constant(6.0), G=Constant(0.9),degree=1)
gg = Expression('x[0] - B*x[1]', B=Constant(4.0),degree=1)"""



# Define variational problem
F  =  (u_1*v_1/dt)*dx \
    + (u_2*v_2/dt)*dx \
    + inner(grad(u_1), grad(v_1))*dx \
    + d*inner(grad(u_2), grad(v_2))*dx \
    - (dot(u_n1,v_1)/dt + dot(u_n2,v_2)/dt)*dx \
    - (dot(ff_gg[0],v_1) + dot(ff_gg[1],v_2))*dx 



#Time parameters:
t = 0.0
T=600*dt
#vtk file
vtkfile_u_1 = File('reaction_system/u_1.pvd')
vtkfile_u_2 = File('reaction_system/u_2.pvd')

#Jacobian, gateau derivative direction (du)
J=derivative(F,u,du)


while (t<T):
# Update current time
    t += dt
    #u_n.assign(u)
    u_n.assign(u)
    solve(F==0, u,J=J)
    
# 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
    
    plot(u_1)