gradient_flow_vortex.py

import numpy as np

def func_initial_vortex(x):
    return np.tanh(0.3*x)

def func_correctsolution_vortex(x):
    return np.tanh(0.54795*x)

def func_gradE_vortex(r,f,df,ddf,m):
    return (ddf+df/r-m**2*f/r**2+(1.0-f**2)*f)

def fwd_gradflow_vortex(f,dr,alpha,N,m,reg,func_gradE):

    df,ddf = deriv_vortex(f,dr,N,func_initial_vortex)

    grad_f=np.array([func_gradE(i*dr+reg,f[i],df[i],ddf[i],m) for i in range(N)])

    nextgrad_f=f+alpha*grad_f
    nextgrad_f[0]=0.0; nextgrad_f[N-1]=1.0

    return nextgrad_f,df,ddf,grad_f

def deriv_vortex(f,dr,N,func):

    # Dirichlet BC
    ftemp = np.insert(f,0,func(-1.0*dr+reg))
    ftemp = np.insert(ftemp,N+1,func(N*dr+reg))

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

    return df,ddf


###############
N=200; dr=0.05
reg=0.001

m=1
Nsep_print=100

alpha=0.001

Niter=4000
###############

f_all=[]
#df_all=[]; ddf_all=[]; grad_f_all=[]

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

f=finit

f_all=f_all+[f]

for i in range(Niter):

    nextgrad_f,df,ddf,grad_f = fwd_gradflow_vortex(f,dr,alpha,N,m,reg,func_gradE_vortex)

    if i%Nsep_print == 0:
        f_all=f_all+[nextgrad_f]

#    df_all=df_all+[df]
#    ddf_all=ddf_all+[ddf]
#    grad_f_all=grad_f_all+[grad_f]

    f=nextgrad_f

gradient_flow_vortex_plot.py

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

iter_plot=int(Niter/Nsep_print)
correctsol=np.array([func_correctsolution_vortex(i*dr+reg) for i in range(N)])

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

fig = plt.figure(figsize = (10,7))
plt.rcParams["font.size"] = 16
#fig = plt.figure(figsize = (10,7), facecolor='lightblue')

plt.plot(r, finit,color="black",label="initial configuration")
plt.plot(r, f_all[5], label="U, iteration: 500", linestyle="--", color=cm.hsv(0.4))
plt.plot(r, f_all[10], label="U, iteration: 1000", linestyle="--", color=cm.hsv(0.5))
plt.plot(r, f_all[20], label="U, iteration: 2000", linestyle="--", color=cm.hsv(0.6))
plt.plot(r, f_all[30], label="U, iteration: 3000", linestyle="--", color=cm.hsv(0.7))
plt.plot(r, f_all[iter_plot], label="U, iteration: 4000", linestyle="--", color=cm.hsv(0.8))
plt.plot(r, correctsol,color="red",label="correct solution")
#plt.plot(r, f_all[Niter],color="red",label="grad_f")
#plt.plot(r, omega_all[Niter],color="blue",label="grad_omega")

plt.xlabel(r"$r$")
plt.ylabel(r"$U$")
plt.grid()

plt.xlim(-0.01,10)
#plt.ylim(0,1.01)

plt.legend(loc = 'lower right')

plt.show()