Homogeneous Neumann not applied

3 months ago by

Hello guys!
I don't understand why my code does not implement the homogeneous Neumann condition that comes from the following variational form, which is correct:

Here u1 and u2 are concentrations of two species which are reacting and diffusing within a square domain. From the plot i see, as it time-steps, that the concentration of one of the species increases on the left boundary of the domain. That shouldn't happen.
I don't understand if i have to specify that the Neumann boundaries are actually the boundaries of the domain, even though in all examples-code I checked it was not necessary.
Please let me know! Thank you!

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

# Data

mesh = UnitSquareMesh(nx,ny)

#Function space
P2 = FiniteElement('P',triangle,2)
ME = MixedElement([P2, P2])
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)


#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)
u_n1, u_n2 = split(u_n)

#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')

#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 progress bar
    progress.update(t / T)
Community: FEniCS Project

1 Answer

3 months ago by
Hi there,
pure Neumann boundary might be the reason. Try to test the problem by applying a Dirichlet boundary in order to understand if the code is working as expected.
Best, Emek
Hi! I am followed you suggestion and added some simple DBC, which I added in between the space function and the Test functions:
#Dirichlet Boundary 

def boundary_D(x,on_boundary):
    if on_boundary:
        if near(x[0],0,tol)or near(x[0],1,tol):
            return True
            return False
        return False
u_d = Constant(5)
bc = DirichletBC(V.sub(0),u_d,boundary_D)​

...and nothing happens. The plot stays exactly the same. I really don't know what to try more, the problem itself should be very simple to implement...

written 3 months ago by Giovanna Alati  
The constants D1 and D2 are very very sensitive. Check them out carefully and try to play with them as well. Moreover, the whole problem is difficult to solve, especially if you have very small numbers. You may solve the problem with 0.1 as initial concentration but not with 0.0001.
written 3 months ago by Emek  
This should not happen. Afaik Dirichlet BC are hardly enforced. Did you forget to give your solve command the boundary condition?
written 3 months ago by Lukas O.  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »