Thursday, November 10, 2011

SNOW: Robust simulation of strongly anisotropic Cahn-Hilliard systems


Purpose
  • The package "Snow" solves the strongly anisotropic Cahn-Hilliard system (regularized by the nonlinear Willmore functional) on a 2-D square domain. It uses a Legendre-Galerkin method for systems of coupled elliptic equations in space and a first-order energy stabilized finite difference scheme in time.

Specifications
  • Name: Snow.
  • Author: Feng Chen.
  • Finishing date: 11/10/2011.
  • Languages: Fortran 90, MATLAB.
  • Required libraries: BLAS, LAPACK.

Simple Example
  • Equation: \begin{equation} \left\{ \begin{aligned} & \phi_{t} =\frac{1}{\epsilon}\Delta\mu, && \quad \text{in } \Omega, \\ & \mu =\frac{\delta\mathcal E}{\delta\phi}, && \quad \text{in } \Omega, \\ & \omega =\frac{1}{\epsilon}f'(\phi)-\epsilon\Delta\phi, && \quad \text{in } \Omega, \\ & \frac{\partial \phi}{\partial \boldsymbol n} = \frac{\partial \mu}{\partial \boldsymbol n} = \frac{\partial \omega}{\partial \boldsymbol n} = 0, && \quad \text{on } \partial\Omega , \end{aligned} \right. \end{equation} where \begin{equation} \begin{aligned} & \mathcal E(\phi)= \int_{\Omega}(\mathcal F+ \frac{\beta}{2} \mathcal G)\text{d} \Omega,\\ & \mathcal F(\phi)= \frac{\gamma(\boldsymbol n)}{\epsilon}(f(\phi)+\frac{\epsilon^{2}}{2}|\nabla\phi|^{2}) ,\\ &\mathcal G (\phi) = \omega^2,\\ &f(\phi) = \frac{1}{4} (\phi^2-1)^2,\\ &\omega=\frac{1}{\epsilon}f'(\phi)-\epsilon\Delta\phi, \\ & \gamma ( \boldsymbol n) =1+\alpha\cos(4\theta)= 1 + \alpha (4 \sum_{i=1}^d n_i^4 -3). \end{aligned} \end{equation}
  • Initial conditions: \begin{equation} \begin{aligned} &\phi_0 = \tanh (\frac{\sqrt{x^2+y^2}-r}{\epsilon}), \\ &\omega_0 = f'(\phi_0). \end{aligned} \end{equation}
  • Parameters: \begin{equation} \Omega = (-1,1)^2 , \, T=0.15, \, r=0.5, \, \epsilon=0.02, \, \beta=5\times10^{-4}, \, \alpha=0.3. \end{equation}

Quick Start
  • Compiling and running:
    cd ./Lion
    make library
    ifort tst_Snow.f90 -llibrary -llapack -lblas
    ./a.out
    matlab plot_Snow.m
    
  • Graphics: The left panel shows the evolution of the anisotropic dynamics in $\Omega$, and the right one shows the decreasing energy against time steps.

References
Code Highlights
    
    !! Assign initial conditions
    do i = 0, nx
     do j = 0, ny
      ! quadrature points
      x = (zx(i)*(x2-x1) + (x2+x1))/2d0
      y = (zy(j)*(y2-y1) + (y2+y1))/2d0
      select case (SnowAnisoCase) 
       case (1)  ! one circle
        phi(i,j) = tanh((sqrt(x*x+y*y)-r)/SnowEps)
       case (2)  ! two circles
        phi(i,j) = tanh((sqrt((x+1d0/2d0)**2+                   &
                              (y-1d0/2d0)**2)-1d0/4d0)/SnowEps) &
                 + tanh((sqrt((x-1d0/5d0)**2+                   &
                              (y+1d0/5d0)**2)-1d0/2d0)/SnowEps)
        phi(i,j) = phi(i,j) - 1d0
      end select 
     end do
    end do     
    omega =- (phi*phi*phi-phi)/SnowEps      
    
    

No comments:

Post a Comment