Using the pyYeti CLA tools
==========================

Introduction
------------

Many thanks to **Thomas Weber** for his initial draft of this tutorial!

This and other jupyter notebooks are available here:
https://github.com/twmacro/pyyeti/tree/master/docs/tutorials.

This tutorial will use a simple model to guide through the CLA process
using the pyYeti tools. The model is a simple space station model
contained in outboard.blk and inboard.blk. The outboard model will be
used to create an external superelement, as if this model was delivered
as a Craig-Bampton model by an independent group. The inboard model is
the core of the residual structure. In the rocket launching business,
the outboard model is analogous to a spacecraft model or an upstream
component on the rocket, and the inboard model is analogous to the
non-upstream part of the launch vehicle model.

Outboard:

.. figure:: outboard.png
   :alt: outboard.png

   outboard.png

Inboard:

.. figure:: inboard.png
   :alt: inboard.png

   inboard.png

Outline
-------

The general outline for this tutorial is:

-  Outboard superelement creation
-  Prepare outboard model for CLA
-  Run modes on assembled model
-  Simulate loading events
-  Summarize results
-  Compare results

Special note about this tutorial
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In order to speed up the automatic documentation creation step on “Read
the Docs”, a number of lines below are commented out. These all start
with ``# SPEED:``. These lines all produce output that is not contained
within this notebook and would therefore only cost execution time while
creating the documentation. However, if you decide to run this tutorial
for yourself, it is recommended that these lines be uncommented so you
can take a look at all the output. Note that some screenshots of the
output is included below.

Nastran runs
------------

The outboard model will be used to create an external superelement to
simulate a delivery from an independent team. The residual structure run
will bring in the outboard model as an external superelement as if we
did not have the bulk data. The inboard bulk data is included directly.
Both .dat files are shown below for convenience; they exist in the
``srcdir`` directory that is shown below. These files were run in
Siemens Nastran version 2020.

.. code:: ipython3

    import inspect
    import os
    from pathlib import Path
    import pyyeti
    
    srcdir = Path(inspect.getfile(pyyeti)).parent / "tests" / "cla_test_data_2020"
    print(f"{srcdir = }")
    
    wrkdir = Path.cwd()
    print(f"{wrkdir = }")


.. code-block:: none

    srcdir = PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/pyyeti/envs/stable/lib/python3.12/site-packages/pyyeti/tests/cla_test_data_2020')
    wrkdir = PosixPath('/home/docs/checkouts/readthedocs.org/user_builds/pyyeti/checkouts/stable/docs/tutorials')


outboard.dat
~~~~~~~~~~~~

The following Nastran input file creates the “outboard” external
superelement (SE 101). This results in 3 files for use in the residual
run: outboard.op4, outboard.asm, and outboard.pch.

In cases where a Craig-Bampton model is provided but the .op4, .asm, and
.pch files are not provided, the routine
`pyyeti.nastran.bulk.wtextseout <../modules/nastran/generated/pyyeti.nastran.bulk.wtextseout.html#pyyeti.nastran.bulk.wtextseout>`__
can be used to create these three files. This is typically done in the
“prepare_4_cla.py” script shown below.

File “outboard.dat”:

.. code:: text

       INIT MASTER(S)  $ delete .MASTER and .DBALL files on exit
       ASSIGN OUTPUT4='outboard.op4' UNIT=101,DELETE
       SOL 103
       CEND
       
       TITLE = Outboard
       ECHO = SORT
       WEIGHTCHECK(SET=ALL) = YES
       GROUNDCHECK(SET=ALL,DATAREC=YES) = YES
       METHOD=1
       SET 1 = 1 THRU 48
       DISPLACEMENT(PLOT) = 1
       FORCE(PLOT) = all
       EXTSEOUT(ASMBULK,EXTBULK,EXTID=101,MATOP4=101)
       
       BEGIN BULK
       EIGRL          1            50.0
       SPOINT   1995001    THRU 1995022
       QSET1            1995001    THRU 1995022
       BSET1,123456,3,11,19,27
       include 'outboard.blk'
       ENDDATA

residual.dat
~~~~~~~~~~~~

Create the residual structure. This file uses inboard.blk and brings in
SE 101 via the outboard.op4, outboard.asm and outboard.pch files. This
run creates the nas2cam.op2 and nas2cam.op4 files that will be used for
analysis in Python.

File “residual.dat”:

