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 \
- (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!

There is definitely something wrong in the assignment of the initial conditions.

This should do what you are expecting

#initial condition initializzation
u1init_interp = interpolate(u1init, V)
u2init_interp = interpolate(u2init, V)
u.sub(0).assign(u1init_interp)
u.sub(1).assign(u2init_interp)


Note that I avoided using u_1 and u_2 on purpose, since they have already been used when you wrote split(u).

Output of the first 3 time steps is

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) = 2.475e-14 (tol = 1.000e-10) r (rel) = 5.999e-13 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
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) = 3.480e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.582e-14 (tol = 1.000e-10) r (rel) = 7.420e-13 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
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) = 3.408e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.737e-14 (tol = 1.000e-10) r (rel) = 8.032e-13 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.


i.e. the residual at the initial guess now changes from one time to the next. I havenâ€™t plotted the solution, since I have no intuition of what solution to expect.

Finally, for future posts please make an effort to clean up the code. There are a lot of commented lines that are just noise that a person who is willing to help you have to move through.