Merge pull request #139 from kejacobson/time_domain

First step of generalizing the time domain integrator component to mphys builders

.gitignore

@@ -61,3 +61,5 @@ tests/constant
 tests/FFD
 tests/system
 
+# OM output
+reports

examples/time_domain/README.md

@@ -1,6 +1,26 @@
-This directory contains exploratory work looking at different ways of doing time integration in OpenMDAO
+# Time-domain Prototyping
 
-timestep_groups: each time step is a OpenMDAO group that connected within the model.
-This approach requires all of the states to be allocated and stored in memory
+Note: these demonstrations are initial proofs of concept and have not been exercised with actual high-fidelity codes. Time-domain infrastructure and standards in MPhys may need to be changed as issues appear when actual codes are integrated.
 
-timestep_loop: time steps are solved in a single explicit component which loops over an OpenMDAO problem that represents the time steps
+This directory contains exploratory work looking at two different ways of doing time integration in OpenMDAO.
+
+
+1. `timestep_groups`: The first mode is the default mode of OpenMDAO where each time step of the temporal problem is a subsystem and connections are used to indicate the dependency of a particular time step on previous time
+This requires OpenMDAO to allocate and store all states in the system in memory at all times;
+this is not practical with large high-fidelity codes where there could be hundreds or thousands of time steps.
+
+2. `timestep_loop`: The alternative approach is to use nested OpenMDAO problems.
+The outer OpenMDAO problem is the optimization problem, and the inner problem represents a single coupled time step of the time-domain problem.
+There is an explicit component in the outer problem that calls the inner problem in a time step loop in order to solve the time-domain problem.
+This explicit component does not expose the temporally varying states to the outer problem, and it manages the temporal data in such a way that it doesn't need to keep all of it in memory at the same time, i.e., it can write and read states to/from disk as needed in order to solve the temporal nonlinear, forward, or reverse problems.
+
+These two modes are demonstrated in the `prototype` directory which models a simplified structural dynamics problem with a weakly coupled forcing component.
+The purpose of these prototypes is to demonstrate that the `timestep_loop` method will get the same answer as the default `timestep_loop` approach while not having to expose every temporal state to the outer OpenMDAO problem.
+
+While the prototype shows that the `timestep_loop` approach is feasible, the temporal integrator component is specific to the structural problem being solved.
+Therefore, the `builder_version` looks to generalize this timestep loop approach by introducing unsteady versions of the MPhys Builders and Scenarios called the Integrator.
+The `TimeDomainBuilder` extends the standard MPhys Builder class to provide additional information about the temporal integration such as how many temporal backplanes of data are required for a particular state.
+The Integrator component performs the time step loops and manages the temporal state data including shuffling backplanes of states as the time is advanced.
+In this simplified example, the aerodynamic loads are computed on the same coordinate locations as the structural mesh.
+Since a normal load and displacement transfer scheme is not required,
+the example's "load and displacement transfer components" are the modal force and displacement conversions in order to fit a aerostructural coupling group pattern.

examples/time_domain/builder_version/aero_solver.py

