Homogeneous Neumann not applied
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 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) #Functions u = Function(V) u_n = Function(V) #Splitting of components u_1, u_2 = split(u) #Constants ... #Initial conditions g = Expression(('A+B+0.001*exp(-100*pow(x-1.0/3.0,2)-100*pow(x-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') 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) # Update progress bar progress.update(t / T)
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.