https://github.com/JoeMcEwen/FAST-PT
Tip revision: 68c2df21063d59bb39d5a1e5f1747e463ba00d29 authored by Jonathan Blazek on 21 September 2016, 01:19:09 UTC
updated arxiv number
updated arxiv number
Tip revision: 68c2df2
RG_STS.py
'''
Code to calcualte Renormalization group resutls using
a super time step (STS) integrator. The integrator we use is Eq. 2.9 of
'SUPER-TIME-STEPPING ACCELERATION OF
EXPLICIT SCHEMES FOR PARABOLIC PROBLEMS
Vasilios Alexiades* Genevieve Amiezy and Pierre-Alain Gremaudz '
The STS integrator works best for k_max greater than 1 and for over 500 k-grid points.
Please see the paper for more details.
J. E. McEwen (c) 2016
mcewen.24@osu.edu
'''
import numpy as np
import matplotlib.pyplot as plt
from fastpt_extr import p_window
import FASTPT
import time, sys
def RG_STS(STS_params,name,k,P,d_lambda,max,n_pad,P_window,C_window):
stage, mu, dlambda_CFL=STS_params
#save name
name='RG_STS_'+name
print('save name=', name)
# The spacing in log(k)
Delta=np.log(k[1])-np.log(k[0])
# This window function tapers the edge of the power spectrum.
# It is applied to STS
# You could change it here.
W=p_window(k,P_window[0],P_window[1])
#W=1
t1=time.time()
# windowed initial power spectrum
P_0=P*W
# standard 1-loop parameters
nu=-2
fastpt=FASTPT.FASTPT(k,nu=nu,n_pad=n_pad)
P_spt=fastpt.one_loop(P_0,C_window=C_window)
P_spt=P_0+P_spt
# initial lambda
Lambda=0
d_out=np.zeros((3,k.size+1))
d_out[0,:]=np.append(Lambda,k)
d_out[2,:]=np.append(Lambda,P_0)
d_out[1,:]=np.append(Lambda,P_spt)
dt_j=np.zeros(stage)
for j in range(stage):
arg =np.pi*(2*j-1)/(2*stage)
dt_j[j]=(dlambda_CFL)*((1.+mu)-(1.-mu)*np.cos(arg))**(-1.)
i=0
while (Lambda <= max):
t_step =0
for j in range(stage):
K=fastpt.one_loop(P,C_window=C_window)
K=K*W
P=P+ dt_j[j]*K
t_step=t_step+dt_j[j]
d_lambda=t_step
# check for failure.
if (np.any(np.isnan(P))):
print('RG flow has failed. It could be that you have not chosen a step size well.')
print('You may want to consider a smaller step size.')
print('iteration number and lambda', i, Lambda)
sys.exit()
if (np.any(np.isinf(P))):
print('RG flow has failed. It could be that you have not chosen a step size well.')
print('You may want to consider a smaller step size.')
print('iteration number and lambda', i, Lambda)
sys.exit()
#update lambda and the iteration
i=i+1
Lambda+=d_lambda
# break if the next step is already passed lambda max
if (Lambda >=max):
break
#print('lambda', Lambda )
# update data for saving
d_update=np.append(Lambda,P)
d_out=np.row_stack((d_out,d_update))
# set to True to check at each step
if(False):
ax=plt.subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('k')
ax.plot(k,P)
ax.plot(k,P_0, color='red')
plt.grid()
plt.show()
''' Since STS method computes its own Delta Lambda.
So,this routine will not end exaclty at max.
Therefore, take Euler steps to reach the end
'''
last_step=np.absolute(Lambda-d_lambda-max)
#print('Lambda', Lambda-d_lambda)
#print('last step', last_step)
k1=fastpt.one_loop(P,C_window=C_window)
k1=k1*W
k2=fastpt.one_loop(k1*last_step/2.+ P,C_window=C_window)
k2=k2*W
k3=fastpt.one_loop(k2*last_step/2. +P,C_window=C_window)
k3=k3*W
k4=fastpt.one_loop(k3*last_step +P,C_window=C_window)
k4=k4*W
# full step
P=P+1/6.*(k1+2*k2+2*k3+k4)*last_step
Lambda=Lambda-d_lambda+last_step
# update data for saving
d_update=np.append(Lambda,P)
d_out=np.row_stack((d_out,d_update))
# save the data
t2=time.time()
print('time to run seconds', t2-t1)
print('time to run minutes', (t2-t1)/60.)
print('iteration', i )
print('save name', name)
print('shape', d_out.shape)
np.save(name,d_out)
print('number of iterations and output shape', i, d_out.shape)
# return last step
return P
if __name__ == "__main__":
# set the STS parameters here
# this combo seems to work well for k_max=10, 2000 grid points
# if you encounter instabilities you may want to fiddle with these values.
stage=10
mu=.1
# this is \Delta \lambda_{CFL}
dlambda_CFL=1e-3
STS_params=[stage, mu, dlambda_CFL]
V=sys.version_info[0]
if (V < 3):
from ConfigParser import SafeConfigParser
if (V >=3 ):
from configparser import SafeConfigParser
parser = SafeConfigParser()
name='kmax10_example.ini'
parser.read(name)
k_max=parser.getfloat('floats', 'k_max')
k_min=parser.getfloat('floats', 'k_min')
step=parser.getfloat('floats', 'step')
max=parser.getfloat('floats', 'max')
P_right=parser.getfloat('floats', 'P_w_right')
P_left=parser.getfloat('floats', 'P_w_left')
C_window=parser.getfloat('floats', 'C_window')
n_pad=parser.getint('integers', 'n_pad')
down_sample=parser.getint('integers', 'down_sample')
read_name=parser.get('files', 'in_file')
name=parser.get('files', 'out_file')
d=np.loadtxt(read_name) # load data
k=d[:,0]
P=d[:,1]
id=np.where( (k >= k_min) & (k <= k_max) )[0]
k=k[id]
P=P[id]
k=k[::down_sample]
P=P[::down_sample]
# if your array is not even in size, FAST-PT will not work-
# trim if so.
if (k.size % 2 != 0):
k=k[:-1]
P=P[:-1]
print('Details of run.')
print('save name :', name)
print('k min and max:', k_min, k_max)
print('grid size : ', k.size)
print('d log k: ', np.log(k[1])-np.log(k[0]) )
print('down sample factor:', down_sample)
P_window=np.array([P_left,P_right])
P_rg=RG_STS(STS_params,name,k,P,step,max,n_pad,P_window,C_window)
ax=plt.subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('k')
ax.set_ylabel(r'$P(k)$', size=30)
ax.set_xlabel(r'$k$', size=30)
ax.plot(k,P, label='linear power')
ax.plot(k,P_rg, label='RG' )
plt.legend(loc=3)
plt.grid()
plt.show()