@@ -0,0 +1,117 @@
+#!/usr/bin/env python
+import numpy as np
+import openmdao.api as om
+from mphys.time_domain.time_domain_builder import TimeDomainBuilder
+from mphys.time_domain.time_domain_variables import TimeDomainInput
+
+
+class HarmonicForcer(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare("number_of_nodes")
+        self.options.declare("c1", default=1e-3)
+
+    def setup(self):
+        self.nnodes = self.options["number_of_nodes"]
+        self.c1 = self.options["c1"]
+
+        self.add_input(
+            "amplitude_aero",
+            shape_by_conn=True,
+            tags=["mphys_input"],
+            desc="amplitude_aero",
+        )
+        self.add_input(
+            "freq_aero", shape_by_conn=True, tags=["mphys_input"], desc="frequency"
+        )
+        self.add_input("time", tags=["mphys_input"], desc="current time")
+        self.add_input("x_aero", shape_by_conn=True, tags=["mphys_coupling"])
+        self.add_output("f_aero", shape=self.nnodes * 3, tags=["mphys_coupling"])
+
+    def compute(self, inputs, outputs):
+        amp = inputs["amplitude_aero"]
+        freq = inputs["freq_aero"]
+        time = inputs["time"]
+
+        outputs["f_aero"] = (
+            amp * np.sin(freq * time) * np.ones(self.nnodes * 3)
+            - self.c1 * inputs["x_aero"]
+        )
+
+    def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
+        amp = inputs["amplitude_aero"]
+        freq = inputs["freq_aero"]
+        time = inputs["time"]
+
+        if mode == "fwd":
+            if "f_aero" in d_outputs:
+                if "amplitude_aero" in d_inputs:
+                    d_outputs["f_aero"] += (
+                        np.sin(freq * time)
+                        * np.ones(self.nnodes * 3)
+                        * d_inputs["amplitude_aero"]
+                    )
+                if "freq_aero" in d_inputs:
+                    d_outputs["f_aero"] += (
+                        time
+                        * np.cos(freq * time)
+                        * np.ones(self.nnodes * 3)
+                        * d_inputs["freq_aero"]
+                    )
+                if "time" in d_inputs:
+                    d_outputs["f_aero"] += (
+                        freq
+                        * np.cos(freq * time)
+                        * np.ones(self.nnodes * 3)
+                        * d_inputs["time"]
+                    )
+                if "x_aero" in d_inputs:
+                    d_outputs["f_aero"] -= self.c1 * d_inputs["x_aero"]
+
+        if mode == "rev":
+            if "f_aero" in d_outputs:
+                if "amplitude_aero" in d_inputs:
+                    d_inputs["amplitude_aero"] += np.sin(freq * time) * np.sum(
+                        d_outputs["f_aero"]
+                    )
+                if "freq_aero" in d_inputs:
+                    d_inputs["freq_aero"] += (
+                        time * np.cos(freq * time) * np.sum(d_outputs["f_aero"])
+                    )
+                if "time" in d_inputs:
+                    d_inputs["time"] += (
+                        freq * np.cos(freq * time) * np.sum(d_outputs["f_aero"])
+                    )
+                if "x_aero" in d_inputs:
+                    d_inputs["x_aero"] -= self.c1 * d_outputs["f_aero"]
+
+
+class FakeAeroBuilder(TimeDomainBuilder):
+    def __init__(self, root_name, nmodes, dt):
+        self.nmodes = nmodes
+        self.dt = dt
+        self._read_number_of_nodes(root_name)
+
+    def _read_number_of_nodes(self, root_name):
+        filename = f"{root_name}_mode1.dat"
+        fh = open(filename)
+        while True:
+            line = fh.readline()
+            if "zone" in line.lower():
+                self.nnodes = int(line.split("=")[2].split(",")[0])
+                return
+
+    def get_number_of_nodes(self):
+        return self.nnodes
+
+    def get_ndof(self):
+        return 1
+
+    def get_coupling_group_subsystem(self, scenario_name=None):
+        return HarmonicForcer(number_of_nodes=self.nnodes)
+
+    def get_timestep_input_variables(self, scenario_name=None) -> list[TimeDomainInput]:
+        return [
+            TimeDomainInput("amplitude_aero", (1)),
+            TimeDomainInput("freq_aero", (1)),
+            TimeDomainInput("x_aero0", (self.nnodes * 3), True),
+        ]

examples/time_domain/builder_version/run.py

@@ -0,0 +1,62 @@
+import numpy as np
+import openmdao.api as om
+
+from aero_solver import FakeAeroBuilder
+from xfer_scheme import ModalXferBuilder
+from struct_solver import ModalStructBuilder
+from mphys.time_domain.integator_aerostructural import IntegratorAerostructural
+
+
+class Model(om.Group):
+    def setup(self):
+        nsteps = 10
+        dt = 1.0
+        nmodes = 2
+        root_name = "sphere_body1"
+
+        aero_builder = FakeAeroBuilder(root_name, nmodes, dt)
+        aero_builder.initialize(self.comm)
+
+        struct_builder = ModalStructBuilder(nmodes, dt)
+        struct_builder.initialize(self.comm)
+
+        ldxfer_builder = ModalXferBuilder(root_name, nmodes)
+        ldxfer_builder.initialize(self.comm)
+
+        dvs = om.IndepVarComp()
+        dvs.add_output("amplitude_aero", 1.0)
+        dvs.add_output("freq_aero", 0.1)
+        dvs.add_output("m", [1.0, 1.0])
+        dvs.add_output("c", [0.0, 0.0])
+        dvs.add_output("k", [1.0, 1.5])
+        self.add_subsystem("design_variables", dvs, promotes=["*"])
+
+        inputs = om.IndepVarComp()
+        inputs.add_output("u_struct|0", [1.0, 2.0], distributed=True)
+        inputs.add_output("x_aero0", np.ones(aero_builder.nnodes * 3), distributed=True)
+        self.add_subsystem("analysis_inputs", inputs, promotes=["*"])
+
+        self.add_subsystem(
+            "integrator",
+            IntegratorAerostructural(
+                nsteps=nsteps,
+                dt=dt,
+                aero_builder=aero_builder,
+                struct_builder=struct_builder,
+                ldxfer_builder=ldxfer_builder,
+                nonlinear_solver = om.NonlinearRunOnce(),
+                linear_solver = om.LinearRunOnce()
+            ),
+            promotes=["*"],
+        )
+
+
+def main():
+    prob = om.Problem()
+    prob.model = Model()
+    prob.setup()
+    prob.run_model()
+
+
+if __name__ == "__main__":
+    main()