.. code:: text

       NASTRAN SYSTEM(402) = 0  $ AUTOMATICALLY DELETE DUPLICATE CARDS
       NASTRAN NLINES = 10000
       ASSIGN INPUTT4='outboard.op4',UNIT=101
       INIT MASTER(S)  $ delete .MASTER and .DBALL files on exit
       $ NAS2CAM op2/op4 files:
       assign output2 = 'nas2cam.op2', status=new, unit=29,delete $
       assign output4 = 'nas2cam.op4', status=new, unit=30,
            form=unformatted,delete $
       
       DIAG    8,47
       SOL 111
       echooff
       include '../nas2cam/nas2cam_111.v9'
       include '../nas2cam/nas2cam_subdmap_2023.v9'
       
       $ bug fix (?) alter for v2020:
       COMPILE PHASE0
       $ - delete line that prevents RVDOF from being used for resvecs
       $ - must also include RESVEC(DYNRSP)=YES if need damping on these resvecs
       ALTER 'IF ( NOT(RESVEC0) )'(2),'IF ( NOT(RESVEC0) )'(2) $ DELETE
       $ Note: an alternative to RVDOF is to use the older PARAM,RESVEC,YES and
       $ USET,U6 approach. That works without this little alter.
       ENDALTER $
       
       echoon
       CEND
       
       TITLE = System Modes
       ECHO = Sort
       METHOD = 1
       FREQ = 1
       DLOAD = 1
       DISPLACEMENT(PLOT) = ALL
       FORCE(PLOT) = ALL
       WEIGHTCHECK(SET=ALL) = YES
       GROUNDCHECK(SET=ALL,DATAREC=YES) = YES
       RESVEC(NOAPPL,RVDOF,NORVEL,NOINRL,NODAMP,NODYNRSP)=YES
       
       $-------------------------------------------------------------------
       $ nas2cam params:
       PARAM,PRFMODES,1
       $
       $ TO GENERATE GRAVITY FORCE, SET GRAVDIR EQUAL TO GRAVITY DIRECTION
       $  AND SET THE GRAVITY FIELD:
       $
       PARAM,GRAVDIR,3
       PARAM,GRAVFELD,9806.65  $ mm/sec**2
       $-------------------------------------------------------------------
       
       SUBCASE 1
           LABEL = Modes run with BHH matrix
           ANALYSIS = MODES
       
       BEGIN BULK
       $-------------------------------------------------------------------
       $ NAS2CAM input:
       PARAM,DBDICT,0
       DTI,TMAA,1,101,0
       DTI,TKAA,1,101,0
       DTI,TGM,1,0
       DTI,TPHG,1,0
       DTI,TPHA,1,0
       DTI,TBHH,1,0
       $-------------------------------------------------------------------
       PARAM,POST,-1
       PARAM,SESDAMP,YES
       EIGRL          1           150.0
       $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
       FREQ           1      2.
       RLOAD2         1       1                       1
       DAREA          1      11       1    1.0
       TABLED1        1
                   0.01     1.0   150.0    1.0     ENDT
       $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
       $ for residual flexibility vectors:
       RVDOF1,123,8,22,24
       INCLUDE 'outboard.asm'
       include 'inboard.blk'
       INCLUDE 'outboard.pch'
       ENDDATA

Prepare Outboard for CLA
------------------------

The prepare_4_cla.py file below prepares the outboard model for CLA. The
primary goal is to setup data recovery.

This example also runs
`pyyeti.cb.cbcheck <../modules/generated/pyyeti.cb.cbcheck.html#pyyeti.cb.cbcheck>`__,
though this could be done in a separate model checkout run. This is not
discussed further here, but is the subject of a different tutorial:
`Using cb.cbcheck to check mass and
stiffness <../tutorials/cbcheck.html>`__.

A special data recovery category: “cglf”
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The CG load factor category is special in that row 11 (at Python index
10) is a time-consistent RSS (root-sum-square) of rows 2 and 3. Rows 2
and 3 are the two 90 degrees apart shear-based lateral CG load factors
in the SC coordinate system, so row 11 becomes the RSS shear-based
lateral load factor. Similarly, row 12 is the RSS of 4 and 5, 13 is the
RSS of 7 and 8, and 14 is the RSS of 9 and 10. Together, these rows
cover both the shear-based and moment-based lateral CG load factors in
the SC and LV coordinate systems. The corresponding RSS values for both
coordinate systems should match each other. For PSD (power spectral
density) analyses, an analogous RSS is handled with the help of
`pyyeti.cla.PSD_consistent_rss <../modules/generated/pyyeti.cla.PSD_consistent_rss.html#pyyeti.cla.PSD_consistent_rss>`__.
The data recovery matrix for this category is created by
`pyyeti.cb.mk_net_drms <../modules/generated/pyyeti.cb.mk_net_drms.html#pyyeti.cb.mk_net_drms>`__
and the last four rows (11-14) have only zeros.

To handle the RSS’ing during data recovery, the file “dr_file.py” is
created for the “cglf” data recovery category with the routines “cglf”
and “cglf_psd”. All the other data recovery categories are simple and
handled directly in the “prepare_4_cla.py” file. See also
`pyyeti.cla.DR_Def.add <../modules/edited/pyyeti.cla.DR_Def.add.html#pyyeti.cla.DR_Def.add>`__
for more details on creating data recovery categories.

For reference, the contents of “dr_file.py” are shown here:

::

       import numpy as np
       from pyyeti import cla
       
       
       def get_xyr():
           # return the xr, yr, and rr indexes for the "cglf" data recovery
           xr = np.array([1, 3, 6, 8])  # 'x' row(s)
           yr = xr + 1  # 'y' row(s)
           rr = np.arange(4) + 10  # rss  rows
           return xr, yr, rr
       
       
       def cglf(sol, nas, Vars, se):
           resp = Vars[se]["cglf"] @ sol.a
           xr, yr, rr = get_xyr()
           resp[rr] = np.sqrt(resp[xr] ** 2 + resp[yr] ** 2)
           return resp
       
       
       def cglf_psd(sol, nas, Vars, se, freq, forcepsd, drmres, case, i):
           resp = Vars[se]["cglf"] @ sol.a
           cla.PSD_consistent_rss(resp, *get_xyr(), freq, forcepsd, drmres, case, i)

Miscellaneous notes regarding the “prepare_4_cla.py” file:

-  The dynamic uncertainty factor is set to 1.25
-  The SRS (shock response spectrum) is computed for some rows of some
   categories:

   -  Q = 10 and 33
   -  Frequency range from 0.1 to 50.0 Hz with step size 0.1 Hz
   -  The SRS calculations are performed by
      `pyyeti.srs.srs <../modules/generated/pyyeti.srs.srs.html#pyyeti.srs.srs>`__,
      `pyyeti.srs.srs_frf <../modules/generated/pyyeti.srs.srs_frf.html#pyyeti.srs.srs_frf>`__,
      and
      `pyyeti.srs.vrs <../modules/generated/pyyeti.srs.vrs.html#pyyeti.srs.vrs>`__
      for time-domain, frequency-domain, and PSD analyses, respectively.

-  All categories are added by using locally defined functions (all
   named "_") and the
   `pyyeti.cla.DR_Def.addcat <../modules/edited/pyyeti.cla.DR_Def.addcat.html#pyyeti.cla.DR_Def.addcat>`__
   decorator.
-  An elastic mode only “net_ifatm” category is added via
   `pyyeti.cla.DR_Def.add_0rb <../modules/edited/pyyeti.cla.DR_Def.add_0rb.html#pyyeti.cla.DR_Def.add_0rb>`__.
   The name of the category will be “net_ifatm_0rb” (for zero
   rigid-body).
-  A summary of all categories is printed to an excel file
   “dr_summary.xlsx”. This is meant for visual checking and is shown
   below for reference.
-  Finally, the critical data is saved to “cla_params.pgz” via
   `pyyeti.ytools.save <../modules/generated/pyyeti.ytools.save.html#pyyeti.ytools.save>`__.
   (Both
   `pyyeti.ytools.save <../modules/generated/pyyeti.ytools.save.html#pyyeti.ytools.save>`__
   and
   `pyyeti.ytools.load <../modules/generated/pyyeti.ytools.load.html#pyyeti.ytools.load>`__
   are imported into the “cla” module for convenience.)

.. code:: ipython3

    # prepare_4_cla.py
    
    from pathlib import Path
    import numpy as np
    import re
    from pyyeti import cla, nastran, cb
    from pyyeti.nastran import op4
    
    
    def getlabels(lbl, id_dof):
        return ["{} {:4d}-{:1d}".format(lbl, g, i) for g, i in id_dof]
    
    
    if __name__ == "__main__":
        se = 101
        uset, coords, bset = nastran.asm2uset(srcdir / "outboard.asm")
        bset = bset.nonzero()[0]
        dct = op4.read(srcdir / "outboard.op4")
        maa = dct["mxx"]
        kaa = dct["kxx"]
        atm = dct["mug1"]
        ltm = dct["mef1"]
        pch = srcdir / "outboard.pch"
    
        atm_labels = getlabels("Grid", nastran.rddtipch(pch, "tug1"))
        ltm_labels = getlabels("CBAR", nastran.rddtipch(pch, "tef1"))
        iflabels = getlabels("Grid", uset.index[bset])
    
        # setup CLA parameters:
        mission = "Micro Space Station"
    
        ref = [600.0, 150.0, 150.0]
        g = 9806.65
        net = cb.mk_net_drms(maa, kaa, bset, uset=uset, ref=ref, g=g)
    
        # run cbcheck:
        # SPEED: chk = cb.cbcheck(
        # SPEED:     "outboard_cbcheck.out", maa, kaa, bset, bref=np.arange(6), uset=uset, uref=ref
        # SPEED: )
    
        # define some defaults for data recovery:
        defaults = dict(
            se=se,
            uf_reds=(1, 1, 1.25, 1),
            srsfrq=np.arange(0.1, 50.1, 0.1),
            srsQs=(10, 33),
            drfile=srcdir / "dr_file.py"
        )
    
        drdefs = cla.DR_Def(defaults)
    
        @cla.DR_Def.addcat
        def _():
            name = "scatm"
            desc = "Outboard Internal Accelerations"
            units = "mm/sec^2, rad/sec^2"
            labels = atm_labels
            drms = {"atm": atm}
            drfunc = "Vars[se]['atm'] @ sol.a"
            # want translation histories and srs curves for nodes 35 & 36:
            prog = re.compile(" 3[56]-[123]")
            i = [i for i, s in enumerate(atm_labels) if prog.search(s)]
            histpv = np.zeros(len(labels), bool)
            histpv[i] = True
            srspv = histpv
            srsopts = dict(eqsine=1, ic="steady")
            drdefs.add(**locals())
    
        @cla.DR_Def.addcat
        def _():
            name = "scltm"
            desc = "Outboard Internal Loads"
            units = "mN, mN-mm"
            labels = ltm_labels
            drms = {"ltm": ltm}
            drfunc = "Vars[se]['ltm'] @ sol.d"
            drdefs.add(**locals())
    
        @cla.DR_Def.addcat
        def _():
            name = "ifatm"
            desc = "S/C Interface Accelerations"
            units = "mm/sec^2, rad/sec^2"
            labels = iflabels
            ifatm = np.eye(bset.shape[0], maa.shape[0])
            drms = {"ifatm": ifatm}
            drfunc = "Vars[se]['ifatm'] @ sol.a"
            srsopts = dict(eqsine=1, ic="steady")
            histpv = "all"
            drdefs.add(**locals())
    
        @cla.DR_Def.addcat
        def _():
            name = "ifltm"
            desc = "I/F Loads"
            units = "mN, mN-mm"
            labels = iflabels
            ifltma = maa[bset]
            ifltmd = kaa[bset][:, bset]
            drms = {"ifltma": ifltma, "ifltmd": ifltmd}
            drfunc = "Vars[se]['ifltma'] @ sol.a + Vars[se]['ifltmd'] @ sol.d"
            drdefs.add(**locals())
    
        @cla.DR_Def.addcat
        def _():
            name = "cglf"
            desc = "S/C CG Load Factors"
            units = "G"
            labels = net.cglf_labels
            drms = {"cglf": net.cglfa}
            histpv = slice(5)
            drdefs.add(**locals())
    
        @cla.DR_Def.addcat
        def _():
            name = "net_ifatm"
            desc = "NET S/C Interface Accelerations"
            units = "g, rad/sec^2"
            labels = net.ifatm_labels
            drms = {"net_ifatm": net.ifatm}
            drfunc = "Vars[se]['net_ifatm'] @ sol.a"
            srsopts = dict(eqsine=1, ic="steady")
            histpv = "all"
            drdefs.add(**locals())
    
        @cla.DR_Def.addcat
        def _():
            name = "net_ifltm"
            desc = "NET I/F Loads"
            units = "mN, mN-mm"
            labels = net.ifltm_labels
            drms = {"net_ifltm": net.ifltma}
            drfunc = "Vars[se]['net_ifltm'] @ sol.a"
            srsopts = dict(eqsine=0, ic="steady")
            histpv = "all"
            drdefs.add(**locals())
    
        # add a 0rb version of the NET ifatm:
        drdefs.add_0rb("net_ifatm")
    
        # make excel summary file for visual checking:
        # SPEED: drdefs.excel_summary("dr_summary.xlsx")
    
        # save data to gzipped pickle file:
        sc = dict(mission=mission, drdefs=drdefs)
        cla.save("cla_params.pgz", sc)
    
        print("Done!")


