1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
"""
Viscous free fall: metamodel of a field function
================================================
"""
# %%
#
# In this example, we present how to create the metamodel of a field function.
# This example considers the :ref:`free fall model <use-case-viscous-fall>`.
# We create the metamodel automatically using :class:`openturns.experimental.PointToFieldFunctionalChaosAlgorithm`
# and then also with a manual approach:
# We first compute the Karhunen-Loève decomposition of a sample of trajectories.
# Then we create a create a polynomial chaos which takes the inputs and returns
# the KL decomposition modes as outputs.
# Finally, we create a metamodel by
# combining the KL decomposition and the polynomial chaos.
# %%
# Define the model
# ----------------
# %%
import openturns as ot
import openturns.experimental as otexp
import openturns.viewer as otv
from openturns.usecases import viscous_free_fall
ot.Log.Show(ot.Log.NONE)
# %%
# Load the viscous free fall example.
vff = viscous_free_fall.ViscousFreeFall()
distribution = vff.distribution
model = vff.model
# %%
# Generate a training sample.
size = 2000
ot.RandomGenerator.SetSeed(0)
inputSample = distribution.getSample(size)
outputSample = model(inputSample)
# %%
# Compute the global metamodel
algo = otexp.PointToFieldFunctionalChaosAlgorithm(
inputSample, outputSample, distribution
)
algo.run()
result = algo.getResult()
metaModel = result.getPointToFieldMetaModel()
# %%
# Validate the metamodel
# ----------------------
# %%
# Create a validation sample.
# %%
size = 10
validationInputSample = distribution.getSample(size)
validationOutputSample = model(validationInputSample)
# %%
graph = validationOutputSample.drawMarginal(0)
graph.setColors(["red"])
graph2 = metaModel(validationInputSample).drawMarginal(0)
graph2.setColors(["blue"])
graph.add(graph2)
graph.setTitle("Model/metamodel comparison")
graph.setXTitle(r"$t$")
graph.setYTitle(r"$z$")
view = otv.View(graph)
# %%
# We see that the blue trajectories (i.e. the metamodel) are close to the red
# trajectories (i.e. the validation sample).
# This shows that the metamodel is quite accurate.
# However, we observe that the trajectory singularity that occurs when the object
# touches the ground (i.e. when :math:`z` is equal to zero), makes the metamodel less accurate.
# %%
# Sensitivity analysis
# --------------------
# %%
# Compute the sensitivity indices
sensitivity = otexp.FieldFunctionalChaosSobolIndices(result)
s1 = sensitivity.getFirstOrderIndices()
st = sensitivity.getTotalOrderIndices()
# %%
# We can notice that `v0` and `m` are the most influencial parameters
# and that there are almost no interactions (total indices being close to first order indices)
print(s1, st)
# %%
# Draw the sensitivity indices
graph = sensitivity.draw()
view = otv.View(graph)
# %%
# Manual approach
# ---------------
# %%
# Step 1: compute the KL decomposition of the output
# %%
algo = ot.KarhunenLoeveSVDAlgorithm(outputSample, 1.0e-6)
algo.run()
klResult = algo.getResult()
scaledModes = klResult.getScaledModesAsProcessSample()
# %%
graph = scaledModes.drawMarginal(0)
graph.setTitle("KL modes")
graph.setXTitle(r"$t$")
graph.setYTitle(r"$z$")
view = otv.View(graph)
# %%
# We create the lifting function which takes coefficients of the Karhunen-Loève (KL) modes as inputs and returns the trajectories.
klLiftingFunction = ot.KarhunenLoeveLifting(klResult)
# %%
# The `project` method computes the projection of the output sample (i.e. the trajectories) onto the KL modes.
outputSampleChaos = klResult.project(outputSample)
# %%
# Step 2: compute the metamodel of the KL modes
# We create a polynomial chaos metamodel which takes the input sample and returns the KL modes.
algo = ot.FunctionalChaosAlgorithm(inputSample, outputSampleChaos, distribution)
algo.run()
chaosMetamodel = algo.getResult().getMetaModel()
# %%
# The final metamodel is a composition of the KL lifting function and the polynomial chaos metamodel.
# We combine these two functions using the :class:`~openturns.PointToFieldConnection` class.
# %%
metaModel = ot.PointToFieldConnection(klLiftingFunction, chaosMetamodel)
# %%
# Reset ResourceMap
ot.ResourceMap.Reload()
# %%
otv.View.ShowAll()
|
