(sect:randompressure)=
# Random Pressure Loads

In a previous section we discussed random vibration input (see section
{numref}`%s <randomvib>`). That section addresses a loading where the frequency
content (or power spectral density) of the loading is known for a few points
on the structure. In contrast,
for hypersonic vehicles
a random loading may occur at
thousands of points on the surface. Many aspects of the loading are the same,
but the specification is different, and for performance reasons, the
solutions are performed differently.

The starting point for this analysis requires the following.

1.  A surface sideset where the loading will be applied.

2.  A temporal correlation function to apply on the surface. The temporal correlation
    function is the inverse Fourier transform of a power spectral density (PSD).

3.  A spatial correlation relation. Currently, that relation may only be specified
    as a pair of exponential decay constants.

Details of the problem setup may be found in the User's Manual. This section provides
a simple example of the setup and an informal discussion of the sources of
the data.

## Example Problem Set-up

For our example, we consider a cylinder in a flow field as shown in Figure {numref}`%s <fig:cyl2>`
with the [full input](/tests/Salinas_rtest/training/exampleproblem/transient/load/randompressure/cylinder_random.inp) attached.
The structure is a right circular cylinder of diameter 1
unit, and height 2 units. The flow is directed towards this cylinder in the $X$
direction, and the PSD and corresponding temporal correlation function are
shown in Figure {numref}`%s <fig:randpresspsd>`.

:::{figure-md} fig:cyl2

![](mesh.png)

[Example Random Pressure Geometry](/tests/Salinas_rtest/training/exampleproblem/transient/load/randompressure/cylinder_random.exo)
:::

:::{subfigure} AB
:subcaptions: below
:name: fig:randpresspsd

![](PSD.png)
![](cfun.png)

Example Random Pressure PSD and Correlation Functions.
:::

We are interested in this example, in frequencies up to 500 Hz, so the
cutoff frequency is 500 Hz. There is no point in adding energy above the
desired cutoff frequency -- it only complicates the procedure. [^1] The PSD of the input thus controls much of the solution.

The spatial correlation is often more difficult to obtain. For our example,
we require a decay constant of 2.0 units in the flow direction, and 5.0
perpendicular to the flow. One can think of corresponding decay distances of
0.5 and 0.2 respectively. Thus, down the flow, points more than about 1.5
units away will not be well correlated. [^2]
Perpendicular to the flow, correlations decay even faster.

One rarely has much experimental data about the spatial correlation. Some
information is sometimes garnered from the temporal correlation. For example,
if the correlation function has a characteristic time, $\tau$, one would
expect the spatial correlation length to be of the order of $\delta = v\tau$,
where $v$ is the flow velocity. For a structure in a fluid, the dimensions of
the turbulent layer also provide a bound on the spatial correlation.

## Example: Input Specifications

The physical quantities of the previous section can be interpreted and
expressed as Sierra/SD input as follows.

- The temporal correlation function of Figure {numref}`%s <fig:randpresspsd>` can be
  digitized as a Sierra/SD "time only" function. In input {numref}`%s <fig:randpress:fnctn>` we use a
  triangular pulse for simplicity. The correlation function should be symmetric
  about the origin, and it should have the value of 1.0 at the origin. The
  `correlation_function` is used in the load section, as shown
  in input {numref}`%s <fig:randpress:load>`.

- In the loads section, we also define the following quantities.

  :::{input}
  cutoff_freq = 500
  coordinate = 1 # to set flow direction
  ntimes = 5 # varies from 3 to 20. Too small causes poor replication
             # of the temporal correlation function. Too large results
             # in ill conditioning and singularity.
  :::

Recall that the full correlation matrix is a tensor product of the spatial
correlation with temporal components. The "NTimes" parameter controls the
number of samples in the time domain.

All that remains is setting the spatial correlation decay constants in the
loads section. The full text is shown in input {numref}`%s <fig:randpress:load>`.

:::{input}
correlation_length_z = 0.5
correlation_length_r = 0.2
:::

:::{input}
:caption: Random Pressure Correlation Function
:name: fig:randpress:fnctn

FUNCTION 1
  type linear
    data   -0.001909859319285   0
    data   0   1
    data   0.0019098593192856   0
END
:::

:::{input}
:caption: Random Pressure Load Section. Note that the "flow" direction is the $Z$ coordinate direction of coordinate frame 1.
:name: fig:randpress:load

LOADS
  sideset 1
    randompressure
    cutoff_freq = 500
    delta_t = 0.001
    correlation_function = 1
    ntimes   = 5
    correlation_length_z = 0.5
    correlation_length_r = 0.2
    coordinate = 1
END
BEGIN RECTANGULAR COORDINATE SYSTEM 1
  origin 0 0 0
  z point 1 0 0
  xz point 1 0 1
END
:::

## Example: Verifying the Load

This is a fairly complex input, and it is advisable to verify the generated
loads to ensure consistency. We examine four quantities.

1.  average force on a node.

