1""" 2Viscous free fall: metamodel of a field function 3================================================ 4""" 5# %% 6# 7# In this example, we present how to create the metamodel of a field function. This examples considers the :ref:`free fall model <use-case-viscous-fall>`. 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. 8 9# %% 10# Define the model 11# ---------------- 12 13# %% 14from __future__ import print_function 15import openturns as ot 16import numpy as np 17import openturns.viewer as viewer 18from matplotlib import pylab as plt 19ot.Log.Show(ot.Log.NONE) 20 21# %% 22# We first define the time grid associated with the model. 23 24# %% 25tmin = 0.0 # Minimum time 26tmax = 12. # Maximum time 27gridsize = 100 # Number of time steps 28mesh = ot.IntervalMesher([gridsize-1]).build(ot.Interval(tmin, tmax)) 29 30# %% 31vertices = mesh.getVertices() 32 33# %% 34# Creation of the input distribution. 35 36# %% 37distZ0 = ot.Uniform(100.0, 150.0) 38distV0 = ot.Normal(55.0, 10.0) 39distM = ot.Normal(80.0, 8.0) 40distC = ot.Uniform(0.0, 30.0) 41distribution = ot.ComposedDistribution([distZ0, distV0, distM, distC]) 42 43# %% 44dimension = distribution.getDimension() 45dimension 46 47 48# %% 49# Then we define the Python function which computes the altitude at each time value. In order to compute all altitudes with a vectorized evaluation, we first convert the vertices into a Numpy `array` and use the Numpy functions`exp` and `maximum`: this increases the evaluation performance of the script. 50 51# %% 52def AltiFunc(X): 53 g = 9.81 54 z0 = X[0] 55 v0 = X[1] 56 m = X[2] 57 c = X[3] 58 tau = m / c 59 vinf = - m * g / c 60 t = np.array(vertices) 61 z = z0 + vinf * t + tau * (v0 - vinf) * (1 - np.exp(- t / tau)) 62 z = np.maximum(z, 0.) 63 return [[zeta[0]] for zeta in z] 64 65 66# %% 67# In order to create a `Function` from this Python function, we use the `PythonPointToFieldFunction` class. Since the altitude is the only output field, the third argument `outputDimension` is equal to `1`. If we had computed the speed as an extra output field, we would have set `2` instead. 68 69# %% 70outputDimension = 1 71alti = ot.PythonPointToFieldFunction( 72 dimension, mesh, outputDimension, AltiFunc) 73 74# %% 75# Compute a training sample. 76 77# %% 78size = 2000 79ot.RandomGenerator.SetSeed(0) 80inputSample = distribution.getSample(size) 81outputSample = alti(inputSample) 82 83# %% 84# Compute the KL decomposition of the output 85# ------------------------------------------ 86 87# %% 88algo = ot.KarhunenLoeveSVDAlgorithm(outputSample, 1.0e-6) 89algo.run() 90KLResult = algo.getResult() 91scaledModes = KLResult.getScaledModesAsProcessSample() 92 93# %% 94graph = scaledModes.drawMarginal(0) 95graph.setTitle('KL modes') 96graph.setXTitle(r'$t$') 97graph.setYTitle(r'$z$') 98view = viewer.View(graph) 99 100# %% 101# We create the `postProcessingKL` function which takes coefficients of the the K.-L. modes as inputs and returns the trajectories. 102 103# %% 104karhunenLoeveLiftingFunction = ot.KarhunenLoeveLifting(KLResult) 105 106# %% 107# The `project` method computes the projection of the output sample (i.e. the trajectories) onto the K.-L. modes. 108 109# %% 110outputSampleChaos = KLResult.project(outputSample) 111 112# %% 113# We limit the sampling size of the Lilliefors selection in order to reduce the computational burden. 114 115# %% 116ot.ResourceMap.SetAsUnsignedInteger( 117 "FittingTest-LillieforsMaximumSamplingSize", 1) 118 119# %% 120# We create a polynomial chaos metamodel which takes the input sample and returns the K.-L. modes. 121 122# %% 123algo = ot.FunctionalChaosAlgorithm(inputSample, outputSampleChaos) 124algo.run() 125chaosMetamodel = algo.getResult().getMetaModel() 126 127# %% 128# The final metamodel is a composition of the KL lifting function and the polynomial chaos metamodel. In order to combine these two functions, we use the `PointToFieldConnection` class. 129 130# %% 131metaModel = ot.PointToFieldConnection( 132 karhunenLoeveLiftingFunction, chaosMetamodel) 133 134# %% 135# Validate the metamodel 136# ---------------------- 137 138# %% 139# Create a validation sample. 140 141# %% 142size = 10 143validationInputSample = distribution.getSample(size) 144validationOutputSample = alti(validationInputSample) 145 146# %% 147graph = validationOutputSample.drawMarginal(0) 148graph.setColors(['red']) 149graph2 = metaModel(validationInputSample).drawMarginal(0) 150graph2.setColors(['blue']) 151graph.add(graph2) 152graph.setTitle('Model/metamodel comparison') 153graph.setXTitle(r'$t$') 154graph.setYTitle(r'$z$') 155view = viewer.View(graph) 156plt.show() 157 158# %% 159# 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. 160