Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61096473
active_remeshing.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, May 4, 12:41
Size
5 KB
Mime Type
text/x-python
Expires
Mon, May 6, 12:41 (2 d)
Engine
blob
Format
Raw Data
Handle
17464549
Attached To
R6746 RationalROMPy
active_remeshing.py
View Options
import
numpy
as
np
from
pickle
import
load
from
matplotlib
import
pyplot
as
plt
from
active_remeshing_engine
import
ActiveRemeshingEngine
from
rrompy.reduction_methods
import
(
RationalInterpolantGreedyPivotedNoMatch
as
RIGPNM
,
RationalInterpolantGreedyPivotedPoleMatch
as
RIGPG
)
from
rrompy.parameter.parameter_sampling
import
(
QuadratureSampler
as
QS
,
SparseGridSampler
as
SGS
)
zs
,
ts
=
[
0.
,
100.
],
[
0.
,
.
5
]
z0
,
t0
,
n
=
np
.
mean
(
zs
),
np
.
mean
(
ts
),
150
solver
=
ActiveRemeshingEngine
(
z0
,
t0
,
n
)
solver
.
cutOffPolesRMinRel
,
solver
.
cutOffPolesRMaxRel
=
-
2.5
,
2.5
solver
.
cutOffPolesIMin
,
solver
.
cutOffPolesIMax
=
-.
01
,
.
01
mus
=
[[
z0
,
ts
[
0
]],
[
z0
,
ts
[
1
]]]
for
mu
in
mus
:
u
=
solver
.
solve
(
mu
,
return_state
=
True
)[
0
]
Y
=
solver
.
applyC
(
u
,
mu
)[
0
]
_
=
solver
.
plot
(
u
,
what
=
"REAL"
,
name
=
"u(z={}, t={})"
.
format
(
*
mu
),
is_state
=
True
,
figsize
=
(
12
,
4
))
print
(
"Y(z={}, t={}) = {} (solution average)"
.
format
(
*
mu
,
np
.
real
(
Y
)))
fighandles
=
[]
with
open
(
"./active_remeshing_hf_samples.pkl"
,
"rb"
)
as
f
:
zspace
,
tspace
,
Yex
=
load
(
f
)
# zspace, tspace = np.linspace(*zs, 200), np.linspace(*ts, 50)
# Yex = [[solver.solve([z, t]) for t in tspace] for z in zspace]
# (from a ~2h45m simulation on one node of the EPFL Helvetios cluster)
radius2Err
=
np
.
mean
(
np
.
abs
(
Yex
)
**
2.
)
YCmin
,
YCmax
=
np
.
quantile
(
Yex
,
.
05
),
np
.
quantile
(
Yex
,
.
95
)
YexC
=
np
.
clip
(
Yex
,
YCmin
,
YCmax
)
approx
=
[]
for
match
in
range
(
2
):
params
=
{
"POD"
:
True
,
"S"
:
5
,
"greedyTol"
:
1e-4
,
"nTestPoints"
:
500
,
"polybasis"
:
"LEGENDRE"
,
"samplerTrainSet"
:
QS
(
zs
,
"UNIFORM"
),
"samplerPivot"
:
QS
(
zs
,
"CHEBYSHEV"
),
"SMarginal"
:
5
,
"samplerMarginal"
:
SGS
(
ts
),
"errorEstimatorKind"
:
"LOOK_AHEAD_OUTPUT"
}
if
match
:
print
(
"
\n
Testing output-based matching."
)
params
[
"matchingWeight"
]
=
1.
params
[
"matchingShared"
]
=
.
75
params
[
"polybasisMarginal"
]
=
"PIECEWISE_LINEAR_UNIFORM"
algo
=
RIGPG
else
:
print
(
"
\n
Testing matching-free approach."
)
algo
=
RIGPNM
approx
+=
[
algo
([
0
],
solver
,
mu0
=
[
z0
,
t0
],
approxParameters
=
params
,
verbosity
=
5
)]
if
match
:
approx
[
match
]
.
setTrainedModel
(
approx
[
0
])
else
:
approx
[
match
]
.
setupApprox
()
verb
=
approx
[
match
]
.
verbosity
verbTM
=
approx
[
match
]
.
trainedModel
.
verbosity
approx
[
match
]
.
verbosity
,
approx
[
match
]
.
trainedModel
.
verbosity
=
0
,
0
for
j
,
t
in
enumerate
(
tspace
):
out
=
approx
[
match
]
.
getApprox
(
np
.
pad
(
zspace
.
reshape
(
-
1
,
1
),
[(
0
,
0
),
(
0
,
1
)],
"constant"
,
constant_values
=
t
))
pls
=
approx
[
match
]
.
getPoles
([
None
,
t
])
pls
[
np
.
abs
(
np
.
imag
(
pls
))
>
1e-5
]
=
np
.
nan
if
j
==
0
:
Ys
=
np
.
empty
((
len
(
zspace
),
len
(
tspace
)))
poles
=
np
.
empty
((
len
(
tspace
),
len
(
pls
)))
Ys
[:,
j
]
=
out
.
re
.
data
if
len
(
pls
)
>
poles
.
shape
[
1
]:
poles
=
np
.
pad
(
poles
,
[(
0
,
0
),
(
0
,
len
(
pls
)
-
poles
.
shape
[
1
])],
"constant"
,
constant_values
=
np
.
nan
)
poles
[
j
,
:
len
(
pls
)]
=
np
.
real
(
pls
)
approx
[
match
]
.
verbosity
=
verb
approx
[
match
]
.
trainedModel
.
verbosity
=
verbTM
YsC
=
np
.
clip
(
Ys
,
YCmin
,
YCmax
)
err
=
(
np
.
abs
(
Yex
-
YsC
)
/
(
np
.
abs
(
Yex
)
**
2.
+
radius2Err
)
**
.
5
/
(
np
.
abs
(
Ys
)
**
2.
+
radius2Err
)
**
.
5
)
fighandles
+=
[
plt
.
figure
(
figsize
=
(
15
,
5
))]
ax1
=
fighandles
[
-
1
]
.
add_subplot
(
1
,
2
,
1
)
ax2
=
fighandles
[
-
1
]
.
add_subplot
(
1
,
2
,
2
)
if
match
:
ax1
.
plot
(
poles
,
tspace
)
else
:
ax1
.
plot
(
poles
,
tspace
,
"k."
)
ax1
.
set_ylim
(
ts
)
ax1
.
set_xlabel
(
"z"
)
ax1
.
set_ylabel
(
"t"
)
ax1
.
grid
()
if
match
:
ax2
.
plot
(
poles
,
tspace
)
else
:
ax2
.
plot
(
poles
,
tspace
,
"k."
)
for
mm
in
approx
[
match
]
.
musMarginal
:
ax2
.
plot
(
zs
,
[
mm
[
0
,
0
]]
*
2
,
"k--"
,
linewidth
=
1
)
ax2
.
set_xlim
(
zs
)
ax2
.
set_ylim
(
ts
)
ax2
.
set_xlabel
(
"z"
)
ax2
.
set_ylabel
(
"t"
)
ax2
.
grid
()
plt
.
show
()
print
(
"Approximate poles"
)
fighandles
+=
[
plt
.
figure
(
figsize
=
(
15
,
5
))]
ax1
=
fighandles
[
-
1
]
.
add_subplot
(
1
,
2
,
1
)
ax2
=
fighandles
[
-
1
]
.
add_subplot
(
1
,
2
,
2
)
p
=
ax1
.
contourf
(
np
.
repeat
(
zspace
.
reshape
(
-
1
,
1
),
len
(
tspace
),
axis
=
1
),
np
.
repeat
(
tspace
.
reshape
(
1
,
-
1
),
len
(
zspace
),
axis
=
0
),
YsC
,
vmin
=
YCmin
,
vmax
=
YCmax
,
levels
=
np
.
linspace
(
YCmin
,
YCmax
,
50
))
plt
.
colorbar
(
p
,
ax
=
ax1
)
ax1
.
set_xlabel
(
"z"
)
ax1
.
set_ylabel
(
"t"
)
ax1
.
grid
()
p
=
ax2
.
contourf
(
np
.
repeat
(
zspace
.
reshape
(
-
1
,
1
),
len
(
tspace
),
axis
=
1
),
np
.
repeat
(
tspace
.
reshape
(
1
,
-
1
),
len
(
zspace
),
axis
=
0
),
YexC
,
vmin
=
YCmin
,
vmax
=
YCmax
,
levels
=
np
.
linspace
(
YCmin
,
YCmax
,
50
))
ax2
.
set_xlabel
(
"z"
)
ax2
.
set_ylabel
(
"t"
)
ax2
.
grid
()
plt
.
colorbar
(
p
,
ax
=
ax2
)
plt
.
show
()
print
(
"Approximate and exact output
\n
"
)
fighandles
+=
[
plt
.
figure
(
figsize
=
(
9
,
6
))]
ax1
=
fighandles
[
-
1
]
.
add_subplot
(
1
,
1
,
1
)
p
=
ax1
.
contourf
(
np
.
repeat
(
zspace
.
reshape
(
-
1
,
1
),
len
(
tspace
),
axis
=
1
),
np
.
repeat
(
tspace
.
reshape
(
1
,
-
1
),
len
(
zspace
),
axis
=
0
),
np
.
log10
(
err
),
vmin
=
-
10
,
vmax
=
0
,
levels
=
np
.
linspace
(
-
10
,
0
,
50
))
plt
.
colorbar
(
p
,
ax
=
ax1
)
ax1
.
set_xlabel
(
"z"
)
ax1
.
set_ylabel
(
"t"
)
ax1
.
grid
()
plt
.
show
()
print
(
"Output error (log-chordal)
\n
"
)
Event Timeline
Log In to Comment