2.  variance of the force on a node.

3.  temporal force correlation on a single node.

4.  cross correlation of forces between nodes.

All these quantities require output of the total input force, which is
obtained by specifying `force` in the `OUTPUTS` section of the Sierra/SD
input. We will use MATLAB tools to evaluate many of
the results. Data can be read into MATLAB from the Exodus results using
`exo2mat` or using other methods.

### Average Nodal Force

The average nodal force may be determined either by evaluating the MATLAB
results directly, or using the `statistics` output from Sierra/SD. The built
in statistical output is easiest. Supply the `statistics` keyword to the
`OUTPUTS` section, and results will be written to an additional Exodus
file. This has the added benefit that these results may be easily visualized
using Paraview or Ensight. See Figure {numref}`%s <randpress:out>`.

For long time integration, the average value of the nodal force should
approach zero. Shorter time samples will have greater variation. The random
variables depend on `cutoff_freq`. The number of random samples can be
computed as,


$$
N_{samples} = Time_{analysis} \cdot \mbox{cutoff\_freq}
$$

The fractional mean of the force should be within about $3/\sqrt{N_{samples}}$. Or,

$$
Error_{mean} = \left|\frac{F_{mean}}{F_o}\right| < \frac{3}{\sqrt{N_{samples}}}
$$

Here $F_o$ is the force applied for a correlation function of 1. It
involves the scale factors of the function, the sideset distribution
factors and the effective area for each node. [^3]
See the comments section, {numref}`%s <randompress:comments>` for, discussion on the
effective area.

For the example in Figure {numref}`%s <randpress:out>`, mean forces are of the order of 1/1000.
In this example, we took 10,000 time steps, with each of 0.1 ms for a total
time $Time_{analysis}=1\,s$. With $Delta\_T=1/\mbox{cutoff\_freq}=1\,ms$,
the total number of random samples is $N_{samples} = 1000$.

For nodes in the center of the loading area, the effective area is about 0.0098
square units. Because the sideset distribution factors are all one, we have
$F_o=0.0098$. Then,


$$
Error_{mean} = \frac{0.001}{0.0098} = 0.1
$$

which is greater than $\frac{3}{\sqrt{1000}}= 0.095$. A distribution
of the mean is shown in Figure {numref}`%s <fig:randpress:distmean>`.

:::{subfigure} AB
:subcaptions: below
:name: randpress:out

![](mean.png)
![](std.png)

Variation of Mean and Standard Deviation of Force Magnitude on the Surface.
:::

:::{figure-md} fig:randpress:distmean

![](meandist.png)

Distribution of Mean Forces on Surface.
:::

### Variance of Nodal Force

The standard deviation, which is the square root of the variance, is also
available as an output from the analysis, and may be
plotted on the structure using standard visualization tools. See Figure
{numref}`%s <randpress:out>`.

Again, the standard deviation is a statistical quantity, which is only
meaningful for large numbers of samples. In the limit of large $N$, the
standard deviation should approach $F_o$, as defined above, provided that the
correlation function is 1 at time 0.

The plots show a value of $F_{std} \approx 0.0085$ which is under the
expected value of 0.0098. Because the averaging process tends to round out
the correlation function, the measured values of the standard deviation are
typically somewhat less than $F_o$. The autocorrelation function analysis of
the following section should make this more clear.

### Temporal Nodal Force Autocorrelation

The statistics of the loading on a single node should recover the initial input temporal correlation.
Figure {numref}`%s <fig:tstats>` shows the correlation function extracted from the raw time
response data. The correlation function may be computed as,


$$
f_c(n) = \frac{1}{F_o^2} \sum_i w_i w_{i-n}.
$$

Where $w_i$ is the force on a node at time $t_i$.
This data can only be obtained using MATLAB or another external tool, i.e. it
is not available as part of the statistical output. In MATLAB we get this with,
`C = xcorr(f1,f1)`,
where `f1` is the force time history on a node of the
surface. We recover a correlation that is *similar* to the
original triangle correlation in the input. Because of interpolation and
finite sample length, we do not expect the same curves precisely.

The curves of Figure {numref}`%s <fig:tstats>` should be considered "good enough" in
a statistical sense. A temporal interpolation from multiples of ` Delta_T` to the integration time step is being performed, which smooths the
values of the correlation.

:::{figure-md} fig:tstats

![](autocorrelation.png)

Nodal Force Autocorrelation.
:::

### Spatial Cross Correlation

The previous section discussed the autocorrelation function, i.e. the
temporal correlation of signals on the same spatial location. We can also
examine the cross correlation functions. We will only evaluate the functions
at the peak.

This is more difficult. We use the MATLAB "find" method to get the
indices of the nodes with $x=-0.5$, and $y=0$. We loop through these nodes,
and compute the
"xcorr"
function between the node at the center and the
other nodes. The peak value of this solution is then plotted versus the
distance in Figure {numref}`%s <fig:rp_crosscorrelation>`.

