Page MenuHomec4science
Files F61895969

test_surface.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 9, 16:04

test_surface.py
View Options

#!/usr/bin/env 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 print_function, division
import sys
import numpy as np
import tamaas as tm
def compute_mse(values, mean):
return np.mean((values - mean)**2) / values.size
def main():
tm.initialize(1)
# Spectrum parameters in wavelengths
disctretization_size = 1.
lambda_s = 8.
lambda_r = 64.
lambda_l = 64.
box_size = 512
# Independent parameters
Hurst = 0.8
hrms = 1.
# Spectrum parameters in wavenumbers
k_l = int(box_size // lambda_l)
k_r = int(box_size // lambda_r)
k_s = int(box_size // lambda_s)
N = int(box_size // disctretization_size)
# Analytical moments, cf. Yastrebov et al. (2015)
# "From infinitesimal to full contact between rough surfaces: Evolution
# of the contact area", appendix A
xi = k_l / k_r
zeta = k_s / k_r
H = Hurst
T = {0: 2*np.pi,
2: np.pi,
4: 3*np.pi/4}
m = {}
for q in [0, 2, 4]:
m[q] = T[q] * k_r**(q-2*H)*((1-xi**(q+2))/(q+2) + (zeta**(q-2*H)-1)/(q-2*H))
# Some spectral parameters
alpha = m[0]*m[4]/m[2]**2
rms_slopes = np.sqrt(np.pi * m[2] / 2)
# Output setup
realizations = 100
hrms_out = np.zeros(realizations)
alpha_out = np.zeros(realizations)
rms_slopes_out = {"spectral": np.zeros(realizations),
"geometric": np.zeros(realizations)}
for i in range(realizations):
# Surface generator setup
sg = tm.SurfaceGeneratorFilterFFT()
sg.getGridSize().assign(N)
sg.getHurst().assign(Hurst)
sg.getRMS().assign(hrms)
sg.getQ0().assign(k_l)
sg.getQ1().assign(k_r)
sg.getQ2().assign(k_s)
sg.getRandomSeed().assign(i)
sg.Init()
# Generating surface
surface = sg.buildSurface() / sg.analyticalRMS()
# Computing RMS heights
# We know surface.mean() = 0, so we get an unbiased estimator here
hrms_out[i] = np.sqrt(1 / (N**2-1) * np.sum(surface**2))
# Computing spectrum bandwidth
surface /= box_size**2 # Normalize so that the computed PSD does not contain the box size
moments = tm.SurfaceStatistics.computeMoments(surface)
m0 = hrms_out[i]**2 / box_size**4
m2 = moments[0]
m4 = moments[1]
alpha_out[i] = m0*m4/m2**2
# Computing RMS slopes
rms_slopes_out["spectral"][i] = tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
rms_slopes_out["geometric"][i] = tm.SurfaceStatistics.computeRMSSlope(surface)
# print(rms_slopes_out)
# print(rms_slopes)
# Computing Mean Square Error
print("""========= Analytical values =========
hrms = {}
alpha = {}
m0, m2, m4 = {}
rms_slopes = {}""".format(hrms, alpha, m, rms_slopes))
print("========= MSE Calculations =========")
data_association = [
(hrms_out, hrms, "hrms"),
(alpha_out, alpha, "alpha"),
(rms_slopes_out["spectral"], rms_slopes, "RMS Slopes Spectral"),
(rms_slopes_out["geometric"], rms_slopes, "RMS Slopes Geometric")
]
# Computing and printing MSE values
mse_vals = []
for out, analytical, string in data_association:
mse = np.sqrt(compute_mse(out, analytical)) / analytical
print("Normalized error (MSE) for {}: {}".format(string, mse))
mse_vals.append(mse)
# Checking tolerances
tols = [5e-3, 5e-3, 1e-1, 1e-1]
for mse, tol in zip(mse_vals, tols):
if mse >= tol:
print("{} >= {}".format(mse, tol))
return 1
tm.finalize()
return 0
if __name__ == "__main__":
sys.exit(main())

Event Timeline

c4science · Help