add 2d_polarization

"""

Material: BTO.

Order Parameter: Spontaneous Polarization.

Model: Phase Field with external electric field and strain.

Numerics: Fourier spectral method. 2D.

You can store all data as .h5 (Write_to_HDF5)

"""

from material_data import Mat_Data_BTO_D, initial_polarization

from fourier_frequencies import fourier_space

from solve_piezoelectricity import Solve_Piezoelectricity_2D

from green_operator import Green_Operator_Piezoelectricity_2D

from calculate_derivatives import Calculate_Derivatives_2D, Calculate_Derivatives_Land_Potential_2D

from make_plot import Plot_Polarization, Plot_Potential

from store_hdf5 import Write_to_HDF5

# ----------------------------------------------------------------------------

import numpy as np

import timeit

import itertools

#-----------------------------------------------------------------------------

def Evolve_Sponteneous_Polarization_2D(folder, Nx, Ny, nsteps, nt, dt, domain_type, E_app_y, E_app_x, eps_app_x, eps_app_y):

    # ----------------------------------------------------------------------------

    # import material data

    G, l_wall, P0, mob, e00, k_sep, k_grad, l_coef, e_piezo, C0, C_cub, K0 = Mat_Data_BTO_D()  

    dx, dy = (l_wall/5, l_wall/5)

    # ----------------------------------------------------------------------------

    # initial polarization

    P = initial_polarization(Nx, Ny, domain_type) * P0

    # ----------------------------------------------------------------------------

    # frequencies

    freq = fourier_space(Nx, dx, Ny, dy)

    # ----------------------------------------------------------------------------

    # save inital polarization    

    Plot_Polarization(folder + "/P_", 0, Nx, Ny, P[0], P[1], P0)

    #Write_to_HDF5(folder, 0, P)

    # ----------------------------------------------------------------------------

    # plot the potential

    Plot_Potential(folder, P0, Nx, l_coef[0], l_coef[1], l_coef[2], l_coef[3], l_coef[4])

    # ----------------------------------------------------------------------------

    # assign tensor for each grid

    C0 = np.einsum('ijkl, xy -> ijklxy', C0, np.ones((Nx,Ny)))

    K0 = np.einsum('ij, xy -> ijxy', K0, np.ones((Nx,Ny)))

    # ----------------------------------------------------------------------------

    eps_ext  = np.zeros((2,2,Nx,Ny))            # applied elas. field

    E_ext    = np.zeros((2,Nx,Ny))             # applied elec. field

    # ----------------------------------------------------------------------------

    eps_ext[0,0] = eps_app_x; eps_ext[1,1] = eps_app_y

    E_ext[0] = E_app_x; E_ext[1] = E_app_x

    # ----------------------------------------------------------------------------

    # random electric field, zero for perfect crystal

    E_random = np.zeros_like(E_ext)

    # ----------------------------------------------------------------------------

    # derivatives

    eP_val, deP_dPx_val, deP_dPy_val, eps0_val, deps0_P_dPx_val, deps0_P_dPy_val = Calculate_Derivatives_2D()

    Psi_val, dPsi_dPx_val, dPsi_dPy_val = Calculate_Derivatives_Land_Potential_2D()

    # ----------------------------------------------------------------------------

    # time evolution

    for step in range(nsteps):

        # ------------------------------------------------------------------------

        print('\n---------------------- Time step:\t' + str(step)+'\t----------------------')

        # -------------------------------------------------------------------------

        # calculate derivatives, convert symbolic derivatives into numpy

        eP = np.array(eP_val(P[0], P[1], P0, e_piezo[0], e_piezo[1], e_piezo[2]))

        deP_dPx = np.array(deP_dPx_val(P[0], P[1], P0, e_piezo[0], e_piezo[1], e_piezo[2]))

        deP_dPy = np.array(deP_dPy_val(P[0], P[1], P0, e_piezo[0], e_piezo[1], e_piezo[2]))

        # -------------------------------------------------------------------------

        eps0 = np.array(eps0_val(P[0], P[1], P0, e00))

        deps0_P_dPx = np.array(deps0_P_dPx_val(P[0], P[1], P0, e00))

        deps0_P_dPy = np.array(deps0_P_dPy_val(P[0], P[1], P0, e00))

        # -------------------------------------------------------------------------

        e0 = np.zeros((2,2,2,Nx,Ny))    # homogeneous piez. tensor, volume average of e(P)

        for i,j,k in itertools.product(range(2), repeat=3):

            e0[i,j,k] = np.mean(eP[i,j,k])

        # -------------------------------------------------------------------------

        # greens functions

        G_elas, G_piezo, G_dielec = Green_Operator_Piezoelectricity_2D(C0, e0, K0, freq)

        # ------------------------------------------------------------------------

        # Solve piezoelectricity

        sigma, D, eps_elas, E = Solve_Piezoelectricity_2D(C0, K0, e0, C0, K0, eP, G_elas, G_piezo, G_dielec, 

                                                          eps0, eps_ext, E_ext, E_random, Nx, Ny, P)

        # ------------------------------------------------------------------------

        del G_elas, G_piezo, G_dielec

        # Print energies

        # ------------------------------------------------------------------------

        H_elas = 0.5*np.einsum('ijkl...,ij...,kl...',C0, eps_elas, eps_elas)

        # ------------------------------------------------------------------------

        print('\nElastic Energy:\t\t\t\t\t\t'+ str(np.mean(H_elas)) + '\t[N/m2]')

        # ------------------------------------------------------------------------

        H_piezo = -1*np.einsum('ijk...,ij...,k...', eP, eps_elas, E)

        # ------------------------------------------------------------------------

        print('Piezoelectric Energy:\t\t\t\t\t'+ str(np.mean(H_piezo)) + '\t[N/m2]')

        # ------------------------------------------------------------------------

        H_elec = - 0.5*np.einsum('ij...,i...,j...', K0, E, E) - np.einsum('i...,i...', P, E)

        # ------------------------------------------------------------------------

        print('Electric Energy:\t\t\t\t\t'+ str(np.mean(H_elec)) + '\t[N/m2]')

        # ------------------------------------------------------------------------

        # Derivatives 

        deps_elas_dPx = -1*deps0_P_dPx; deps_elas_dPy = -1*deps0_P_dPy

        # ------------------------------------------------------------------------

        dH_elas_dPx = 0.5*(np.einsum('ijkl...,ij...,kl...', C0, deps_elas_dPx, eps_elas) + 

                           np.einsum('ijkl...,ij...,kl...', C0, eps_elas, deps_elas_dPx) )

        dH_elas_dPy = 0.5*(np.einsum('ijkl...,ij...,kl...', C0, deps_elas_dPy, eps_elas) + 

                           np.einsum('ijkl...,ij...,kl...', C0, eps_elas, deps_elas_dPy) )

        # ------------------------------------------------------------------------

        dH_piezo_dPx = -1*(np.einsum('ijk...,ij...,k...', deP_dPx, eps_elas, E) + 

                           np.einsum('ijk...,ij...,k...', eP, deps_elas_dPx, E) )

        dH_piezo_dPy = -1*(np.einsum('ijk...,ij...,k...', deP_dPy, eps_elas, E) + 

                           np.einsum('ijk...,ij...,k...', eP, deps_elas_dPy, E) )

        # ------------------------------------------------------------------------

        dH_elec_dPx = -1*E[0]

        dH_elec_dPy = -1*E[1]

        # ------------------------------------------------------------------------

        Psi = np.array(Psi_val(P[0], P[1], P0, l_coef[0], l_coef[1], l_coef[2], l_coef[3], l_coef[4]))

        # ------------------------------------------------------------------------

        H_sep = k_sep * G/l_wall * Psi

        # ------------------------------------------------------------------------

        print('Separation Energy:\t\t\t\t\t'+ str(np.mean(H_sep)) + '\t[N/m2]')

        # ------------------------------------------------------------------------

        # ------------------------------------------------------------------------

        dPsi_dPx = np.array(dPsi_dPx_val(P[0], P[1], P0, l_coef[1], l_coef[2], l_coef[3], l_coef[4]))

        dPsi_dPy = np.array(dPsi_dPy_val(P[0], P[1], P0, l_coef[1], l_coef[2], l_coef[3], l_coef[4]))

        # ------------------------------------------------------------------------

        dH_dPx = k_sep * G/l_wall * dPsi_dPx + dH_elas_dPx + dH_piezo_dPx + dH_elec_dPx 

        dH_dPy = k_sep * G/l_wall * dPsi_dPy + dH_elas_dPy + dH_piezo_dPy + dH_elec_dPy 

        # ------------------------------------------------------------------------

        dH_dPx_k = np.fft.fft2(dH_dPx)

        dH_dPy_k = np.fft.fft2(dH_dPy)

        # ------------------------------------------------------------------------

        Px_k = np.fft.fft2(P[0]); Py_k = np.fft.fft2(P[1])

        # ------------------------------------------------------------------------

        dPx_dx = np.fft.ifft2(1j*freq[0]*Px_k).real

        dPx_dy = np.fft.ifft2(1j*freq[1]*Px_k).real

        dPy_dx = np.fft.ifft2(1j*freq[0]*Py_k).real

        dPy_dy = np.fft.ifft2(1j*freq[1]*Py_k).real

        # ------------------------------------------------------------------------

        H_grad = 0.5*k_grad*G*l_wall/P0**2*(dPx_dx**2 + dPx_dy**2 + dPy_dx**2 + dPy_dy**2)

        print('Gradient Energy:\t\t\t\t'+str(np.mean(H_grad)))

        # ------------------------------------------------------------------------

        del dH_elas_dPx, dH_elas_dPy, dH_piezo_dPx, dH_piezo_dPy, dH_elec_dPx, dH_elec_dPy

        del dH_dPx, dH_dPy, dPx_dx, dPx_dy, dPy_dx, dPy_dy, dPsi_dPx, dPsi_dPy

        # ------------------------------------------------------------------------

        Px_k = (Px_k - dt * mob * dH_dPx_k) / (1 + dt * mob * k_grad *G*l_wall/P0**2 * (freq[0]**2 + freq[1]**2))

        Py_k = (Py_k - dt * mob * dH_dPy_k) / (1 + dt * mob * k_grad *G*l_wall/P0**2 * (freq[0]**2 + freq[1]**2))

        # ------------------------------------------------------------------------

        Px = np.fft.ifft2(Px_k).real

        Py = np.fft.ifft2(Py_k).real

        # ------------------------------------------------------------------------

        P = np.array([Px, Py])

        # ------------------------------------------------------------------------

        if ((step + 1) % nt == 0 ):

            # --------------------------------------------------------------------

            Plot_Polarization(folder + "/P_", step + 1, Nx, Ny, P[0], P[1], P0)

            """

            Use the following if you want to save any of them:

                Stress=sigma

                Elas_strain=eps_elas

                Elec_disp=D

                Elec_field=E

                Elas_en=H_elas

                Grad_en=H_grad

                Sep_en=H_sep

                Elec_en=H_elec

                Piezo_en=H_piezo 

            """

            #Write_to_HDF5(folder, step + 1, P)#, Grad_en=H_grad)

            # --------------------------------------------------------------------

        del sigma, D, eps_elas, E 

        del H_elas, H_piezo, H_elec, H_grad, Psi, H_sep

    # ----------------------------------------------------------------------------

def main():

    start_tm = timeit.default_timer()

    # ----------------------------------------------------------------------------

    # grid parameters

    Nx, Ny = (360, 360)

    # time step parameters

    nsteps, nt, dt = (50000, 100, 1e-12)   # dt is in sec

    # applied fields in each direction

    E_app_x, E_app_y = (0, 0)  # must be in V/m

    eps_app_x, eps_app_y = (0, 0)

    # ----------------------------------------------------------------------------

    # folder to save results

    import os

    folder = os.getcwd() 

    # ----------------------------------------------------------------------------

    domain_type = "random" # domain_type = "ramdom" OR "90" OR "180"

    # ----------------------------------------------------------------------------

    # Evolve

    Evolve_Sponteneous_Polarization_2D(folder, Nx, Ny, nsteps, nt, dt, domain_type, E_app_y, E_app_x, eps_app_x, eps_app_y)

    # ----------------------------------------------------------------------------

    stop_tm = timeit.default_timer()

    print('\nExecution time: ' + str( format((stop_tm - start_tm)/60,'.5f')) + ' min\n' ) 

if __name__ == "__main__":

    main()

 No newline at end of file