Cahn-Hilliard

Cahn-Hilliard.py

   @overrides
    def solve(self): #-> np.ndarray:
        """
        The sequential version of Cahn-Hilliard solver.

        Update a conserved order parameter, in our case, the concetration field $c(r, t)$.

        Calculate free energy derivative at domain nodes.

        The eight neighbour nodes of the central node (C) are north (N), south (S),
        east(E), west (W), nortwest (NW), northeast (NE),
        southwest (SW), southeast (SE).

            NW---N---EN
            |   |   |
            W---C---E
            |   |   |
            SW---S---SE

        """
        c0 = 0.5  # average composition of B atom [atomic fraction]

        from microtex.modeling.quantities import R

        dx, dy, dt, La, T, ac, Da, Db, nx, ny = (
            self.dx,
            self.dy,
            self.dt,
            self.La,
            self.T,
            self.ac,
            self.Da,
            self.Db,
            self.nx,
            self.ny,
        )

        c = self.state[-1]  # The last field alias.

        for j in range(ny):
            for i in range(nx):

                ip = i + 1
                im = i - 1
                jp = j + 1
                jm = j - 1
                ipp = i + 2
                imm = i - 2
                jpp = j + 2
                jmm = j - 2

                if ip > nx - 1:  # periodic boundary condition
                    ip = ip - nx
                if im < 0:
                    im = im + nx
                if jp > ny - 1:
                    jp = jp - ny
                if jm < 0:
                    jm = jm + ny
                if ipp > nx - 1:
                    ipp = ipp - nx
                if imm < 0:
                    imm = imm + nx
                if jpp > ny - 1:
                    jpp = jpp - ny
                if jmm < 0:
                    jmm = jmm + ny

                cc = c[i, j]  # at (i,j) "centeral point"
                ce = c[ip, j]  # at (i+1.j) "eastern point"
                cw = c[im, j]  # at (i-1,j) "western point"
                cs = c[i, jm]  # at (i,j-1) "southern point"
                cn = c[i, jp]  # at (i,j+1) "northern point"
                cse = c[ip, jm]  # at (i+1, j-1)
                cne = c[ip, jp]
                csw = c[im, jm]
                cnw = c[im, jp]
                cee = c[ipp, j]  # at (i+2, j+1)
                cww = c[imm, j]
                css = c[i, jmm]
                cnn = c[i, jpp]

                mu_chem_c = R * T * (np.log(cc) - np.log(1.0 - cc)) + La * (
                    1.0 - 2.0 * cc
                )  # chemical term of the diffusion potential
                mu_chem_w = R * T * (np.log(cw) - np.log(1.0 - cw)) + La * (
                    1.0 - 2.0 * cw
                )
                mu_chem_e = R * T * (np.log(ce) - np.log(1.0 - ce)) + La * (
                    1.0 - 2.0 * ce
                )
                mu_chem_n = R * T * (np.log(cn) - np.log(1.0 - cn)) + La * (
                    1.0 - 2.0 * cn
                )
                mu_chem_s = R * T * (np.log(cs) - np.log(1.0 - cs)) + La * (
                    1.0 - 2.0 * cs
                )

                mu_grad_c = -ac * (
                    (ce - 2.0 * cc + cw) / dx / dx + (cn - 2.0 * cc + cs) / dy / dy
                )  # gradient term of the diffusion potential
                mu_grad_w = -ac * (
                    (cc - 2.0 * cw + cww) / dx / dx
                    + (cnw - 2.0 * cw + csw) / dy / dy
                )
                mu_grad_e = -ac * (
                    (cee - 2.0 * ce + cc) / dx / dx
                    + (cne - 2.0 * ce + cse) / dy / dy
                )
                mu_grad_n = -ac * (
                    (cne - 2.0 * cn + cnw) / dx / dx
                    + (cnn - 2.0 * cn + cc) / dy / dy
                )
                mu_grad_s = -ac * (
                    (cse - 2.0 * cs + csw) / dx / dx
                    + (cc - 2.0 * cs + css) / dy / dy
                )

                mu_c = mu_chem_c + mu_grad_c  # total diffusion potental
                mu_w = mu_chem_w + mu_grad_w
                mu_e = mu_chem_e + mu_grad_e
                mu_n = mu_chem_n + mu_grad_n
                mu_s = mu_chem_s + mu_grad_s

                nabla_mu = (mu_w - 2.0 * mu_c + mu_e) / dx / dx + (
                    mu_n - 2.0 * mu_c + mu_s
                ) / dy / dy
                dc2dx2 = ((ce - cw) * (mu_e - mu_w)) / (4.0 * dx * dx)
                dc2dy2 = ((cn - cs) * (mu_n - mu_s)) / (4.0 * dy * dy)

                DbDa = Db / Da
                mob = (Da / R / T) * (cc + DbDa * (1.0 - cc)) * cc * (1.0 - cc)
                dmdc = (Da / R / T) * (
                    (1.0 - DbDa) * cc * (1.0 - cc)
                    + (cc + DbDa * (1.0 - cc)) * (1.0 - 2.0 * cc)
                )

                # right-hand side of Cahn-Hilliard equation
                dcdt = mob * nabla_mu + dmdc * (
                    dc2dx2 + dc2dy2
                )

                # Update order parameter c.

                self.__empty[i, j] = c[i, j] + dcdt * dt

        self.__state.append(self.__empty.copy())