Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_autocorrelation.py
dd37e82c03d3features/netcdf-mpi

Last Change
Raw File

test_autocorrelation.py
View Options

# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program 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 Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
from __future__ import print_function
import tamaas as tm
import numpy as np
from scipy.special import jv
import scipy.integrate as integrate
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
def test_autocorrelation(tamaas_fixture):
tm.initialize()
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
n = grid_size // 2 # points on curve
x = np.linspace(0, L/2, n)
bessel_acf = np.zeros(n)
integrate_bessel_v = np.vectorize(integrate_bessel, otypes=[tm.dtype])
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 = 1000
etages = np.zeros((grid_size, grid_size, nsamples))
# Computing ACF with Hu & Tonder algo reimplemented
spec = tm.Isopowerlaw2D()
spec.q0, spec.q1, spec.q2 = kl, kr, ks
spec.hurst = 0.8
SG = tm.SurfaceGeneratorFilter2D([grid_size] * 2)
SG.spectrum = spec
for i in range(nsamples):
SG.random_seed = i
a = SG.buildSurface()
etages[:, :, i] = a[:, :]
cov = np.zeros(n)
for i in np.arange(n, grid_size, 1):
cov[i-n] = np.sum((etages[n, i, :]-np.mean(etages[n, i, :])) *
(etages[n, n, :]-np.mean(etages[n, n, :])))/nsamples
error = np.sum((cov-bessel_acf)*(cov-bessel_acf)) \
/ np.sum(bessel_acf*bessel_acf)
print(error)
assert error < 1e-1
if __name__ == "__main__":
test_autocorrelation()
c4science · Help