2
$\begingroup$

So, I have written both a C and python code to solve the 2D Cahn-Hilliard equation: \begin{equation} \frac{\partial c}{\partial t} = \nabla^2\left(c^3 - c - \kappa\nabla^2c\right) \end{equation} Taking its Fourier transform yields the following ODE: \begin{equation} \frac{\partial \hat c}{\partial t} = \underbrace{-k^2\{c^3-c\}_\mathbf{k}}_{F\left(\hat c\right)} \overbrace{- \kappa k^4\hat c}^{G\left(\hat c\right)} \end{equation} Where G is linear and "fast-varying" (high eigenvalues) and F "slow-varying" (low eigenvalues) and $\mathbf k = \sqrt{k_x^2 + k_y^2}$ the Fourier "wavelength". To overcome the highly-constraining 4th-order term in terms of timestep size, I decided to use the following second order implicit-explicit "IMEX" scheme where $\hat f \equiv \{c^3 - c\}_\mathbf k$: \begin{equation} \left(3+2\Delta t\kappa k^4\right)\hat c^{n+1} = 4\hat c^n - \hat c^{n-1} - 2\Delta t k^2\left(2\hat f^n - \hat f^{n-1}\right) \end{equation} Where the iterations are started using the first order IMEX-BDF scheme. Now, this works perfectly well and exhibits second-order error convergence as expected. I then wanted to try and use the fourth-order IMEX-BDF scheme to obtain even better accuracy at basically no extra CPU cost (it requires more memory to store the previous steps, but memory is not the issue in this case): \begin{equation} \left(25+12\Delta t\kappa k^4\right)\hat c^{n+1} = 48\hat c^n - 36\hat c^{n-1} + 16\hat c^{n-2} - 3\hat c^{n-3} - 12\Delta t k^2\left(4\hat f^n - 6\hat f^{n-1} + 4\hat f^{n-2} - \hat f^{n-3}\right) \end{equation} However, surprisingly to me, this scheme yields a slightly worst second order error convergence behavior than the second order scheme (slightly worst in that the error constant looks to be bigger for the IMEX-BDF4 scheme). I have observed the same behavior for the IMEX-BDF3 scheme (second order error convergence instead of third; slightly worst error constant than IMEX-BDF2 - but a bit better than IMEX-BDF4).

To validate the convergence order, I took as reference solution a fourth-order classical Runge-Kutta scheme with very low $\Delta t$.

Thus my question is the following: is there any reason why the 4th- and 3rd-order IMEX-BDF schemes would be somehow limited to second order in the particular use case of the Cahn-Hilliard equation ? Or is it me who implemented something the wrong way ? Since I'm using Fourier spectral methods for physical differentiation, I'm assuming the physical discretisation is not a limiting factor in terms of simulation precision since I'm using a fairly high number of mesh points - 256x256.

My python code w/ IMEX2, IMEX4 & RK4: pastebin

References:

Willem Hundsdorfera, Steven J. Ruuth, IMEX extensions of linear multistep methods with general monotonicity and boundedness properties

Jingzhi Zhu, Long-Qing Chen, Jie Shen, and Veena Tikare, Coarsening kinetics from a variable-mobility Cahn-Hilliard equation: Application of a semi-implicit Fourier spectral method

$\endgroup$
2
$\begingroup$

Based on the references you provided, I assume your initial conditions are random numbers (see second reference, Sect. IV)? If so, it's probably no surprise that you do not see higher-order convergence in time, at least not from the beginning. Your initial conditions are not smooth and it may take some time to get to a smooth solution/state. There, however, higher-order convergence should be seen. I recommend starting with a smooth and periodic initial condition to check convergence order.

In case your initial conditions are different (i.e. smooth and periodic), please specify them and share the code, ideally as a minimal working example.

$\endgroup$
6
  • $\begingroup$ Indeed, my initial conditions are random. And I've actually managed to obtain fourth-order convergence behavior using classical Runge-Kutta (explicit) scheme as well as using (a new one; had not tried that out when I wrote my original question) ETDRK4 (exponential time differencing RK4) scheme - which has the same "great stability" properties as the IMEX scheme I described. All this leads me to believe I've simply made some mistake somewhere when implementing the IMEX schemes and that I somehow can't manage to get it working properly. $\endgroup$ – Gilles Poncelet Mar 24 at 22:45
  • $\begingroup$ I've updated my original question with my code for solving the C-H equation in case some people are interested. I've made a few edits about my comments about convergence: I do have second-order convergence for both 2nd- and 4th-order IMEX schemes. Maybe you are right and it is an issue about waiting for a longer "steadier" state for IMEX schemes, but it still surprises me compared to RK4 & ETDRK4. $\endgroup$ – Gilles Poncelet Mar 24 at 23:05
  • $\begingroup$ Is there a k missing in line 77 of your code? Your formula says k**2, but the code has only k. $\endgroup$ – Robert Speck Mar 25 at 6:44
  • $\begingroup$ No there (normally) is not: that’s a lack of clarity on my part in my code. I computed k as k^2 in my code (you can see that line 120) so to avoid repetitive operations. (So when I write *k it’s actually *k^2 and *k^2 => *k^4). $\endgroup$ – Gilles Poncelet Mar 25 at 6:47
  • $\begingroup$ OK, got it. Sure your signs in line 69 are correct? I'd assume a +6*dt*k in front of the brackets. $\endgroup$ – Robert Speck Mar 25 at 7:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.