### Homogeneous Neumann BC not fulfilled?

350
views
1
8 months ago by
In my opinion, in order to achieve homogeneous Neumann boundary conditions in most cases the "do nothing" approach is fine. However, I found out that in my code the boundary conditions are 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] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 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
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')

​

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?

Community: FEniCS Project
Maybe I'm missing something; but it looks like you haven't set any boundary conditions at all. Is this a well posed problem?
written 8 months ago by Alexander G. Zimmerman
This is a slightly modified version of the standard Cahn-Hilliard Demo
https://fenicsproject.org/olddocs/dolfin/2016.2.0/python/demo/documented/cahn-hilliard/python/demo_cahn-hilliard.py.html

I can't tell you by heart but I think with the boundary 0-Neumann for  $c$c  and  $\mu$μ this should be well-posed. There is some work by Elliott and Songmu on this.
written 8 months ago by Lukas O.

1
8 months ago by

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.

Unfortunately the problem persists for every theta value I tried. [including 1 and 0.51 ;) ]
written 8 months ago by Lukas O.
Do the boundary values oscillate in time? If so, then this still may be a related issue. If not, then my hunch is probably unrelated.
written 8 months ago by Alexander G. Zimmerman
I can't find any oscillatory behavior. If you have 0-Neumann in your FEniCS code, did you check if they hold?
written 8 months ago by Lukas O.
1
Yes, I have a couple of tests with adiabatic walls (homogeneous Neumann BC on the temperature equation). Here's a snapshot of the output from [my natural convection of air test](https://github.com/geo-fluid-dynamics/phaseflow-fenics/blob/master/tests/test_wang2010.py) visualized in ParaView:
written 8 months ago by Alexander G. Zimmerman