.. code-block:: none

    Done!


For reference, here is a screenshot of “dr_summary.xlsx”:

.. figure:: clapng/dr_summary.png
   :alt: dr_summary.png

   dr_summary.png

Run Events
----------

There are 4 events run in this CLA: transfer orbit engine start (TOES),
oil and water mixing experiment (OWLab), transfer orbit burn (TOBurn),
and transfer orbit engine cutoff (TOECO). After these events are run,
the results are summarized and compared.

First, for this tutorial, we’ll define a small convience function to
ease running each event in it’s own subdirectory:

.. code:: ipython3

    def set_event_dir(event):
        os.chdir(wrkdir)
        event_dir = Path(event)
        event_dir.mkdir(exist_ok=True)
        os.chdir(event_dir)

TOES
~~~~

Outline of Transfer Orbit Engine Start run script:

-  Load data recovery data
-  Load Nastran data
-  Form ULVS for the outboard model (the SC)

   -  ULVS is a row partition of the system modes to the superelement
      external DOF (typically the b-set and q-set DOF)

-  Prepare spacecraft data recovery matrices
-  Initialize results (ext, mnc, mxc for all drms)
-  Set rfmodes. This typically defines which modes are
   residual-flexibility modes, but really defines which modes are to be
   treated statically. See
   `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__
   for more information.
-  Setup modal mass, damping, and stiffness

   -  Damping is diagonal, 2% modal damping

-  Load in forcing functions
-  Form force transform
-  Do ODE pre-calcs
-  Loop over all cases, solving the ODEs and performing data recovery
-  Compute the P99/90 statistical maximum and minimum values for each
   item
-  Save results and create tables and plots

   -  Some screen shots of the tables and shown below for reference

.. code:: ipython3

    set_event_dir("toes")

.. code:: ipython3

    # Simulate event and recover responses
    import numpy as np
    from scipy.io import matlab
    from pyyeti import stats, ode, cla
    from pyyeti.nastran import n2p, op2
    
    
    if __name__ == "__main__":
        # event name:
        event = "TOES"
    
        # load data recovery data:
        sc = cla.load(wrkdir / "cla_params.pgz")
        cla.PrintCLAInfo(sc["mission"], event)
    
        # load nastran data:
        nas = op2.rdnas2cam(srcdir / "nas2cam")
    
        # form ulvs for some SEs:
        SC = 101
        n2p.addulvs(nas, SC)
    
        # prepare spacecraft data recovery matrices
        DR = cla.DR_Event()
        DR.add(nas, sc["drdefs"])
    
        # initialize results (ext, mnc, mxc for all drms)
        results = DR.prepare_results(sc["mission"], event)
    
        # set rfmodes:
        rfmodes = nas["rfmodes"][0]
    
        # setup modal mass, damping and stiffness
        m = None  # None means identity
        k = nas["lambda"][0]
        assert nas["nrb"] == 6
        k[: nas["nrb"]] = 0.0
        b = 2 * 0.02 * np.sqrt(k)
        mbk = (m, b, k)
    
        # load in forcing functions:
        toes = matlab.loadmat(srcdir / "toes_ffns.mat", squeeze_me=True, struct_as_record=False)
    
        # form force transform:
        T = n2p.formdrm(nas, 0, [[8, 12], [24, 13]])[0].T
    
        # do pre-calcs and loop over all cases:
        ts = ode.SolveUnc(*mbk, 1 / toes["sr"], rf=rfmodes)
        LC = toes["ffns"].shape[0]
        t = toes["t"]
        for j, force in enumerate(toes["ffns"]):
            print("Running {} case {}".format(event, j + 1))
            genforce = T @ ([[1], [0.1], [1], [0.1]] * force[None, :])
            # solve equations of motion
            sol = ts.tsolve(genforce, static_ic=1)
            sol.t = t
            sol = DR.apply_uf(sol, *mbk, nas["nrb"], rfmodes)
            caseid = "{} {:2d}".format(event, j + 1)
            # perform data recovery:
            results.time_data_recovery(sol, nas["nrb"], caseid, DR, LC, j)
    
        # compute P99/90 statistical extreme values:
        results.calc_stat_ext(stats.ksingle(0.99, 0.90, LC))
    
        # save results:
        cla.save("results.pgz", results)
    
        # make some srs plots and tab files:
        # SPEED: results.rpttab()
        # SPEED: results.srs_plots()
        # SPEED: results.resp_plots()
    
        print("Done!")


.. code-block:: none

    Mission:  Micro Space Station
    Event:    TOES
    Running TOES case 1
    Running TOES case 2
    Running TOES case 3
    Running TOES case 4
    Running TOES case 5
    Running TOES case 6
    Running TOES case 7
    Done!


For reference, here are some screenshots created by the above run:

The
`pyyeti.cla.DR_Results.rpttab <../modules/edited/pyyeti.cla.DR_Results.rpttab.html#pyyeti.cla.DR_Results.rpttab>`__
routine writes tables of results. Here is a the “net_ifatm.tab” file
(with some lines deleted for brevity). The extrema count table at the
bottom shows that case 5 drove most of the “net_ifatm” extreme values.

.. figure:: clapng/toes_net_ifatm_tab.png
   :alt: toes_net_ifatm_tab.png

   toes_net_ifatm_tab.png

The
`pyyeti.cla.DR_Results.srs_plots <../modules/edited/pyyeti.cla.DR_Results.srs_plots.html#pyyeti.cla.DR_Results.srs_plots>`__
routine plots all requested SRS curves to a file or files. Here is the
first page of the “TOES_srs.pdf” file:

|image0|

.. |image0| image:: clapng/toes_srs_pg1.png

