https://github.com/JoeMcEwen/FAST-PT
Revision c56761670c0e4c154b7f2794c79c06c840ba5b0a authored by Jonathan Blazek on 01 March 2019, 22:19:49 UTC, committed by Jonathan Blazek on 01 March 2019, 22:19:49 UTC
1 parent 137418e
Tip revision: c56761670c0e4c154b7f2794c79c06c840ba5b0a authored by Jonathan Blazek on 01 March 2019, 22:19:49 UTC
merging changes from develop branch to prepare for merge with dev. develop will be deprecated.
merging changes from develop branch to prepare for merge with dev. develop will be deprecated.
Tip revision: c567616
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()
Computing file changes ...