### Can anyone understand why the Newton iteration does not proceed?

294

views

0

Hello guys again, I am really sorry for my indelicate behavior.

Here i have reformulated the problem. I am practically reproducing the case of advection-diffusion that can be found in the Fenics tutorial with some little changes. ( Only 2 species, no source function, imposed initial conditions different from zero).

As for the platform i am using, it is Jupiter Notebook.

I keep on getting the same error:

"Newton solver finished in 0 iterations and 0 linear solver iterations.

No Jacobian form specified for nonlinear variational problem.

Differentiating residual form F to obtain Jacobian J = F'."

Let me know, thank you.

UPDATE: I think that the iteration is not taking place because it is not taking the right initial condition. The plot that i get with imshow does not coincide with the function written in g.

Here i have reformulated the problem. I am practically reproducing the case of advection-diffusion that can be found in the Fenics tutorial with some little changes. ( Only 2 species, no source function, imposed initial conditions different from zero).

As for the platform i am using, it is Jupiter Notebook.

I keep on getting the same error:

"Newton solver finished in 0 iterations and 0 linear solver iterations.

No Jacobian form specified for nonlinear variational problem.

Differentiating residual form F to obtain Jacobian J = F'."

Let me know, thank you.

UPDATE: I think that the iteration is not taking place because it is not taking the right initial condition. The plot that i get with imshow does not coincide with the function written in g.

```
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
# Data
D1 = 0.05
D2 = 1
K = 100
A = 0.1305
B = 0.7695
T = 3.0
num_steps = 10000
dt = T/num_steps
nx = 100
ny = 100
#Mesh
mesh = UnitSquareMesh(nx,ny)
#Function space
P1 = FiniteElement('P',triangle,1)
ME = MixedElement([P1, P1])
V = FunctionSpace(mesh,ME)
#Test Functions
v_1, v_2 = TestFunctions(V)
u = Function(V)
u_n = Function(V)
#Splitting of components
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)
#Initial conditions
g= Expression(('A+B+pow(10,-3)*exp(-100*pow((x[0]-1/3),2)-100*pow((x[1]-1/2),2))','B/pow(A+B,2)'),element=ME,A=A,B=B)
u_n = project(g,V)
#Constants
dt = Constant(dt)
K = Constant(K)
D1 = Constant(D1)
D2 = Constant(D2)
A = Constant(A)
B = Constant(B)
#Variational Problem
F = ((u_1 -u_n1)/dt)*v_1*dx + D1*dot(grad(u_1),grad(v_1))*dx + ((u_2 -u_n2)/dt)*v_2*dx + D2*dot(grad(u_2),grad(v_2))*dx - K*(A-u_1+((u_1)**2)*u_2)*v_1*dx-K*(B-((u_1)**2)*u_2)*v_2*dx
# Create progress bar
progress = Progress('Timestepping')
set_log_level(PROGRESS)
#Time stepping
t=0
u_n1, u_n2 = u_n.split(True) #n
for n in range(num_steps):
t = t + dt
solve(F == 0, u)
#Splitting
u_1, u_2 = u.split(True) #n+1
#Update
u_n.assign(u)
# Update progress bar
progress.update(t / T)
```

Community: FEniCS Project

### 1 Answer

0

Hello Martina Cotti

Errors in your code is as follows:-

time T should be defined as T=3.0

If u define T=3

dt = T/num_steps

dt will become Zero.

so define T as float, not as integer,

Thanks Try this code below , its Working

Errors in your code is as follows:-

time T should be defined as T=3.0

If u define T=3

dt = T/num_steps

dt will become Zero.

so define T as float, not as integer,

Thanks Try this code below , its Working