The
`pyyeti.cla.DR_Results.resp_plots <../modules/edited/pyyeti.cla.DR_Results.resp_plots.html#pyyeti.cla.DR_Results.resp_plots>`__
routine plots all requested time-history curves to a file or files. Here
is the first page of the “TOES_hist.pdf” file:

.. figure:: clapng/toes_hist_pg1.png
   :alt: toes_hist_pg1.png

   toes_hist_pg1.png

OWLab
~~~~~

Outline for oil and water mixing experiment:

-  Load data recovery data
-  Load Nastran data
-  Form ULVS for the outboard model (the SC)
-  Prepare spacecraft data recovery matrices
-  Initialize Results
-  Set rfmodes
-  Setup modal mass, damping, and stiffness

   -  Damping is diagonal, 2% modal damping

-  Form force transform
-  Define the PSD forces
-  Calculate the PSD responses
-  Save results and make plots

Note:
^^^^^

We’ll see two warnings from this run, one about the frequency step being
too large for accuracy, and another about a division by zero.

The first warning happens during the calculation of the SRS curves
within the routine
`pyyeti.srs.vrs <../modules/generated/pyyeti.srs.vrs.html#pyyeti.srs.vrs>`__.
For the integration frequency vector, this routine merges the
frequencies from the forcing function, which range from 25 to 45 Hz by
0.5 Hz, with the frequencies at which to compute the SRS, which range
from 0.1 to 50 Hz by 0.1 Hz. At the lowest frequencies then, the
delta-frequency is 0.1 Hz which is larger than 0.1 / Q, so we get the
warning. We could refine the SRS frequency vector in the prepare_4_cla
step above to get rid of this warning. However, as noted in
`pyyeti.srs.vrs <../modules/generated/pyyeti.srs.vrs.html#pyyeti.srs.vrs>`__,
the resulting SRS should be conservative and for this tutorial, this is
acceptable. Additionally, in this case, we probably only care about the
SRS in the 25 to 45 Hz range, which should be accurate.

The second warning happens during the calculation of the apparent
frequency inside the
`pyyeti.cla.DR_Results.psd_data_recovery <../modules/edited/pyyeti.cla.DR_Results.psd_data_recovery.html#pyyeti.cla.DR_Results.psd_data_recovery>`__
routine. The “scltm” has 4 zero rows and each of them will cause a
divide-by-zero. For those rows, the “x” coordinate (which is normally
the apparent frequency for PSDs) is set to NaN, which is perfectly fine.

.. code:: ipython3

    set_event_dir("owlab")

.. code:: ipython3

    # Simulate event and recover responses
    import numpy as np
    import scipy.interpolate as interp
    
    from pyyeti import ode, cla
    from pyyeti.nastran import n2p, op2
    
    
    if __name__ == "__main__":
        # event name:
        event = "OWLab"
    
        # load data recovery data:
        sc = cla.load(wrkdir / "cla_params.pgz")
        cla.PrintCLAInfo(sc["mission"], event)
    
        # load nastran data:
        nas = op2.rdnas2cam(srcdir / "nas2cam")
    
        # form ulvs for some SEs:
        SC = 101
        n2p.addulvs(nas, SC)
    
        # prepare spacecraft data recovery matrices
        DR = cla.DR_Event()
        DR.add(nas, sc["drdefs"])
    
        # initialize results (ext, mnc, mxc for all drms)
        results = DR.prepare_results(sc["mission"], event)
    
        # set rfmodes:
        rfmodes = nas["rfmodes"][0]
    
        # setup modal mass, damping and stiffness
        m = None  # None means identity
        k = nas["lambda"][0]
        assert nas["nrb"] == 6
        k[: nas["nrb"]] = 0.0
        b = 2 * 0.02 * np.sqrt(k)
        mbk = (m, b, k)
    
        # form force transform:
        T = n2p.formdrm(nas, 0, [[22, 123]])[0].T
    
        # random part:
        freq = cla.freq3_augment(np.arange(25.0, 45.1, 0.5), nas["lambda"][0])
        rnd = [
            np.array(
                [
                    # freq     x      y      z
                    [1.0, 90.0, 110.0, 110.0],
                    [30.0, 90.0, 110.0, 110.0],
                    [31.0, 200.0, 400.0, 400.0],
                    [40.0, 200.0, 400.0, 400.0],
                    [41.0, 90.0, 110.0, 110.0],
                    [50.0, 90.0, 110.0, 110.0],
                ]
            ),
            np.array(
                [
                    # freq     x      y      z
                    [1.0, 90.0, 110.0, 110.0],
                    [20.0, 90.0, 110.0, 110.0],
                    [21.0, 200.0, 400.0, 400.0],
                    [30.0, 200.0, 400.0, 400.0],
                    [31.0, 90.0, 110.0, 110.0],
                    [50.0, 90.0, 110.0, 110.0],
                ]
            ),
        ]
    
        fs = ode.SolveUnc(*mbk, rf=rfmodes)
        for j, ff in enumerate(rnd):
            caseid = "{} {:2d}".format(event, j + 1)
            print("Running {} case {}".format(event, j + 1))
            F = interp.interp1d(ff[:, 0], ff[:, 1:].T, axis=1, fill_value=0.0)(freq)
            results.solvepsd(nas, caseid, DR, fs, F, T, freq)
            results.psd_data_recovery(caseid, DR, len(rnd), j)
    
        # save results:
        cla.save("results.pgz", results)
    
        # make some srs plots and tab files:
        # SPEED: results.rpttab()
        # SPEED: results.srs_plots(Q=10, direc="srs_cases", showall=True, plot="semilogy")
        # SPEED: results.resp_plots()
    
        print("Done!")


.. code-block:: none

    Mission:  Micro Space Station
    Event:    OWLab
    Running OWLab case 1


