Merge branch 'develop' into 'main'

"""

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

"""

Symbolic differentiation of the spontaneous strain and piezoelectric tensor wrt P 

"""

import sympy as sym

import numpy as np

import itertools

def Calculate_Derivatives_2D():

    ndim = 2

    Px = sym.Symbol('Px')

    Py = sym.Symbol('Py')

    P0 = sym.Symbol('P0')

    e33 = sym.Symbol('e33')

    e31 = sym.Symbol('e31')

    e15 = sym.Symbol('e15')

    e00 = sym.Symbol('e00')

    norm = sym.sqrt(Px**2 + Py**2)

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

    I = np.eye(ndim)

    p = P/norm

    # Piezoelectric tensor

    eP = np.zeros((ndim,ndim,ndim),object)

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

        eP[k,i,j] = (norm/P0)**3*(e33*p[i]*p[j]*p[k] + e31*(I[i,j]-p[i]*p[j])*p[k]

                                      + e15*0.5*( (I[k,i]-p[k]*p[i])*p[j] + (I[k,j]-p[k]*p[j])*p[i] ))

    # derivative of eP wrt Px

    deP_dPx = np.zeros_like(eP, object)

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

        deP_dPx[k,i,j] = sym.diff(eP[k,i,j], Px)

    # derivative of eP wrt Py

    deP_dPy = np.zeros_like(eP, object)

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

        deP_dPy[k,i,j] = sym.diff(eP[k,i,j], Py)

    # calculate spon. strain

    eps0_P = np.zeros((ndim,ndim), object)

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

        eps0_P[i,j] = 3/2*e00*(norm/P0)**2*(p[i]*p[j] - 1/3*I[i,j])

    # get its deivative wrt Px

    deps0_P_dPx = np.zeros_like(eps0_P, object)

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

        deps0_P_dPx[i,j] = sym.diff(eps0_P[i,j], Px)

    # get its deivative wrt Py

    deps0_P_dPy = np.zeros_like(eps0_P, object)

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

        deps0_P_dPy[i,j] = sym.diff(eps0_P[i,j], Py)

    eP_val = sym.lambdify([Px, Py, P0, e33, e31, e15], eP, 'numpy')

    deP_dPx_val = sym.lambdify([Px, Py, P0, e33, e31, e15], deP_dPx, 'numpy')

    deP_dPy_val = sym.lambdify([Px, Py, P0, e33, e31, e15], deP_dPy, 'numpy')

    eps0_val = sym.lambdify([Px, Py, P0, e00], eps0_P, 'numpy')

    deps0_P_dPx_val = sym.lambdify([Px, Py, P0, e00], deps0_P_dPx, 'numpy')

    deps0_P_dPy_val = sym.lambdify([Px, Py, P0, e00], deps0_P_dPy, 'numpy')

    del eP, deP_dPx, deP_dPy, eps0_P, deps0_P_dPy, deps0_P_dPx

    return eP_val, deP_dPx_val, deP_dPy_val, eps0_val, deps0_P_dPx_val, deps0_P_dPy_val

def Calculate_Derivatives_Land_Potential_2D():

    Px = sym.Symbol('Px')

    Py = sym.Symbol('Py')

    P0 = sym.Symbol('P0')

    a0 = sym.Symbol('a0')

    a1 = sym.Symbol('a1')

    a2 = sym.Symbol('a2')

    a3 = sym.Symbol('a3')

    a4 = sym.Symbol('a4')

    Psi = a0 + a1/P0**2 * (Px**2 + Py**2) + a2/P0**4 * (Px**4 + Py**4) \

             + a3/P0**4 * Px**2 * Py**2 + a4/P0**6 * (Px**6 + Py**6)

    dPsi_dPx = sym.diff(Psi, Px)

    dPsi_dPy = sym.diff(Psi, Py)

    Psi_val = sym.lambdify([Px, Py, P0, a0, a1, a2, a3, a4], Psi, 'numpy')

    dPsi_dPx_val = sym.lambdify([Px, Py, P0, a1, a2, a3, a4], dPsi_dPx, 'numpy')

    dPsi_dPy_val = sym.lambdify([Px, Py, P0, a1, a2, a3, a4], dPsi_dPy, 'numpy')

    del Psi, dPsi_dPx, dPsi_dPy

    return Psi_val, dPsi_dPx_val, dPsi_dPy_val

"""

Calculates frequencies in Fourier space

Returns the vector of frequencies

"""

import numpy as np

def fourier_space(Nx, dx, Ny, dy, Nz=None, dz=None):

    if Nz == None:

        kx = 2.0*np.pi*np.fft.fftfreq(Nx, dx)

        ky = 2.0*np.pi*np.fft.fftfreq(Ny, dy)

        kx_grid, ky_grid = np.meshgrid(kx,ky, indexing='ij')

        freq = np.array([kx_grid, ky_grid])

        del kx, ky, kx_grid, ky_grid

        return freq

    else:

        kx = 2.0*np.pi*np.fft.fftfreq(Nx, dx)

        ky = 2.0*np.pi*np.fft.fftfreq(Ny, dy)

        kz = 2.0*np.pi*np.fft.fftfreq(Nz, dz)

        kx_grid, ky_grid, kz_grid = np.meshgrid(kx, ky, kz, indexing='ij')

        freq = np.array([kx_grid, ky_grid, kz_grid])

        del kx, ky, kz, kx_grid, ky_grid, kz_grid

        return freq

 No newline at end of file

"""

