gradient_flow_omega-stab-Skyrmion.py

import numpy as np
from numba import jit
from numba.typed import Dict

@jit(nopython=True, cache=True, fastmath=True)
def func_solution_f_1st(r):
    return np.pi/(1.0+0.561769*r+14.1245*r**2)

@jit(nopython=True, cache=True, fastmath=True)
def func_solution_f_2nd(r):
    return 3.04269/(1.0+0.0691*r+14.8899*r**2)

@jit(nopython=True, cache=True, fastmath=True)
def func_solution_omega(r):
    return -0.921259/(1.0+0.0694971*r+7.69719*r**2+2.66569*r**3+67.413*r**4)

@jit(nopython=True, cache=True, fastmath=True)
def func_gradflow_f_flat(f,df,ddf,domega,r,mpi,fpi,beta):
    return (ddf+2.0*df/r-(np.sin(2.0*f)/r**2+mpi**2*np.sin(f))+2.0*beta/np.pi**2/fpi**2*np.sin(f)**2*domega/r**2)

@jit(nopython=True, cache=True, fastmath=True)
def func_gradflow_omega_flat(f,df,omega,domega,ddomega,r,momega,beta):
    return (ddomega+2.0*domega/r-momega**2*omega+beta/2.0/np.pi**2*np.sin(f)**2*df/r**2)

@jit(nopython=True, cache=True, fastmath=True)
def func_initconf_f(r):
    return np.pi/(1.0+0.561769*r+14.1245*r**2)

@jit(nopython=True, cache=True, fastmath=True)
def func_initconf_omega(r):
    return -0.921259/(1.0+67.413*r**4)

##########################

@jit(nopython=True, cache=True, fastmath=True)
def fwd_gradflow(f,omega,const_phys,const_calc):

    #####
    mpi=const_phys["mpi"]; momega=const_phys["momega"]; fpi=const_phys["fpi"]; beta=const_phys["beta"]
    N=int(const_calc["N"]); dr=const_calc["dr"]; reg=const_calc["reg"]; alpha_f=const_calc["alpha_f"]; alpha_omega=const_calc["alpha_omega"]
    #####

    df,ddf,domega,ddomega=deriv(f,omega,const_calc)

    grad_f=np.array([func_gradflow_f_flat(f[i],df[i],ddf[i],domega[i],i*dr+reg,mpi,fpi,beta) for i in range(N)])
    grad_omega=np.array([func_gradflow_omega_flat(f[i],df[i],omega[i],domega[i],ddomega[i],i*dr+reg,momega,beta) for i in range(N)])

    next_f = f + alpha_f*grad_f
    next_omega = omega + alpha_omega*grad_omega

    return next_f, next_omega

@jit(nopython=True, cache=True, fastmath=True)
def deriv(f,omega,const_calc):

    #####
    N=int(const_calc["N"]); dr=const_calc["dr"]; reg=const_calc["reg"]
    #####

    f_temp = np.array([0.0 for i in range(N+2)])
    f_temp[1:N+1]=f[0:N]
    omega_temp = np.array([0.0 for i in range(N+2)])
    omega_temp[1:N+1]=omega[0:N]

    ##### Dirichlet BC for f
    f_temp[0] = np.pi
    #    f_temp[N+1] = func_solution_f_1st(N*dr+reg)

    ##### Neumann BC for f
    #    f_temp[0] = f[0]
    f_temp[N+1] = f[N-1]

    ##### Dirichlet BC for omega
    #    omega_temp[0] = func_solution_omega(reg)
    #    omega_temp[N+1] = func_solution_omega(N*dr+reg)

    ##### Neumann BC for lfs of omega
    omega_temp[0] = omega[0]
    omega_temp[N+1] = omega[N-1]

    df=[(f_temp[i+1]-f_temp[i-1])/2.0/dr for i in range(1,N+1)]
    ddf=[(f_temp[i+1]+f_temp[i-1]-2.0*f_temp[i])/dr**2 for i in range(1,N+1)]
    domega=[(omega_temp[i+1]-omega_temp[i-1])/2.0/dr for i in range(1,N+1)]
    ddomega=[(omega_temp[i+1]+omega_temp[i-1]-2.0*omega_temp[i])/dr**2 for i in range(1,N+1)]

    return df,ddf,domega,ddomega


dict_barbeta={"barbeta_01":1.0,"barbeta_02":5.0,"barbeta_03":10.0}


#####
N=250; dr=0.016; reg=dr

# unit: powers of fm
hbarc=0.1973

fpi=0.186/hbarc
mpi=0.138/hbarc; momega=0.782/hbarc

alpha_f = 0.000001; alpha_omega = 0.000001

Niter=200000; Nsep_print=10000
#####

#####
const_phys=Dict(); const_calc=Dict()

for k,v in {"N":float(N),"dr":dr,"reg":reg,"alpha_f":alpha_f,"alpha_omega":alpha_omega}.items():
    const_calc[k]=v

for k,v in {"mpi":mpi,"momega":momega,"fpi":fpi}.items():
    const_phys[k]=v
#####

finit = np.array([func_initconf_f(i*dr+reg) for i in range(N)])
omegainit = np.array([func_initconf_omega(i*dr+reg) for i in range(N)])

dict_f={}; dict_omega={}

for key_barbeta, val_barbeta in dict_barbeta.items():

    barbeta = val_barbeta; beta = 2.0*np.pi**2*fpi/momega*barbeta
    const_phys["beta"]=beta
    print("barbeta: ",barbeta)

    #####
    f_all=[]; omega_all=[]
    f=finit; omega=omegainit
    f_all=f_all+[f]; omega_all=omega_all+[omega]
    #####

    for i in range(Niter):

        next_f, next_omega = fwd_gradflow(f,omega,const_phys,const_calc)

        if i % Nsep_print == 0:
            f_all=f_all+[next_f]; omega_all=omega_all+[next_omega]

        f=next_f; omega=next_omega

    dict_f[key_barbeta]=f_all
    dict_omega[key_barbeta]=omega_all

gradient_flow_omega-stab-Skyrmion_plot.py

import matplotlib.pyplot as plt
import matplotlib.collections as mc
import matplotlib.cm as cm

r=[(i*dr+reg) for i in range(N)]

iter_print=int(Niter/Nsep_print)

#solution_f=np.array([func_solution_f_1st(i*dr+reg) for i in range(N)])
#solution_omega=np.array([func_solution_omega(i*dr+reg) for i in range(N)])

fig = plt.figure(figsize = (10,8))
plt.rcParams["font.size"] = 18

for key_barbeta, val_barbeta in dict_barbeta.items():
    plt.plot(r, dict_f[key_barbeta][iter_print],color=cm.hsv(val_barbeta/10.0),lw=2.6,label=rf"f, $\bar\beta=${val_barbeta}")
    plt.plot(r, dict_omega[key_barbeta][iter_print],color=cm.hsv(val_barbeta/10.0),lw=1.4,label=rf"$\omega$, $\bar\beta=${val_barbeta}")

plt.plot(r,[func_initconf_f(i) for i in r],":",color="lightgreen")
plt.plot(r,[func_initconf_omega(i) for i in r],"-.",color="lightgreen")

plt.xlabel(r"$r$ [fm]")
plt.ylabel(r"$f \ [-], \omega[{\rm fm}^{-1}]$")

plt.title(rf"$f,\omega$, Niter={Niter}")
plt.grid()

plt.xlim(0,2.5)
plt.ylim(-1.0,3.2)

plt.legend(loc = 'upper right')

plt.savefig("omega_stab_Skyrmion_bbeta-1-5-10.png")
plt.show()