.. code-block:: none

    /home/docs/checkouts/readthedocs.org/user_builds/pyyeti/envs/stable/lib/python3.12/site-packages/pyyeti/srs.py:1353: RuntimeWarning: Integration frequency vector may produce inaccurate results; refine the step. See the documentation for more information.
      warn(
    /home/docs/checkouts/readthedocs.org/user_builds/pyyeti/envs/stable/lib/python3.12/site-packages/pyyeti/cla/dr_results.py:1416: RuntimeWarning: invalid value encountered in divide
      pk_freq = vrms / rms


.. code-block:: none

    Running OWLab case 2
    Done!


TOBURN
~~~~~~

TOBurn is unique among the events analyzed in this tutorial because it
uses a combination equation. There are two components: a “noise”
component (solved in the PSD domain), and a steady-state thrust
component (solved in the time-domain). The combination equation is
simply the sum of these two components, noting that the noise component
can be positive or negative.

Outline for Transfer Orbit Burn:

-  Define a function to combine the steady state burn with the noise of
   the burn
-  Load data recovery data
-  Load Nastran data
-  Form ULVS for the outboard model (the SC)
-  Prepare spacecraft data recover matrices
-  Initialize results
-  Set rfmodes
-  Setup modal mass, damping, and stiffness

   -  Damping is diagonal, 2% modal damping

-  Form force transform
-  Calculate steady-state part
-  Calculate random part
-  Calculate combined results
-  Save results and Make plots

.. code:: ipython3

    set_event_dir("toburn")

.. code:: ipython3

    # Simulate event and recover responses
    import numpy as np
    import scipy.interpolate as interp
    from pyyeti import ode, cla
    from pyyeti.nastran import n2p, op2
    
    
    def toburn_combine(event, results):
        """
        TOBurn combination equation
    
        Parameters
        ----------
        event: string
            Name to use for combined results; eg "TOBurn" (stored in, for
            example, ``results['combined']['SC_atm'].event``)
        results : instance of :class:`cla.DR_Results`
            Contains 'ss' and 'noise' instants of :class:`cla.DR_Results`.
            For example, if there is an "atm" category::
    
               results['ss']['atm']
               results['noise']['atm']
    
        Returns
        -------
        None
    
        Notes
        -----
        Adds the 'combined' instance of :class:`cla.DR_Results` to
        `results`.
    
        The combination equation is::
    
            mx = ss + noise
            mn = ss - noise
        """
        # use form_extreme to help form the combined results:
        results.form_extreme(ext_name=event)
        results["combined"] = results["extreme"]
        del results["extreme"]
    
        # now, just fix the "ext" members:
        for cat, sns in results["combined"].items():
            sns.domain = "combination"
            noise = results["noise"][cat]
            ss = results["ss"][cat]
            term = abs(noise.ext).max(axis=1)
            sns.ext[:, 0] = ss.ext[:, 0] + term
            sns.ext[:, 1] = ss.ext[:, 1] - term
            sns.exttime = None
            sns.maxcase = ["Combination"] * sns.ext.shape[0]
            sns.mincase = sns.maxcase
            # srs:
            if getattr(sns, "srs", None):
                _srs = sns.srs
                for Q in _srs.ext:
                    _srs.ext[Q][:] = ss.srs.ext[Q] + noise.srs.ext[Q]
    
    
    if __name__ == "__main__":
        # event name:
        event = "TOBurn"
    
        # load data recovery data:
        sc = cla.load(wrkdir / "cla_params.pgz")
        cla.PrintCLAInfo(sc["mission"], event)
    
        # load nastran data:
        nas = op2.rdnas2cam(srcdir / "nas2cam")
    
        # form ulvs for some SEs:
        SC = 101
        n2p.addulvs(nas, SC)
    
        # prepare spacecraft data recovery matrices
        DR = cla.DR_Event()
        DR.add(nas, sc["drdefs"])
    
        # initialize results (ext, mnc, mxc for all drms)
        results = cla.DR_Results()
        results["ss"] = DR.prepare_results(sc["mission"], event)
        results["noise"] = DR.prepare_results(sc["mission"], event)
    
        # set rfmodes:
        rfmodes = nas["rfmodes"][0]
    
        # setup modal mass, damping and stiffness
        m = None  # None means identity
        k = nas["lambda"][0]
        assert nas["nrb"] == 6
        k[: nas["nrb"]] = 0.0
        b = 2 * 0.02 * np.sqrt(k)
        mbk = (m, b, k)
    
        # form force transform:
        T = n2p.formdrm(nas, 0, [[8, 12], [24, 13]])[0].T
    
        # steady state part:
        case = "ss"
        ts = ode.SolveUnc(*mbk, rf=rfmodes)
        genforce = T @ [[7000.0], [0.0], [7000.0], [0.0]]
        sol = ts.tsolve(genforce, static_ic=1)
        sol = DR.apply_uf(sol, *mbk, nas["nrb"], rfmodes)
        results[case].time_data_recovery(sol, nas["nrb"], case, DR, 1, 0)
    
        # random part:
        case = "noise"
        freq = cla.freq3_augment(np.arange(5.0, 35.1, 0.5), nas["lambda"][0])
        F = interp.interp1d(
            [1.0, 50.0],
            [[300.0, 300.0], [30.0, 30.0], [350.0, 350.0], [35.0, 35.0]],
            axis=1,
            fill_value=0.0,
        )(freq)
    
        results[case].solvepsd(nas, case, DR, ts, F, T, freq)
        results[case].psd_data_recovery(case, DR, 1, 0)
    
        # combine results:
        toburn_combine(event, results)
    
        # save combined results:
        cla.save("results.pgz", results["combined"])
    
        # make some reports, plots:
        # SPEED: results["combined"].rpttab(excel=event.lower())
        # SPEED: results["combined"].srs_plots()
    
        # Plot SRS for Q=10 for ss, noise and combined:
        # SPEED: results["combined"].srs_plots(Q=10, direc="srs_cases", showboth=True)
    
        # Plot PSD response curves for the noise case
        # SPEED: results["noise"].resp_plots(plot="semilogy")
    
        print("Done!")


.. code-block:: none

    Mission:  Micro Space Station
    Event:    TOBurn


.. code-block:: none

    /home/docs/checkouts/readthedocs.org/user_builds/pyyeti/envs/stable/lib/python3.12/site-packages/pyyeti/srs.py:1353: RuntimeWarning: Integration frequency vector may produce inaccurate results; refine the step. See the documentation for more information.
      warn(
    /home/docs/checkouts/readthedocs.org/user_builds/pyyeti/envs/stable/lib/python3.12/site-packages/pyyeti/cla/dr_results.py:1416: RuntimeWarning: invalid value encountered in divide
      pk_freq = vrms / rms


.. code-block:: none

    Done!


For reference, here is the first page of the “TOBurn_psd.pdf” file (from
`pyyeti.cla.DR_Results.resp_plots <../modules/edited/pyyeti.cla.DR_Results.resp_plots.html#pyyeti.cla.DR_Results.resp_plots>`__):

.. figure:: clapng/toburn_pg1.png
   :alt: toburn_pg1.png

   toburn_pg1.png

TOECO
~~~~~

Outline for Transfer Orbit Engine Cutoff:

-  Load data recovery data
-  Load nastran data
-  Form ULVS for the outboard model (the SC)
-  Prepare spacecraft data recovery matrices
-  Initialize results
-  Set rfmodes
-  Setup modal mass, damping, and stiffness

   -  Damping is diagonal, 2% modal damping

-  Load in forcing functions
-  Form force transform
-  Do pre-calcs and loop over all cases
-  While looping, solve equations of motion
-  Save results and make plots

TOECO involves an acceleration recovery called “alphajoint”. The
following file “alphajoint.py” facilitates setting up this category:

::

   # alphajoint.py

   import os
   from pyyeti import cla
   from pyyeti.nastran import n2p


   def alphajoint(sol, nas, Vars, se):
       return Vars[se]["alphadrm"] @ sol.a


   def get_drdefs(nas, sc):
       drdefs = cla.DR_Def(sc["drdefs"].defaults)

       @cla.DR_Def.addcat
       def _():
           se = 0
           name = "alphajoint"
           desc = "Alpha-Joint Acceleration"
           units = "mm/sec^2, rad/sec^2"
           labels = ["Alpha-Joint {:2s}".format(i) for i in "X,Y,Z,RX,RY,RZ".split(",")]
           drms = {"alphadrm": n2p.formdrm(nas, 0, 33)[0]}
           srsopts = dict(eqsine=1, ic="steady")
           histpv = 1  # second row
           srspv = [1]
           drfile = os.path.abspath(__file__)
           drdefs.add(**locals())

       return drdefs

.. code:: ipython3

    set_event_dir("toeco")

.. code:: ipython3

    # Simulate event and recover responses
    import numpy as np
    from scipy.io import matlab
    from pyyeti import stats, ode, cla
    from pyyeti.nastran import n2p, op2
    
    import sys
    sys.path.insert(0, os.path.abspath(srcdir))
    import alphajoint
    
    
    if __name__ == "__main__":
        # event name:
        event = "TOECO"
    
        # load data recovery data:
        sc = cla.load(wrkdir / "cla_params.pgz")
        cla.PrintCLAInfo(sc["mission"], event)
    
        # load nastran data:
        nas = op2.rdnas2cam(srcdir / "nas2cam")
    
        # form ulvs for some SEs:
        SC = 101
        n2p.addulvs(nas, SC)
    
        # prepare spacecraft and alphajoint data recovery matrices
        DR = cla.DR_Event()
        DR.add(nas, sc["drdefs"])
        DR.add(nas, alphajoint.get_drdefs(nas, sc))
    
        # initialize results (ext, mnc, mxc for all drms)
        results = DR.prepare_results(sc["mission"], event)
    
        # set rfmodes:
        rfmodes = nas["rfmodes"][0]
    
        # setup modal mass, damping and stiffness
        m = None  # None means identity
        k = nas["lambda"][0]
        assert nas["nrb"] == 6
        k[: nas["nrb"]] = 0.0
        b = 2 * 0.02 * np.sqrt(k)
        mbk = (m, b, k)
    
        # load in forcing functions:
        toeco = matlab.loadmat(srcdir / "toeco_ffns.mat", squeeze_me=True, struct_as_record=False)
    
        # form force transform:
        T = n2p.formdrm(nas, 0, [[8, 12], [24, 13]])[0].T
    
        # do pre-calcs and loop over all cases:
        ts = ode.SolveUnc(*mbk, 1 / toeco["sr"], rf=rfmodes)
        LC = toeco["ffns"].shape[0]
        t = toeco["t"]
        for j, force in enumerate(toeco["ffns"]):
            print("Running {} case {}".format(event, j + 1))
            genforce = T @ ([[1], [0.1], [1], [0.1]] * force[None, :])
            # solve equations of motion
            sol = ts.tsolve(genforce, static_ic=1)
            sol.t = t
            sol = DR.apply_uf(sol, *mbk, nas["nrb"], rfmodes)
            caseid = "{} {:2d}".format(event, j + 1)
            results.time_data_recovery(sol, nas["nrb"], caseid, DR, LC, j)
    
        # save results:
        cla.save("results.pgz", results)
    
        # make some srs plots and tab files:
        # SPEED: results.rpttab()
        # SPEED: results.srs_plots()
        # SPEED: results.resp_plots()
    
        print("Done!")


.. code-block:: none

    Mission:  Micro Space Station
    Event:    TOECO
    Running TOECO case 1
    Running TOECO case 2
    Running TOECO case 3
    Done!


Take the data recovery information, spacecraft model, and forcing
function to apply the event and save the results.

Summarize Results
-----------------

.. code:: ipython3

    os.chdir(wrkdir)

.. code:: ipython3

    import numpy as np
    from pyyeti import cla
    
    
    if __name__ == "__main__":
        event = "Envelope"
    
        # load data in desired order:
        results = cla.DR_Results()
        results.merge(
            (
                cla.load(fn)
                for fn in [
                    "toes/results.pgz",
                    "owlab/results.pgz",
                    "toburn/results.pgz",
                    "toeco/results.pgz",
                ]
            ),
            {"OWLab": "O&W Lab"},
        )
    
        results.strip_hists()
        results.form_extreme(event, doappend=2)
    
        # save overall results:
        cla.save("summary_results.pgz", results)
    
        # write extrema reports:
        # SPEED: results['extreme'].rpttab()
        # SPEED: results['extreme'].rpttab(excel=event.lower())
        # SPEED: results['extreme'].srs_plots(Q=33, showall=True)
    
        print("Done!")


.. code-block:: none

    Done!


Here is the “net_ifatm.tab” file (with some lines deleted for brevity).
The extrema count table at the bottom shows that TOBurn drove most of
the “net_ifatm” extreme values.

.. figure:: clapng/ext_net_ifatm_tab.png
   :alt: ext_net_ifatm_tab.png

   ext_net_ifatm_tab.png

The
`pyyeti.cla.DR_Results.rpttab <../modules/edited/pyyeti.cla.DR_Results.rpttab.html#pyyeti.cla.DR_Results.rpttab>`__
routine can also write the results tables to an excel file. In that
case, extrema count pie charts are included. Here are the “scltm” pie
charts:

.. figure:: clapng/scltm_event_pie_chart.png
   :alt: scltm_event_pie_chart.png

   scltm_event_pie_chart.png

Grouping results:
~~~~~~~~~~~~~~~~~

The following is a brief demonstration of working with results contained
in an
`pyyeti.cla.DR_Results <../modules/edited/pyyeti.cla.DR_Results.html#pyyeti.cla.DR_Results>`__
instance.

This code groups the events into time domain events and frequency domain
events, and plots some SRS’s for assessment:

.. code:: ipython3

    # group results together to facilitate investigation:
    Grouped_Results = cla.DR_Results()
    
    # put these in the order you want:
    groups = [
       ('Time Domain', ('TOES', 'TOECO')),
       ('Freq Domain', ('O&W Lab', 'TOBurn')),
    ]
    
    for key, names in groups:
       Grouped_Results[key] = cla.DR_Results()
       for name in names:
           Grouped_Results[key][name] = results[name]
    
    Grouped_Results.form_extreme()
    
    # plot just time domain srs:
    # SPEED: Grouped_Results['Time Domain']['extreme'].srs_plots(
    # SPEED:    direc='timedomain_srs', Q=33, showall=True
    # SPEED: )
    
    # plot the srs of the two groups together:
    # SPEED: Grouped_Results['extreme'].srs_plots(
    # SPEED:    direc='grouped_srs', Q=33, showall=True
    # SPEED: )
    
    print("Done!")


.. code-block:: none

    Done!


Here is a page from the “grouped_srs/Envelope_srs.pdf” file:

.. figure:: clapng/grouped_srs_ifatm_pg.png
   :alt: grouped_srs_ifatm_pg.png

   grouped_srs_ifatm_pg.png

Compare Results
---------------

Compare results generated here to those from the contractor:

.. code:: ipython3

    import numpy as np
    from pyyeti import cla
    import pandas as pd
    
    
    if __name__ == "__main__":
        results = cla.load("summary_results.pgz")
        lvc = pd.read_excel(srcdir / "contractor_results.xlsx", sheet_name=None, index_col=0)
        # SPEED: results['extreme'].rptpct(lvc, names=("Ours", "Theirs"), direc='compare')
    
        print("Done!")


.. code-block:: none

    Done!


For each category,
`pyyeti.cla.DR_Results.rptpct <../modules/edited/pyyeti.cla.DR_Results.rptpct.html#pyyeti.cla.DR_Results.rptpct>`__
creates three files in the “compare” subdirectory: a *.cmp file,
a*.cmp.histogram.png file, and a \*.cmp.magpct.png file. Examples from
the “net_ifatm” category follow.

Here is a the top of the “net_ifatm.cmp” file. The table at the top has
the detailed comparison of maximums, minimums, and absolute-maximums.
There are accompanying comparison histogram counts with statistics for
each of these at the bottom (here, for brevity, only the maximums
histogram data is shown). All “net_ifatm” maximums are within 3%, and
75% of maximums are within 1%.

.. figure:: clapng/cmp_net_ifatm_tab.png
   :alt: cmp_net_ifatm_tab.png

   cmp_net_ifatm_tab.png

Here is the “net_ifatm.cmp.histogram.png” file which shows the histogram
data in a bar plot format. Bars within 5% are shown in blue, bars from 5
to 10% are purple, and bars above 10% are red.

.. figure:: clapng/net_ifatm.cmp.histogram.png
   :alt: net_ifatm.cmp.histogram.png

   net_ifatm.cmp.histogram.png

Here is the “net_ifatm.cmp.magpct.png” file which is a scatter plot
showing the comparisons of each data point versus the magnitude of
“Theirs”. It uses the same color scheme as listed for the histogram data
in a bar plot format. This plot gives a quick view of how the comparison
looks over the range of small numbers to large numbers.

.. figure:: clapng/net_ifatm.cmp.magpct.png
   :alt: net_ifatm.cmp.magpct.png

   net_ifatm.cmp.magpct.png

Since the above “magpct” plot is not that interesting, here is the
“scltm.cmp.magpct.png” file. The shaded out areas are for values smaller
than the filter (which was left at the default +/- 1e-6 and shown by the
dashed vertical lines), so our focus should be on the white, “Filtered
region”. We can quickly see that there are values greater than 1000 that
are off by more than 10%. More investigation would be needed to
disposition these items. The specific rows can be determined from
looking at the “scltm.cmp” file. In this case, all these differences
were for CBAR bending moments and the only difference in the two sets of
results is the version of Nastran used. (A deeper cause was not sought.)

.. figure:: clapng/scltm.cmp.magpct.png
   :alt: scltm.cmp.magpct.png

   scltm.cmp.magpct.png