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
F = ((u_1 -u_n1)/k)*v_1*dx + D1*dot(grad(u_1),grad(v_1))*dx + ((u_2 -u_n2)/k)*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
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

1 Answer


2
3 months ago by
Davide  
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)
<ipython-input-16-0db79a191adb> in <module>()
     60 
     61 #Initial conditions
---> 62 g = Expression(('1.0*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)
     63 u_n = project(g, V)
     64     #u_n.assign(project(g,V))
​
written 3 months ago by Giovanna Alati  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »