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)