Homogeneous Neumann BC not fulfilled?
So I checked the same in the Cahn-Hilliard Demo and it seems that they are not fulfilled there either.
from __future__ import print_function import random from dolfin import * set_log_level(WARNING) # Class representing the intial conditions class InitialConditions(Expression): def __init__(self, **kwargs): random.seed(2 + MPI.rank(mpi_comm_world())) def eval(self, values, x): values = 0.63 + 0.02*(0.5 - random.random()) values = 0.0 def value_shape(self): return (2,) # Model parameters lmbda = 1.0e-02 # surface parameter dt = 5.0e-06 # time step theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson # Form compiler options parameters["form_compiler"]["optimize"] = True parameters["form_compiler"]["cpp_optimize"] = True # Create mesh and build function space mesh = UnitSquareMesh(48, 48) P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1) ME = FunctionSpace(mesh, P1*P1) n = FacetNormal(mesh) # Define trial and test functions du = TrialFunction(ME) q, v = TestFunctions(ME) # Define functions u = Function(ME) # current solution u0 = Function(ME) # solution from previous converged step # Split mixed functions dc, dmu = split(du) c, mu = split(u) c0, mu0 = split(u0) # Create intial conditions and interpolate u_init = InitialConditions(degree=1) u.interpolate(u_init) u0.interpolate(u_init) # Compute the chemical potential df/dc c = variable(c) f = 100*c**2*(1-c)**2 dfdc = diff(f, c) # mu_(n+theta) mu_mid = (1.0-theta)*mu0 + theta*mu # Weak statement of the equations L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx L = L0 + L1 # Compute directional derivative about u in the direction of du (Jacobian) a = derivative(L, u, du) # Create nonlinear problem and Newton solver problem = NonlinearVariationalProblem(L, u, J=a) solver = NonlinearVariationalSolver(problem) # Output file file = File("output.pvd", "compressed") # Step in time t = 0.0 T = 10*dt while (t < T): t += dt u0.assign(u) solver.solve() c, mu = u.split() c.rename('c','phase') print(str(assemble(inner(grad(c),n)*ds))) file << (project(grad(c)), t)
I expected that the output should be something close to zero, which is not the case (sth. up to -5.0). If I write the projected gradient into the output file and inspect it in paraview, I would expect that e.g. on the lower boundary the second component should be close to zero, which is also not the case. I also tried different polynomials orders without success.
I am aware that there are errors involved in projecting the gradient and that the BC are only weakly enforced and hence maybe not that close to zero, however I think values in the order of 1 should not be there.
Did I check the conditions not correctly? Is something wrong with the implementation?
I've run into this problem (unrelated to FEniCS, actually with a code I implemented with another library and also with something written from scratch) when using the Crank-Nicolson (exactly theta = 0.5) time discretization scheme with Neumann boundary conditions. Off the top of my head, I don't remember what my conclusion was; so I'll have to dig up some notes. For now, what happens when you set theta = 0.51? Also as long as we're at it, maybe try a fully implicit run.