```
from dolfin import *
from fenics import *
import matplotlib.pyplot as plt
D1 = 0.05
D2 = 1
K = 100
A = 0.1305
B = 0.7695
T = 3.0
num_steps = 10000
dt = T/num_steps
nx = 100
ny = 100
#print dt
#raise SystemExit
#Mesh
mesh = UnitSquareMesh(nx,ny)
#Function space
P1 = FiniteElement('P',triangle,1)
ME = MixedElement([P1, P1])
V = FunctionSpace(mesh,ME)
#Test Functions
v_1, v_2 = TestFunctions(V)
u = Function(V)
u_n = Function(V)
#Splitting of components
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)
#Initial conditions
#g1 = Expression('A+B+pow(10,-3)*exp(-100*pow((x[0]-1/3),2)-100*pow((x[1]-1/2),2))', A=A,B=B, degree=1)
#g2 = Expression('B/pow(A+B,2)', A=A,B=B, degree=1)
g = Expression(('A+B+pow(10,-3)*exp(-100*pow((x[0]-1/3),2)-100*pow((x[1]-1/2),2))','B/pow(A+B,2)'),element=ME,A=A,B=B)
u_n = project(g,V)
#Constants
dt = Constant(dt)
K = Constant(K)
D1 = Constant(D1)
D2 = Constant(D2)
A = Constant(A)
B = Constant(B)
#Variational Problem
F = ((u_1 -u_n1)/dt)*v_1*dx + D1*dot(grad(u_1),grad(v_1))*dx + ((u_2 -u_n2)/dt)*v_2*dx + D2*dot(grad(u_2),grad(v_2))*dx - K*(A-u_1+((u_1)**2)*u_2)*v_1*dx-K*(B-((u_1)**2)*u_2)*v_2*dx
#F = -dt*K*u_1**2*u_2*v_1*dx + (1+dt*K)*u_1*v_1*dx + dt*D1*dot(grad(u_1),grad(v_1))*dx + \
# + (1+dt*K*u_1**2)*u_2*v_2*dx + dt*D2*dot(grad(u_2),grad(v_2))*dx + \
# - (u_n1*v_1+dt*K*A*v_1)*dx - (u_n2*v_2+dt*K*B*v_2)*dx
# Create progress bar
progress = Progress('Timestepping')
set_log_level(PROGRESS)
#plot
plt.ion()
plt.clf()
xv = mesh.coordinates()
X = xv[:,0].reshape(nx+1,nx+1)
Y = xv[:,1].reshape(ny+1,ny+1)
#plot function
def plot_current(u,t):
uv = u.compute_vertex_values(mesh)
uv = uv.reshape(nx+1,nx+1)
#zmax = np.max(np.abs(v.compute_vertex_values(mesh)))
#ax.cla()
plt.imshow(uv,cmap = 'jet',extent=[0,1,1,0],interpolation = 'none')
plt.gca().invert_yaxis()
plt.colorbar()
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.title('u at $t=$'+str(t))
plt.pause(0.05)
return
#Time stepping
t=0
u_n1, u_n2 = u_n.split(True) #n
fig1 = plt.figure(1)
plot_current(u_n1,t)
fig1.savefig('Problem3_u_n1.pdf')
fig2 = plt.figure(2)
plot_current(u_n2,t)
fig2.savefig('Problem3_u_n2.pdf')
for n in range(num_steps):
t = t + dt
solve(F == 0, u)
#Splitting
u_1, u_2 = u.split(True) #n+1
#plot at t + dt
#fig3 = plt.figure(3)
#plot_current(u_1,t)
#fig3.savefig('Problem3_u_1.pdf')
#fig4 = plt.figure(4)
#plot_current(u_2,t)
#fig4.savefig('Problem3_u_2.pdf')
#Update
u_n.assign(u)
# Update progress bar
progress.update(t / T)
fig5 = plt.figure(5)
plot_current(u_1,2)
```

Hello Md Ayub Ansari,

I still get the same error as above unfortunately :(

I still get the same error as above unfortunately :(

written
6 months ago by
Martina Cotti

Why you have not specified any boundary condition? neither Dritchlet nor Newman, is it not required?

written
6 months ago by
Md Ayub Ansari

Yeah, i didn’t mention that there are homogeneous Neumann BC because the integral only appears in the weak formulation and it is 0

EDIT: for a moment i thought that i had forgotten to specify the boundaries, but they are actually not present in the fenics tutorial code as well

EDIT: for a moment i thought that i had forgotten to specify the boundaries, but they are actually not present in the fenics tutorial code as well

written
6 months ago by
Martina Cotti

Please login to add an answer/comment or follow this question.

Neither did you explain what you want to do, what the problem is you're solving or any hint at the tutorial you just mentioned in your comment.

Plus you didn't even use code formatting for your rather long code so even just copy-pasting your example would need manually fixing the indentations to render the code executable.

I guess your chances of getting an answer are slim.

I think in the way you implemented it does not affect the form and therefore the solution is simply 0!

"Newton solver finished in 0 iterations and 0 linear solver iterations."

If it is not iterating, which solution am i looking at on the plots?

> UPDATE: I think that the iteration is not taking place because it is not taking the right initial condition. The plot that i get with imshow does not coincide with the function written in g.

One thing to consider is that the string passed to the Expression constructor is C++ code. In C++, "1/3" evaluates to zero, because "1" and "3" are interpreted as ints, and integer division and throws away the remainder. (Python 2.x does this as well, but in Python 3.x, there is a separate integer division operator, "//".) You can get around this by using floating-point literals in fractions instead of integer literals, e.g., "1.0/3.0".