new_cross_sectional_properties.py

import unittest

import numpy as np
import openmdao.api as om

from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials


class CrossSectionProperties(om.ExplicitComponent):
    """ Calculates the cross-section area and area moments of inertia. """

    def initialize(self):
        self.options.declare(
            "n_keypoints",
            default=2,
            desc="Number of wing reference axis keypoints",
        )

    def setup(self):
        n = self.options["n_keypoints"]

        # Inputs
        self.add_input(
            "h_spar",
            val=np.ones(n),
            units="m",
            desc="Wingbox spar height distribution",
        )
        self.add_input(
            "t_spar",
            val=1e-3 * np.ones(n),
            units="m",
            desc="Wingbox spar thickness distribution",
        )
        self.add_input(
            "w_skin",
            val=np.ones(n),
            units="m",
            desc="Wingbox skin width distribution",
        )
        self.add_input(
            "t_skin",
            val=1e-3 * np.ones(n),
            units="m",
            desc="Wingbox skin thickenss distribution",
        )

        # Outputs
        self.add_output("A", val=np.ones(n), units="m**2", desc="Wingbox area")
        self.add_output(
            "Iyy",
            val=np.ones(n),
            units="m**4",
            desc="Wingbox area moment of inertia around y axis",
        )
        self.add_output(
            "Izz",
            val=np.ones(n),
            units="m**4",
            desc="Wingbox area moment of inertia around z axis",
        )
        self.add_output(
            "J", val=np.ones(n), units="m**4", desc="Wingbox torsion constant"
        )

        # Partials
        row_col = np.arange(n)
        self.declare_partials("*", "*", rows=row_col, cols=row_col)

    def compute(self, inputs, outputs):
        """ Calculates the module outputs. """

        # get inputs
        h_spar = inputs["h_spar"]
        t_spar = inputs["t_spar"]
        w_skin = inputs["w_skin"]
        t_skin = inputs["t_skin"]

        # calculate area
        outputs["A"] = 2.0 * (h_spar * t_spar + w_skin * t_skin)

        # calculate moments of inertia
        outputs["Iyy"] = (
            h_spar ** 3 * t_spar / 3.0
            + w_skin * t_skin ** 3 / 3.0
            + w_skin * t_skin * h_spar ** 2
        ) / 2.0
        outputs["Izz"] = (
            t_spar ** 3 * h_spar / 3.0
            + t_skin * w_skin ** 3 / 3.0
            + h_spar * t_spar * w_skin ** 2
        ) / 2.0

        # calculate torsion constant
        outputs["J"] = (
            2.0
            * w_skin ** 2
            * h_spar ** 2
            * t_skin
            * t_spar
            / (w_skin * t_spar + h_spar * t_skin)
        )

    def compute_partials(self, inputs, partials):
        """ Calculates the module partials. """

        # get inputs
        h_spar = inputs["h_spar"]
        t_spar = inputs["t_spar"]
        w_skin = inputs["w_skin"]
        t_skin = inputs["t_skin"]

        # calculate area partials
        partials["A", "h_spar"] = 2.0 * t_spar
        partials["A", "t_spar"] = 2.0 * h_spar
        partials["A", "w_skin"] = 2.0 * t_skin
        partials["A", "t_skin"] = 2.0 * w_skin

        # calculate moment of inertia around y partials
        partials["Iyy", "h_spar"] = (t_spar * h_spar ** 2) / 2.0 + t_skin * w_skin * h_spar
        partials["Iyy", "t_spar"] = h_spar ** 3 / 6.0
        partials["Iyy", "w_skin"] = t_skin * (3.0 * h_spar ** 2 + t_skin ** 2) / 6.0
        partials["Iyy", "t_skin"] = w_skin * (h_spar ** 2 + t_skin ** 2) / 2.0

        # calculate moment of inertia around z partials
        partials["Izz", "h_spar"] = t_spar * (3.0 * w_skin ** 2 + t_spar ** 2) / 6.0
        partials["Izz", "t_spar"] = h_spar * (t_spar ** 2 + w_skin ** 2) / 2.0
        partials["Izz", "w_skin"] = (t_skin * w_skin ** 2) / 2.0 + h_spar * t_spar * w_skin
        partials["Izz", "t_skin"] = w_skin ** 3 / 6.0

        # calculate denominator of following partials
        dJ_den = (h_spar * t_skin + t_spar * w_skin) ** 2

        # calculate torsion constant partials
        partials["J", "h_spar"] = 2.0 * h_spar * t_skin * t_spar * w_skin ** 2 \
            * (2.0 * w_skin * t_spar + h_spar * t_skin) \
            / dJ_den

        partials["J", "t_spar"] = 2.0 * h_spar ** 3 * t_skin ** 2 * w_skin ** 2 / dJ_den
        partials["J", "w_skin"] = 2.0 * w_skin * t_skin * t_spar * h_spar ** 2 \
            * (2.0 * h_spar * t_skin + w_skin * t_spar) \
            / dJ_den
        partials["J", "t_skin"] = 2.0 * w_skin ** 3 * t_spar ** 2 * h_spar ** 2 / dJ_den


def assert_partials_match_fd(model, inputs, atol=1e-6, rtol=1e-6, method="fd"):
    prob = om.Problem()
    prob.model = model
    prob.set_solver_print(level=0)
    if method=='cs':
        prob.setup(force_alloc_complex=True)
    else:
        prob.setup()

    # Set inputs
    for key, value in inputs.items():
        prob.set_val(key, value)

    prob.run_model()

    data = prob.check_partials(method=method, out_stream=None)
    assert_check_partials(data, atol, rtol)


class TestCrossSectionProperties(unittest.TestCase):
    def setUp(self):

        self.n_keypoints = 2

        # Create a CrossSectionProperties object
        self.csp = CrossSectionProperties(n_keypoints=self.n_keypoints)

        # Set inputs
        self.inputs = {
            "h_spar": np.asarray([0.1, 0.2]),
            "w_skin": np.asarray([0.5, 0.6]),
            "t_spar": np.asarray([0.002, 0.003]),
            "t_skin": np.asarray([0.001, 0.0005]),
        }

    def test_compute(self):
        outputs = {}
        self.csp.compute(self.inputs, outputs)

        # Unpack inputs
        h_spar = self.inputs["h_spar"]
        t_spar = self.inputs["t_spar"]
        w_skin = self.inputs["w_skin"]
        t_skin = self.inputs["t_skin"]

        # Set reference outputs
        outputs_ref = {
            "A": 2.0 * (h_spar * t_spar + w_skin * t_skin),
            "Iyy": 0.5
            * (
                h_spar ** 3 * t_spar / 3.0
                + w_skin * t_skin ** 3 / 3.0
                + w_skin * t_skin * h_spar ** 2
            ),
            "Izz": 0.5
            * (
                h_spar * t_spar ** 3 / 3.0
                + w_skin ** 3 * t_skin / 3.0
                + h_spar * t_spar * w_skin ** 2
            ),
            "J": 2.0
            * w_skin ** 2
            * h_spar ** 2
            * t_skin
            * t_spar
            / (w_skin * t_spar + h_spar * t_skin),
        }

        # Check outputs
        assert_near_equal(outputs, outputs_ref)

    def test_partials_against_fd(self):
        assert_partials_match_fd(
            self.csp, self.inputs, method="fd", rtol=1e-3, atol=1e-3
        )

    def test_partials_against_cs(self):
        assert_partials_match_fd(
            self.csp, self.inputs, method="cs", rtol=1e-3, atol=1e-5
        )


if __name__ == "__main__":
    unittest.main()