Different results when using divergence(gradient( . )) instead of laplace( . )

[ole-kranzosch]

Hello everyone,

I was trying to implement a more advanced version of the Cahn-Hilliard equation.
There i stumbled across a problem with the imposed zero Neumann boundary conditions, as they were not fullfilled.
The PDE looks like this: c_t = divergence( M(c) gradient(mu(c) - lambda * laplace(c)))
With c being a scalarfield, M and mu being functions of c, lambda being a constant, and c_t the time derivative of c.
After some testing I found the problem in the use of divergence of a gradient.
When using M(c) = 1 using divergence(gradient(mu(c) - lambda*laplace(c))) I get the following 1D result:

When replacing divergence(gradient(...)) with laplace(...) I get the expected (True) results:

The problem now is that I want to use a non-constant M(c) which needs the use of divergence(gradient).

Any help would be appreciated.

Best regards
Ole

My code:

    N = 128  # number of grid points
    gamma = 0.001  # interface width
    epsilon = 0.001  # interface width
    theta = 1.5
    T = 1  # end time
    delT = T/100
    # generate grid
    grid = pde.CartesianGrid([(0, 1)], N, periodic=False)

    # generate initial condition
    state = pde.ScalarField.from_expression(grid, "x")
    eq1 = pde.PDE({
        'c': f'laplace((1 / (2 * {theta})) * (ln(c) - ln(1-c)) + 1 -2 * c - {gamma} * laplace(c))'},
        bc="neumann")  # define the pde
    eq2 = pde.PDE({
        'c': f'divergence(gradient((1 / (2 * {theta})) * (ln(c) - ln(1-c)) + 1 -2 * c - {gamma} * laplace(c)))'},
        bc="neumann")  # define the pde

    storage = pde.MemoryStorage()
    result = eq1.solve(state, t_range=T, dt=1e-3, adaptive=True,
                      tracker=["progress", storage.tracker(delT)])  # solve the pde

[david-zwicker]

This is unfortunately a common problem with finite differences. Another problem you might find is that the total mass (the integral of c) is not conserved, although it obviously should be. A workaround could be to split the divergence operator, i.e., use

gradient(M'(c)) @ gradient(mu(c) - lambda * laplace(c))) +  M(c) * laplace(mu(c) - lambda * laplace(c)))

Note that this is not tested, but I hope it transports the idea.

In our research, we encounter similar problems and we have developed solutions, which are unfortunately not yet part of py-pde. If you encountered this problem for your own research, we could talk about a collaboration!

[ole-kranzosch]

Thank you for the reply. I also thought about splitting the equation and tested it today. The results are looking good so far, at least as expected.
Regarding the problem our aim is to verify/ test a FEM framework and for that i looked into the py-pde package. So avoiding the divergence gradient chain works for me.
Best regards
Ole Kranzosch

[david-zwicker]

That's great to hear! Good luck with your work :)