1"""
2Use case : viscous free fall
3============================
4"""
5from __future__ import print_function
6import openturns as ot
7import numpy as np
8
9
10def AltiFunc(X):
11    g = 9.81
12    z0 = X[0]
13    v0 = X[1]
14    m = X[2]
15    c = X[3]
16    tau = m / c
17    vinf = - m * g / c
18    t = np.array(vertices)
19    z = z0 + vinf * t + tau * (v0 - vinf) * (1 - np.exp(- t / tau))
20    z = np.maximum(z, 0.)
21    return [[zeta[0]] for zeta in z]
22
23
24class ViscousFreeFall():
25    """
26    Data class for the viscous free fall.
27
28
29    Attributes
30    ----------
31
32
33    dim : The dimension of the problem
34          dim=4.
35
36    outputDimension : The output dimension of the problem
37                      outputDimension=1.
38
39    tmin : Constant
40           Minimum time, tmin = 0.0
41
42    tmax : Constant
43           Maximum time, tmax = 12.0
44
45    gridsize : Constant
46               Number of time steps, gridsize = 100.
47
48    mesh : `IntervalMesher`
49
50    vertices : Vertices of the mesh
51
52    distZ0 : `Uniform` distribution of the initial altitude
53             ot.Uniform(100.0, 150.0)
54
55    distV0 : `Normal` distribution of the initial speed
56             ot.Normal(55.0, 10.0)
57
58    distM : `Normal` distribution of the mass
59            ot.Normal(80.0, 8.0)
60
61    distC : `Uniform` distribution of the drag
62            ot.Uniform(0.0, 30.0)
63
64    distribution : `ComposedDistribution`
65                   The joint distribution of the input parameters.
66
67    alti : `PythonPointToFieldFunction`, the exact solution of the fall
68           ot.PythonPointToFieldFunction(dim, mesh, outputDimension, AltiFunc)
69
70
71    Examples
72    --------
73    >>> from openturns.usecases import viscous_free_fall as viscous_free_fall
74    >>> # Load the viscous free fall example
75    >>> vff = viscous_free_fall.ViscousFreeFall()
76    """
77
78    def __init__(self):
79        self.dim = 4  # number of inputs
80        self.outputDimension = 1  # dimension of the output
81
82        self.tmin = 0.0  # Minimum time
83        self.tmax = 12.0  # Maximum time
84        self.gridsize = 100  # Number of time steps
85        self.mesh = ot.IntervalMesher(
86            [self.gridsize-1]).build(ot.Interval(self.tmin, self.tmax))
87        self.vertices = self.mesh.getVertices()
88
89        # Marginals
90        self.distZ0 = ot.Uniform(100.0, 150.0)
91        self.distV0 = ot.Normal(55.0, 10.0)
92        self.distM = ot.Normal(80.0, 8.0)
93        self.distC = ot.Uniform(0.0, 30.0)
94
95        # Joint distribution
96        self.distribution = ot.ComposedDistribution(
97            [self.distZ0, self.distV0, self.distM, self.distC])
98
99        # Exact solution
100        self.alti = ot.PythonPointToFieldFunction(
101            self.dim, self.mesh, self.outputDimension, AltiFunc)
102