Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_autocorrelation.py
611b0aabdd07features/netcdf-mpi

Last Change
Raw File

test_autocorrelation.py
View Options

#!/usr/bin/python
# -*- coding: utf-8 -*-
##*
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
#
# @section LICENSE
#
# Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
# Solides)
#
# Tamaas is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
#
#
################################################################
from tamaas import *
import numpy as np
from scipy.special import jv
import scipy.integrate as integrate
from numpy.fft import *
import sys, imp
def spectrum(r, kl, kr, ks, C, H):
if r > ks:
return 0
elif r > kr:
return C * (r/float(kr))**(-2.*(H+1))
elif r > kl:
return C
else:
return 0
# n is normalization parameter
def integrate_bessel(x, spectrum, n):
result = 2 * np.pi * integrate.quad(lambda y: y*jv(0, y*x) * spectrum(n*y), 0, np.inf)[0]
return result
# n is normalization parameter
def integrate_bessel(x, spectrum, n):
result = 2 * np.pi * integrate.quad(lambda y: y*jv(0, y*x) * spectrum(n*y), 0, np.inf)[0]
return result
def main():
grid_size = 256
lambda_l = 8.
lambda_r = 8.
lambda_s = 1.
L = 4 * lambda_l
kl = int(L / lambda_l)
kr = int(L / lambda_r)
ks = int(L / lambda_s)
# Computing ACF from integral with Bessel function
l = grid_size/2 # points on curve
x = np.linspace(0, L/2, l)
bessel_acf = np.zeros(l)
integrate_bessel_v = np.vectorize(integrate_bessel, otypes=[np.float])
bessel_acf = (L / (2*np.pi) )**2 *integrate_bessel_v(x, lambda x: spectrum(x,
kl, kr, ks,
1., 0.8), L / (2*np.pi))
nsamples = 10000
etages=np.zeros((grid_size,grid_size,nsamples))
# Computing ACF with Hu & Tonder algo reimplemented
for i in range(nsamples):
#generate surface
SG = SurfaceGeneratorFilterFFT()
SG.getGridSize().assign(grid_size)
SG.getHurst().assign(0.8)
SG.getQ0().assign(kl);
SG.getQ1().assign(kr);
SG.getQ2().assign(ks);
SG.getRandomSeed().assign(i);
SG.Init()
a = SG.buildSurface()
etages[:,:,i]=a
cov=np.zeros(grid_size/2)
for i in np.arange(grid_size/2,grid_size,1):
cov[i-grid_size/2]=np.sum((etages[grid_size/2,i,:]-np.mean(etages[grid_size/2,i,:]))*(etages[grid_size/2,grid_size/2,:]-np.mean(etages[grid_size/2,grid_size/2,:])))/nsamples
error=np.sum((cov-bessel_acf)*(cov-bessel_acf))/np.sum(bessel_acf*bessel_acf)
print(error)
if error > 1e-1 :
return 1
return 0
if __name__ == "__main__":
sys.exit(main())
c4science · Help