### For loop does not update the unknown

121

views

0

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:

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

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:

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! )

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)
```

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.

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

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:

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.

Can i be mistaking the plotting?

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