Hi guys i was trying to solve a coupled problem of A-D equations, basically the turning problem, the strong formulation is:

\begin{cases} u_t-\varDelta u=f(u,v)\\ v_t -d\varDelta v =g(u,v)\\ u(x,y,z)=u_0(x,y)\\ v(x,y,z)=v_0(x,y)\\ +\text{Periodic BC} \end{cases}

where basically the forcing terms are:

\begin{align}f(u,v)&=u-\alpha v +\gamma u v - u^3\\ g(u,v)&= u-\beta v\end{align}

and 2 initial conditions:

u_0(x,y),v_0(x,y) that are initialized using a Gaussian distribution in my code, or for simplicity by a constant. \beta, \thinspace d are scalar constants

The code is:

```
from dolfin import *
import numpy as np
# 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.
#Parameters
gamma = Constant(0.9)
alpha = 6
beta = 4
d = 20
T = 5.0 # final time
num_steps = 100 # number of time steps
dt = T / num_steps
# Create mesh and finite element
mesh = UnitSquareMesh(64, 64)
P1 = FiniteElement("P", triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element, constrained_domain=PeriodicBoundary())
# Condensed definition and encapsulation of functions
v_1, v_2 = TestFunctions(V)
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)
#initial conditions: Gaussian functions
u1init = Expression(('exp(-a*pow(x[0], 2) - a*pow(x[1], 2))', 'exp(-a*pow(x[0], 2) - a*pow(x[1], 2))'), degree=2,a=5)
u2init = Expression(('exp(-a*pow(x[0], 2) - a*pow(x[1], 2))', 'exp(-a*pow(x[0], 2) - a*pow(x[1], 2))'),degree=2,a=2)
# Forcing terms
ff = Expression('x[0] - B*x[1] + G*x[0]*x[1] - pow(x[0],3)', A=6.0, B=Constant(4.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 \
+ dot(grad(u_1), grad(v_1))*dx \
+ d*dot(grad(u_2), grad(v_2))*dx \
- (u_n1*v_1/dt + u_n2*v_2/dt)*dx \
- (ff*v_1 + gg*v_2)*dx
#Time loop:
# Create progress bar
#progress = Progress("Time-stepping")
#set_log_level(PROGRESS)
# Time-stepping
t = 0
#initial condition initializzation
u_1= interpolate(u1init, V)
u_2=interpolate(u2init, V)
#u_1,u_2=u1init,u2init
##vtkfile_u_1 = File('reaction_system/u_1.pvd')
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)
#plot(u_1)
# Update progress bar
#progress.update(t / T)
# Save solution to file
#file = File("periodic.pvd")
#file << u
plot(u_1)
# Plot solution
#plot(u_1, interactive=True)
```

The problem I encountered is that for each time I have only one iteration of the *Newton method*, I get a plot but is non-physical.

The compiler says (repeated for each time step):

```
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.125e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.178e-13 (tol = 1.000e-10) r (rel) = 2.856e-12 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
```

I mean the *Weak Formulation* is right, but how to make things work properly? I followed This Tutorial. Thanks in advance!