Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F70279329
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
Sat, Jul 6, 02:32
Size
4 KB
Mime Type
text/x-python
Expires
Mon, Jul 8, 02:32 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18817502
Attached To
rTAMAAS tamaas
test_surface.py
View Options
# -*- coding: utf-8 -*-
#
# Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
# Copyright (©) 2020-2022 Lucas Frérot
#
# 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
,
division
import
pytest
import
numpy
as
np
import
tamaas
as
tm
from
numpy
import
pi
integrate
=
pytest
.
importorskip
(
'scipy.integrate'
)
special
=
pytest
.
importorskip
(
'scipy.special'
)
jv
=
special
.
jv
E
=
special
.
ellipe
def
spectrum
(
r
,
isopowerlaw
):
kl
,
kr
,
ks
=
isopowerlaw
.
q0
,
isopowerlaw
.
q1
,
isopowerlaw
.
q2
H
=
isopowerlaw
.
hurst
C
=
1.
if
r
>
ks
:
return
0
elif
r
>
kr
:
return
C
*
(
r
/
float
(
kr
))
**
(
-
2.
*
(
H
+
1
))
elif
r
>
kl
:
return
C
else
:
return
0
def
integrate_bessel
(
x
,
spectrum
,
n
):
return
2
*
pi
*
integrate
.
quad
(
lambda
y
:
y
*
jv
(
0
,
y
*
x
)
*
spectrum
(
n
*
y
),
0
,
np
.
inf
,
limit
=
200
)[
0
]
integrate_bessel
=
np
.
vectorize
(
integrate_bessel
,
otypes
=
[
tm
.
dtype
])
@pytest.fixture
(
scope
=
"module"
,
params
=
[
tm
.
SurfaceGeneratorFilter2D
,
tm
.
SurfaceGeneratorRandomPhase2D
])
def
stats
(
request
):
iso
=
tm
.
Isopowerlaw2D
()
iso
.
q0
=
16
iso
.
q1
=
32
iso
.
q2
=
64
iso
.
hurst
=
0.8
N
=
243
sg
=
request
.
param
()
sg
.
spectrum
=
iso
sg
.
shape
=
[
N
,
N
]
analytical_stats
=
{
"rms_heights"
:
iso
.
rmsHeights
(),
"rms_slopes_spectral"
:
iso
.
rmsSlopes
(),
"moments"
:
iso
.
moments
(),
"alpha"
:
iso
.
alpha
()
}
stats
=
{
k
:
[]
for
k
in
analytical_stats
}
acf
=
np
.
zeros
(
sg
.
shape
,
dtype
=
tm
.
dtype
)
realizations
=
1000
for
i
in
range
(
realizations
):
sg
.
random_seed
=
i
s
=
sg
.
buildSurface
()
stats
[
'rms_heights'
]
.
append
(
tm
.
Statistics2D
.
computeRMSHeights
(
s
))
stats
[
'rms_slopes_spectral'
]
.
append
(
tm
.
Statistics2D
.
computeSpectralRMSSlope
(
s
))
stats
[
'moments'
]
.
append
(
tm
.
Statistics2D
.
computeMoments
(
s
))
acf
+=
tm
.
Statistics2D
.
computeAutocorrelation
(
s
)
/
realizations
L
=
1
*
sg
.
spectrum
.
q0
n
=
sg
.
shape
[
0
]
//
2
x
=
np
.
linspace
(
0
,
L
/
2
,
n
)
true_acf
=
(
L
/
(
2
*
pi
))
**
2
\
*
integrate_bessel
(
x
,
lambda
x
:
spectrum
(
x
,
sg
.
spectrum
),
L
/
(
2
*
pi
))
return
analytical_stats
,
stats
,
acf
,
true_acf
@pytest.mark.parametrize
(
'key'
,
[
'rms_heights'
,
'rms_slopes_spectral'
])
def
test_statistics
(
stats
,
key
):
analytical_stats
,
stats
,
_
,
_
=
stats
mu
=
np
.
mean
(
stats
[
key
])
error
=
np
.
abs
(
mu
-
analytical_stats
[
key
])
/
analytical_stats
[
key
]
assert
error
<
3e-3
def
test_autocorrelation
(
stats
):
from
tamaas.utils
import
radial_average
_
,
_
,
acf
,
true_acf
=
stats
x
,
y
=
(
np
.
linspace
(
-
1
,
1
,
acf
.
shape
[
i
])
for
i
in
range
(
2
))
r
=
np
.
linspace
(
0
,
1
,
true_acf
.
shape
[
0
])
theta
=
np
.
linspace
(
0
,
2
*
np
.
pi
,
endpoint
=
False
)
r_acf
=
radial_average
(
x
,
y
,
np
.
fft
.
fftshift
(
acf
),
r
,
theta
,
endpoint
=
False
)
acf_error
=
np
.
sum
((
r_acf
-
true_acf
)
**
2
)
/
np
.
sum
(
true_acf
**
2
)
assert
acf_error
<
2e-3
def
test_graph_area
():
N
=
100
x
=
np
.
linspace
(
0
,
1
,
N
,
endpoint
=
False
)
f
=
np
.
sin
(
2
*
np
.
pi
*
x
)
df
=
2
*
np
.
pi
*
np
.
cos
(
2
*
np
.
pi
*
x
)
num_area
=
np
.
sqrt
(
1
+
df
**
2
)
.
mean
()
area
=
2
*
np
.
sqrt
(
1
+
4
*
np
.
pi
**
2
)
/
np
.
pi
*
E
(
4
*
np
.
pi
**
2
/
(
1
+
4
*
np
.
pi
**
2
))
fourier_area
=
tm
.
Statistics1D
.
graphArea
(
f
)
assert
np
.
abs
(
num_area
-
fourier_area
)
/
num_area
<
1e-10
assert
np
.
abs
(
area
-
fourier_area
)
/
area
<
2e-10
Event Timeline
Log In to Comment