Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/examples/surface_1d.py
f52d1730a909features/netcdf-mpi

Last Change
Raw File

surface_1d.py
View Options

#!/usr/bin/python
# -*- coding: utf-8 -*-
#
# @author Lucas Frérot <lucas.frerot@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 __future__ import division
import tamaas as tm
from numpy.fft import fft, fftshift, fftfreq
import matplotlib.pyplot as plt
tm.initialize()
def draw_triangle(pos, w, ax, style, alpha=1.5, base=10):
h = alpha * w
ax.loglog([base**pos[0], base**(pos[0]+w)],
[base**pos[1], base**pos[1]], style)
ax.loglog([base**pos[0], base**pos[0]],
[base**pos[1], base**(pos[1]+h)], style)
ax.loglog([base**(pos[0]+w), base**pos[0]],
[base**pos[1], base**(pos[1]+h)], style)
# ax.text(base**(pos[0]+w/3), base**(pos[1]+h/3), "$-{}$".format(alpha),
# horizontalalignment='center',
# verticalalignment='center')
ax.text(base**(pos[0]+w/2), base**(pos[1]-0.2), '$1$',
horizontalalignment='center',
verticalalignment='top')
ax.text(base**(pos[0]-0.15), base**(pos[1]+h/2), '${}$'.format(alpha),
horizontalalignment='right',
verticalalignment='center')
iso = tm.Isopowerlaw1D()
helper = tm.ParamHelper(iso)
params = {
"Q0": 16,
"Q1": 16,
"Q2": 512,
"Hurst": 0.8
}
helper.set(params)
gen = tm.SurfaceGeneratorFilter1D()
gen.setFilter(iso)
gen.setSizes([2048])
# gen.setRandomSeed(0) # uncomment for fixed seed
surf = gen.buildSurface()
plt.plot(surf)
surf_fft = fft(surf)
psd = surf_fft*surf_fft.conj()
psd = fftshift(psd)
plt.figure()
freqs = fftshift(fftfreq(psd.size, d=1/psd.size))
plt.plot(freqs[psd.size/2+1:], psd.real[psd.size/2+1:])
draw_triangle([5, -4], 2, plt.gca(), '-k', 2*(params["Hurst"]+1), 2)
plt.yscale('log')
plt.xscale('log', basex=2)
plt.ylim([10**(-2), 10**8])
plt.show()
tm.finalize()
c4science · Help