:::{figure-md} fig:rp_crosscorrelation

![](rp_crosscorrelation.png)

Nodal Force Spatial Cross Correlation.
:::

There are obvious differences between the measured loads and the target. The
correlations for close distances are lower. This is understood to be
generated by the temporal interpolation function. At large distances, the
cross correlations never go to zero because of the finite length of the
sample.

(randompress:comments)=
## Random Pressure Comments
### Effective Area

Random pressures are computed as force loads using a consistent pressure
calculation. Pressures at the nodes are spread through the element shape
functions to result in nodal forces. For a uniform mesh, this is similar
to lumping the pressures from a fixed area onto the nodes with $F=P\cdot Area$.
In Figure {numref}`%s <fig:effectivearea>` an element based mesh is shown along a
corresponding effective area for the nodes. For a uniform quadrilateral mesh
like the example above, the nodal effective area is the same as the area of an
element face.

:::{figure-md} fig:effectivearea

![](effectivearea.png)

Nodal Effective Area.
:::

### Temporal Interpolation

To improve performance, the random pressure loading procedure computes random
pressures at multiples of "Delta_T" and then interpolates to integration
time steps. A piecewise linear interpolation introduces unacceptable errors;
sinc interpolation is much better.

Interpolation can be avoided by choosing the integration and sampling times to
be equal. In no case should the integration time be larger than the sampling
time.

### Singularities

To compute the proper temporal and spatial correlations for the forces, we
need to perform a Cholesky factorization of the correlation matrix. This
factor will fail if the matrix is singular. Remember that the correlation
matrix that we factor is a tensor product of temporal and spatial
components, $C=C_{spatial} \otimes C_{temporal}$. If either component is
singular, the matrix $C$ is singular. Several common issues can cause
singularity of this matrix.

1.  Taking NTimes too large or too small. For small Delta_T, NTimes must be
    large enough to sample the time correlation function. However, studies show
    that the condition number of the matrix grows exponentially with NTimes. The
    target value is 5. Values above 20 are not recommended; $C_{temporal}$ is
    numerically singular.

2.  Spatial degeneracy, leading to $C_{spatial}=0$. We have only one means of
    entering the spatial correlation parameters, viz. the `correlation_length`
    variables pair. If either of these quantities are so
    large that $\delta / \mbox{correlation\_length}$ is very close to zero (with
    $\delta$ representing the distance from one node to another on the mesh), then
    the spatial portion of the matrix becomes singular. Effectively, these
    locations are no longer independent, but must apply the same load vector.

3.  Using a
    Delta_T = 1/cutoff_freq and the default sinc function for
    a correlation function may generate a $C_{temporal}$ singularity.
    [^4]
    This is
    because we are now evaluating the correlation function at multiples of $\pi$,
    where it is always zero.

### Time Step

The integration time step specified in the `SOLUTION` section must *always*
be less than or equal to Delta_T.

### Sinc Function

The sinc function defined as $\sin{(x)}/x$ is important in at least two
places in the code. First, it is the only function available for the temporal
interpolation function. Second, by default, we use the sinc function
as the correlation function. In most cases, this use of the function should
probably be replaced by another function. We use it as the default because it
represents the Fourier transform of a flat PSD, which is the simplest
loading.

## Memory, Performance, Parallel and Anything Else of Interest

The matrices generated for these operations are all square and dense. The
matrix order is $d=n_{spatial}\cdot n_{temporal}$. Here
$n_{spatial}$ is the number of points in the surface and $n_{temporal}$
is the value of the `NTimes` parameter.
Because memory requirements grow as the square of these variables,
it is important to manage these carefully. Practically, models up to $d=10^5$
are possible in parallel, but they take a lot of time.

The operation count for Cholesky factorization of a dense matrix is of order
$d^3$. Thus, the computational cost increases much faster than model size.
**The parallel solutions of the Cholesky system are not scalable**. In a
scalable problem, doubling the size of the problem, and also doubling the
number of processors should not change the solution time. Although the sparse linear
solvers for FE solution are scalable, the Cholesky factorization required
to compute random pressure loads is not scalable.

The dense Cholesky factorization from the ScaLAPACK library is used.
The parallel decomposition for this solve is completely different from the FEM
decomposition, and is computed internally without user intervention. The user
input for the parallel solution is exactly the same as the serial input.
However, at this time, parallel solutions are limited to platforms built
under the Intel compiler with MKL libraries. The solution will fail on other
platforms.

[^1]: Although the physics has energy above $500$ Hz, cutting off the PSD at 500 Hz. is
    required because a
    higher cutoff frequency narrows the correlation function with no added
    accuracy.

[^2]: correlation=$\exp{(-3)}=0.0498$

[^3]: A simple way to estimate $F_o$ is to run a very short transient analysis
    after having converted the random pressure load to a constant unit pressure.

[^4]: In this example, we intentionally use the triangular function both for
    simplicity, and to avoid this problem.
