### For loop does not update the unknown

121
views
0
3 months ago by
Hello guys,

I have a system of two partial differential equations simulating the diffusion and reaction of two species 1 and 2.  I would like to check how the concentration of each changes with time in a 2D space. I have initial conditions for both functions and Neumann conditions set to zero.

The code is the following:

# Data....

T = 5
num_steps = 500
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)

#Functions
u = Function(V)
u_n = Function(V)

#Splitting of components
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)

#Constants
k = Constant(dt)
K = Constant(K)
D1 = Constant(D1)
D2 = Constant(D2)
A = Constant(A)
B = Constant(B)

#Initial conditions
g = Expression(('A+B+0.001*exp(-100*pow(x[0]-1.0/3.0,2)-100*pow(x[1]-1.0/2.0,2))','B/pow(A+B,2)'),element=ME,A=A,B=B)
u_n = project(g, V)

#Variational Problem

# Create progress bar
progress = Progress('Timestepping')
set_log_level(PROGRESS)

#Time stepping
t = 0
for n in range(num_steps):

t = t + dt
solve(F == 0, u)

u_p1, u_p2 = u.split(True)

#Update
u_n.assign(u)
#u_n.assign(project(u,V))

# Update progress bar
progress.update(t / T)

​

I think that the problem is in the for loop. When I run i get 15 or so Newton Raphson iterations, where it converges, and then a series "0 iterations". Can anyone spot where the mistake is?

Moreover, i found the line for the progress bar on the tutorial and wanted to tried it out to see what it did, but I really can't see it at the moment. Probably because the loop is not working.

Thanks

Community: FEniCS Project

2
3 months ago by
The error should be this one: the line in which you update the solutions is outside from the loop! Try to write the lines:

u_p1, u_p2 = u.split(True)

#Update
u_n.assign(u)
#u_n.assign(project(u,V))

# Update progress bar
progress.update(t / T)

In the same column as  the previous lines

t = t + dt
solve(F == 0, u),

Like this:
#Time stepping
t = 0
for n in range(num_steps):

t = t + dt
solve(F == 0, u)

u_p1, u_p2 = u.split(True)

#Update
u_n.assign(u)
#u_n.assign(project(u,V))

# Update progress bar
progress.update(t / T)
​

Okay, now i see this time bar in the terminal window, so i suppose it is stepping, but i still don't see any changes in the values on the plots.

Can i be mistaking the plotting?
def plot_current(u,t):
uv = u.compute_vertex_values(mesh)
uv = uv.reshape(nx+1,nx+1)
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​

EDIT: Okay, i think that the plotting is correct but the code is still not stepping. I obtain "num_steps" identical plots.

written 3 months ago by Giovanna Alati
I think i solved the problem, but i would still like someone to tell me why it is so.

I moved the line u_n1, u_n2 = split(u_n) from where it was before (see first post), in between the #Initial Conditions and #Variational Problem sections.
It was not splitting the conditions correctly.
written 3 months ago by Giovanna Alati
Yes! You were splitting too early. Indeed you defined u_n1 and u_n2 to be zero, as you split u_n (which becomes the a function identically zero with the definition u_n = Function(V)). Then you changed u_n projecting the function g, but u_n1 and u_n2 are then not updated (remaining zero). (I think this is the right explanation! )
written 3 months ago by Davide
Yes! You were splitting too early. Indeed you defined u_n1 and u_n2 to be zero, as you split u_n (which becomes the a function identically zero with the definition u_n = Function(V)). Then you changed u_n projecting the function g, but u_n1 and u_n2 are then not updated (remaining zero). (I think this is the right explanation! )
written 3 months ago by Davide
The results still don't make much sense...i have zero Neumanns, but at the boundaries something accumulates.
Moreover, i was playing with the initial conditions, and even when i change just a number, I receive this:

--- Instant: compiling ---
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
​