Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69320616
test_surface.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Jul 1, 06:39
Size
3 KB
Mime Type
text/x-python
Expires
Wed, Jul 3, 06:39 (2 d)
Engine
blob
Format
Raw Data
Handle
18696373
Attached To
rTAMAAS tamaas
test_surface.py
View Options
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
Log In to Comment