Initial and Dirichlet Condition on Chemical Potential for the CH equation

Hi everyone, I am looking at the demo example
https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/cahn-hilliard/python/documentation.html
and noticed that the chemical potential is set to 0 there.

 class InitialConditions(Expression):
     def __init__(self):
         random.seed(2 + MPI.process_number())
     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,)

That seems odd to me because that would be inconsistent with a non constant concentration field, but doesn’t cause any error. Am I missing something? Would it be possible to initialize the chemical potential as

 class InitialConditions(Expression):
     def __init__(self):
         random.seed(2 + MPI.process_number())
     def eval(self, values, x):
         values[0] = 0.63 + 0.02*(0.5 - random.random())
         values[1] = (values[0]-0.9)*(values[0]-0.1)*(values[0]-0.5)-lmbda*div(grad(values[0]))
     def value_shape(self):
         return (2,)

Similarly, is it possible to set the chemical potential with Dirichlet boundary conditions explicitly as a function of the concentration field such as (c-0.9)*(c-0.1)*(c-0.5)-lmbda*div(grad(c)?

You are looking at an ancient version of the demos. Consider: https://fenicsproject.org/olddocs/dolfin/latest/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html

Not with a user expression. You would have to first set c, the project this initial condition into mu.

I would probably look at: A conservative numerical method for the Cahn–Hilliard equation with Dirichlet boundary conditions in complex domains - ScienceDirect and references therein

Thanks for the quick reply! How come it doesn’t seem to matter what initial condition we set the chemical potential to be though? ( Even in the latest demo the chemical potential initial condition is still set to 0)

I am actually not sure. It is quite common in mixed problems that one only has a valid initial state for one variable (for instance navier stokes where the initial condition might just be an inlet condition).

For a problem with an analytical solution this usually throws of the convergence rates a bit.

You could add a projection that brings the initial conditions into a feasible space, and it should make the solution slightly more accurate.

This is a bit outside of my field of expertise.

Thanks! There seems to be quite a strange thing going on here. An even stranger thing though is that for a given initial condition such as

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(MPI.comm_world))
        super().__init__(**kwargs)
    def eval(self, values, x):
        L = pi
        base_val = 0.9
        amp = -0.8
        val_temp = 1.0*sin(L*x[0])*sin(L*x[1])
        values[0] = base_val+amp*val_temp
        values[1] = (values[0]-0.1)*(values[0]-0.9)*(values[0]-0.5)-lmbda*amp*val_temp*(L**2)
    def value_shape(self):
        return (2,)

where the initial value for c at the boundary is 0.9, the dolfin solver has difficulty converging for

bc_c = DirichletBC(ME.sub(0), Constant(0.9), boundary)  # For 'c'
bc_mu = DirichletBC(ME.sub(1), Constant(0), boundary)  # For 'mu'

However, the dolfin solver accepts

bc_c = DirichletBC(ME.sub(0), Constant(0), boundary)  # For 'c'
bc_mu = DirichletBC(ME.sub(1), Constant(0), boundary)  # For 'mu'

and behaves as if the boundary value of c is actually fixed at 0.9 instead of 0. I checked this for various initial conditions by printing out the boundary values at different times of the simulation and the behavior is always the same. Is this some weird implementation issue for legacy fenics?