Calculates the Green's functions in Fourier space

"""

import numpy as np

import itertools

def Green_Operator_Piezoelectricity_2D(C0, e0, K0, freq):

    A_c = np.einsum('pijqx..., px..., qx... -> ijx...', C0, freq, freq)

    A_e = np.einsum('piqx...,  px..., qx... -> ix... ', e0, freq, freq)

    A_d = np.einsum('pqx...,   px..., qx... -> x...  ', K0, freq, freq)

    A_d[0,0] = 1.0

    invA_d = 1/A_d   

    A = A_c + invA_d * np.einsum('ix..., jx... -> ijx... ', A_e, A_e)

    adjG = np.empty_like(A)

    adjG[0, 0] = A[1, 1]

    adjG[1, 1] = A[0, 0]

    adjG[0, 1] = -A[0, 1]

    adjG[1, 0] = -A[1, 0]

    detG = A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0]

    detG[0,0] = 1.0

    invA = adjG/detG

    G_elas = np.zeros_like(C0) # Green op. for elasticity

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

        G_elas[i,j,k,l] = 0.25*(invA[i,l]*freq[k]*freq[j]+invA[j,l]*freq[k]*freq[i]+invA[i,k]*freq[l]*freq[j]+invA[j,k]*freq[i]*freq[l])

    G_piezo = np.zeros_like(e0)     # Green op. for piezoelectricity

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

        G_piezo[k,i,j] =( 1/2 * invA_d *freq[k]*(invA[i,0]*freq[j] + invA[j,0]*freq[i])*A_e[0] + 

                                  1/2 * invA_d * freq[k]*(invA[i,1]*freq[j] + invA[j,1]*freq[i])*A_e[1] )

    G_dielec = invA_d**2*(np.einsum('i...,ij...,j...',A_e,invA,A_e) - A_d)*np.einsum('ix..., jx... -> ijx... ', freq, freq)

    del A_c, A_e, A_d, invA_d, A, adjG, detG, invA

    return G_elas, G_piezo, G_dielec

"""

This script is used for postprocessing:

    - plots polarization magnitude with its directions as an image.

    - plots 6th-order Landau potential

"""

import matplotlib.pyplot as plt

import numpy as np

from scipy import ndimage

def Plot_Polarization(variable, tt, Nx, Ny, P1, P2, max_value):

    x,y = np.meshgrid(np.linspace(0,Nx,Nx),np.linspace(0,Ny,Ny))   

    u = ndimage.rotate(P1, 90)

    v = ndimage.rotate(P2, 90)

    r2 = np.power(np.add(np.power(u,2),np.power(v,2)),0.5)

    widths = np.linspace(0, 1, x.size)

    fig, ax = plt.subplots(figsize=(10, 10))

    P_mag = np.sqrt(P1**2 + P2**2)

    P_mag = ndimage.rotate(P_mag, 90)

    #"""

    im = ax.imshow(P_mag, cmap='jet', vmin=0, vmax=max_value)

    ax.set_xlabel('x')

    ax.set_ylabel('y')

    #"""

    # colorbar

    cax = fig.add_axes([ax.get_position().x1+0.01,ax.get_position().y0,0.02,ax.get_position().height])

    cb = plt.colorbar(im,cax=cax)

    tick_font_size = 20

    cb.ax.tick_params(labelsize=tick_font_size)

    u = u/r2; v= v/r2

    n = 14

    ax.quiver(x[::n,::n],y[::n,::n],u[::n,::n],v[::n,::n],linewidths=widths,scale=3.8, 

              scale_units='inches', color='white')

    ax.axis('off')

    fig.savefig(variable + str(tt), dpi=fig.dpi, bbox_inches='tight')#, pad_inches=0.5)

    plt.close(fig)

    del u, v, P_mag

def Plot_Potential(folder, P0, Nx, a0, a1, a2, a3, a4):

    from matplotlib import cm

    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

    # Make data.

    Px = np.linspace(-P0, P0, Nx)

    Py = np.linspace(-P0, P0, Nx)

    Px, Py = np.meshgrid(Px, Py)

    Psi = a0 + a1/P0**2 * (Px**2 + Py**2) + a2/P0**4 * (Px**4 + Py**4) \

             + a3/P0**4 * (Px**2 * Py**2) + a4/P0**6 * (Px**6 + Py**6)

    # Plot the surface.

    ax.plot_surface(Px, Py, Psi, cmap=cm.coolwarm,

                           linewidth=0, antialiased=False)

    ax.zaxis.set_major_formatter('{x:.02f}')

    fig.savefig(folder + "/Landau_potential", dpi=fig.dpi, bbox_inches='tight')

    plt.close(fig)

 No newline at end of file