# Initial and Dirichlet Condition on Chemical Potential for the CH equation

**URL:** <https://fenicsproject.discourse.group/t/initial-and-dirichlet-condition-on-chemical-potential-for-the-ch-equation/15947>\
**Category:** General\
**Tags:** legacy-dolfin\
**Created:** [October 1, 2024, 8:45pm UTC](https://fenicsproject.discourse.group/t/initial-and-dirichlet-condition-on-chemical-potential-for-the-ch-equation/15947 "2024-10-01T20:45:34Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![andrewli1235](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/andrewli1235/32/7033_2.png) [@andrewli1235](https://fenicsproject.discourse.group/u/andrewli1235)\
**Post date:** [October 1, 2024, 8:45pm UTC](https://fenicsproject.discourse.group/t/initial-and-dirichlet-condition-on-chemical-potential-for-the-ch-equation/15947/1 "2024-10-01T20:45:34Z")

</div>

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](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.

```auto
 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

```auto
 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)`?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 2, 2024, 5:39am UTC](https://fenicsproject.discourse.group/t/initial-and-dirichlet-condition-on-chemical-potential-for-the-ch-equation/15947/2 "2024-10-02T05:39:13Z")

</div>

> [@andrewli1235](#):
>
> 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](https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/cahn-hilliard/python/documentation.html)

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](https://fenicsproject.org/olddocs/dolfin/latest/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html)

> [@andrewli1235](#):
>
> ```auto
> 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,)
> 
> ```

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

> [@andrewli1235](#):
>
> `(c-0.9)*(c-0.1)*(c-0.5)-lmbda*div(grad(c)`?

I would probably look at: [A conservative numerical method for the Cahn–Hilliard equation with Dirichlet boundary conditions in complex domains - ScienceDirect](https://www.sciencedirect.com/science/article/pii/S0898122112006463?via%3Dihub) and references therein

---

<div class="post-metadata">

**Author:** ![andrewli1235](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/andrewli1235/32/7033_2.png) [@andrewli1235](https://fenicsproject.discourse.group/u/andrewli1235)\
**Post date:** [October 2, 2024, 4:37pm UTC](https://fenicsproject.discourse.group/t/initial-and-dirichlet-condition-on-chemical-potential-for-the-ch-equation/15947/3 "2024-10-02T16:37:17Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 2, 2024, 5:52pm UTC](https://fenicsproject.discourse.group/t/initial-and-dirichlet-condition-on-chemical-potential-for-the-ch-equation/15947/4 "2024-10-02T17:52:36Z")

</div>

> [@andrewli1235](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![andrewli1235](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/andrewli1235/32/7033_2.png) [@andrewli1235](https://fenicsproject.discourse.group/u/andrewli1235)\
**Post date:** [October 7, 2024, 5:43pm UTC](https://fenicsproject.discourse.group/t/initial-and-dirichlet-condition-on-chemical-potential-for-the-ch-equation/15947/5 "2024-10-07T17:43:27Z")

</div>

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

```auto
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

```auto
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

```auto
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?