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

294
views
0
6 months ago by
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.
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
I mean, i am really following the tutorial to the letter...
written 6 months ago by Martina Cotti
...and why did i get a -1? What is wrong in my question?
written 6 months ago by Martina Cotti
1
You did not put any effort in reducing your code into a minimal working example that exhibits your problem. The plots for example are unnecessary.
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.
written 6 months ago by klunkean
Are you sure it is an error? The residual is very small always. I think the initial guess might be the solution.
written 6 months ago by Eleonora Piersanti
Also,you should use this syntax!
u_n.assgn(project(g,V))

​

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

written 6 months ago by Eleonora Piersanti
Hey Eleonora! Indeed with this command i get a very different plot, quite similar to my initial conditions. But i keep on getting the same message.
"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?
written 6 months ago by Martina Cotti

> 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".

written 6 months ago by David Kamensky
Thank you, that indeed solved the problem. But it is still not iterating...
written 6 months ago by Martina Cotti

### 1 Answer

0
6 months ago by
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
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 :(
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
written 6 months ago by Martina Cotti
Please login to add an answer/comment or follow this question.

Similar posts:
Search »