add coupled unsteady problem and solver

pterasoftware/__init__.py

@@ -19,7 +19,8 @@
 
 output.py: Contains functions for visualizing geometry and results.
 
-problems.py: Contains the SteadyProblem and UnsteadyProblem classes.
+problems.py: Contains the SteadyProblem, UnsteadyProblem, and CoupledUnsteadyProblem
+classes.
 
 steady_horseshoe_vortex_lattice_method.py: Contains the
 SteadyHorseshoeVortexLatticeMethodSolver class.
@@ -33,6 +34,9 @@
 unsteady_ring_vortex_lattice_method.py: Contains the
 UnsteadyRingVortexLatticeMethodSolver class.
 
+coupled_unsteady_ring_vortex_lattice_method.py: Contains the
+CoupledUnsteadyRingVortexLatticeMethodSolver class.
+
 **Contains the following functions:**
 
 load: Loads a Ptera Software object from a JSON file.
@@ -52,6 +56,7 @@
 # Lazy imports configuration: modules loaded on first access.
 _LAZY_MODULES = {
     "convergence": "pterasoftware.convergence",
+    "coupled_unsteady_ring_vortex_lattice_method": "pterasoftware.coupled_unsteady_ring_vortex_lattice_method",
     "output": "pterasoftware.output",
     "steady_horseshoe_vortex_lattice_method": "pterasoftware.steady_horseshoe_vortex_lattice_method",
     "steady_ring_vortex_lattice_method": "pterasoftware.steady_ring_vortex_lattice_method",

pterasoftware/_core.py

@@ -5,9 +5,13 @@
 import copy
 import math
 from collections.abc import Callable, Sequence
+from typing import TYPE_CHECKING
 
 import numpy as np
 
+if TYPE_CHECKING:
+    from . import problems as problems_mod
+
 from . import _oscillation, _parameter_validation, _transformations, geometry
 from . import operating_point as operating_point_mod
 
@@ -2361,3 +2365,11 @@ def first_results_step(self) -> int:
     @property
     def max_wake_rows(self) -> int | None:
         return self._max_wake_rows
+
+    @property
+    def movement(self) -> CoreMovement:
+        raise NotImplementedError
+
+    @property
+    def steady_problems(self) -> tuple[problems_mod.SteadyProblem, ...]:
+        raise NotImplementedError

pterasoftware/problems.py

@@ -1,23 +1,33 @@
-"""Contains the SteadyProblem and UnsteadyProblem classes.
+"""Contains the SteadyProblem, UnsteadyProblem, and CoupledUnsteadyProblem classes.
 
 **Contains the following classes:**
 
 SteadyProblem: A class used to contain steady aerodynamics problems.
 
 UnsteadyProblem: A class used to contain unsteady aerodynamics problems.
 
+CoupledUnsteadyProblem: A class used to contain coupled unsteady aerodynamics problems
+where SteadyProblems are generated dynamically at each time step.
+
 **Contains the following functions:**
 
 None
 """
 
 from __future__ import annotations
 
+from typing import TYPE_CHECKING
+
 import numpy as np
 
-from . import _core, _transformations, geometry, movements
+from . import _core, _parameter_validation, _transformations, geometry, movements
 from . import operating_point as operating_point_mod
 
+if TYPE_CHECKING:
+    from .unsteady_ring_vortex_lattice_method import (
+        UnsteadyRingVortexLatticeMethodSolver,
+    )
+
 
 class SteadyProblem:
     """A class used to contain steady aerodynamics problems.
@@ -229,3 +239,152 @@ def movement(self) -> movements.movement.Movement:
     @property
     def steady_problems(self) -> tuple[SteadyProblem, ...]:
         return self._steady_problems
+
+
+class CoupledUnsteadyProblem(_core.CoreUnsteadyProblem):
+    """A class used to contain coupled unsteady aerodynamics problems.
+
+    Unlike UnsteadyProblem, which pre-computes all SteadyProblems during initialization,
+    CoupledUnsteadyProblem generates SteadyProblems dynamically at each time step. This
+    enables two-way coupling between the aerodynamic solver and external models (e.g.,
+    structural deformation, rigid body dynamics) where the geometry or operating
+    conditions at step N + 1 depend on the aerodynamic solution at step N.
+
+    The base implementation generates each step's SteadyProblem from the movement's
+    prescribed geometry and operating conditions. Subclasses
+    (AeroelasticUnsteadyProblem, FreeFlightUnsteadyProblem) override
+    initialize_next_problem to incorporate feedback from the solver.
+
+    **Contains the following methods:**
+
+    movement: The CoreMovement that contains this CoupledUnsteadyProblem's
+    OperatingPointMovement and AirplaneMovements.
+
+    steady_problems: A tuple of the SteadyProblems generated so far.
+
+    get_steady_problem: Retrieves the SteadyProblem at a given time step.
+
+    initialize_next_problem: Generates and appends the next time step's SteadyProblem.
+    """
+
+    __slots__ = (
+        "_movement",
+        "_coupled_steady_problems",
+    )
+
+    def __init__(
+        self,
+        movement: _core.CoreMovement,
+        only_final_results: bool | np.bool_ = False,
+    ) -> None:
+        """The initialization method.
+
+        :param movement: The CoreMovement (or subclass) that contains this
+            CoupledUnsteadyProblem's OperatingPointMovement and AirplaneMovements. It
+            must be an instance of CoreMovement.
+        :param only_final_results: Determines whether the solver will only calculate
+            loads for the final time step (for static movements) or will only calculate
+            loads for the time steps in the final complete motion cycle (for non-static
+            movements), which increases simulation speed. Can be a bool or a numpy bool
+            and will be converted internally to a bool. The default is False.
+        :return: None
+        """
+        if not isinstance(movement, _core.CoreMovement):
+            raise TypeError("movement must be a CoreMovement or one of its subclasses.")
+        self._movement = movement
+
+        # Delegate shared initialization (validation, first_averaging_step computation,
+        # load list initialization) to the core class.
+        super().__init__(
+            only_final_results=only_final_results,
+            delta_time=self._movement.delta_time,
+            num_steps=self._movement.num_steps,
+            max_wake_rows=self._movement.max_wake_rows,
+            lcm_period=self._movement.lcm_period,
+        )
+
+        # Generate the initial SteadyProblem for step 0 from the movement's prescribed
+        # geometry and operating conditions.
+        initial_airplanes: list[geometry.airplane.Airplane] = []
+        for airplane_movement in self._movement.airplane_movements:
+            initial_airplanes.append(
+                airplane_movement.generate_airplane_at_time_step(
+                    0, self._movement.delta_time
+                )
+            )
+        initial_operating_point = self._movement.operating_point_movement.generate_operating_point_at_time_step(
+            0, self._movement.delta_time
+        )
+        initial_problem = SteadyProblem(
+            airplanes=initial_airplanes,
+            operating_point=initial_operating_point,
+        )
+
+        # Initialize the mutable list of SteadyProblems with the initial problem.
+        self._coupled_steady_problems: list[SteadyProblem] = [initial_problem]
+
+    # --- Immutable: read only properties ---
+    @property
+    def movement(self) -> _core.CoreMovement:
+        return self._movement
+
+    @property
+    def steady_problems(self) -> tuple[SteadyProblem, ...]:
+        """A tuple of the SteadyProblems generated so far.
+
+        During simulation, this tuple grows as initialize_next_problem is called at each
+        time step. After simulation completes, it contains one SteadyProblem per time
+        step.
+
+        :return: A tuple of SteadyProblems.
+        """
+        return tuple(self._coupled_steady_problems)
+
+    def get_steady_problem(self, step: int) -> SteadyProblem:
+        """Retrieves the SteadyProblem at a given time step.
+
+        :param step: An int representing the time step index. It must be non-negative
+            and less than the number of SteadyProblems generated so far.
+        :return: The SteadyProblem at the given time step.
+        """
+        if step < 0 or step >= len(self._coupled_steady_problems):
+            raise ValueError(
+                f"Step index {step} is out of range. Only"
+                f" {len(self._coupled_steady_problems)} problems have been"
+                " initialized so far."
+            )
+        return self._coupled_steady_problems[step]
+
+    def initialize_next_problem(
+        self, solver: UnsteadyRingVortexLatticeMethodSolver
+    ) -> None:
+        """Generates and appends the next time step's SteadyProblem.
+
+        The base implementation generates the next SteadyProblem from the movement's
+        prescribed geometry and operating conditions. Subclasses override this method to
+        incorporate feedback from the solver (e.g., structural deformation angles,
+        dynamics-integrated operating conditions).
+
+        :param solver: The solver instance, which provides access to the current
+            aerodynamic solution state. Subclasses use this to extract forces, moments,
+            or other quantities needed to compute the next step's geometry or operating
+            conditions.
+        :return: None
+        """
+        next_step = len(self._coupled_steady_problems)
+
+        next_airplanes: list[geometry.airplane.Airplane] = []
+        for airplane_movement in self._movement.airplane_movements:
+            next_airplanes.append(
+                airplane_movement.generate_airplane_at_time_step(
+                    next_step, self._movement.delta_time
+                )
+            )
+        next_operating_point = self._movement.operating_point_movement.generate_operating_point_at_time_step(
+            next_step, self._movement.delta_time
+        )
+        next_problem = SteadyProblem(
+            airplanes=next_airplanes,
+            operating_point=next_operating_point,
+        )
+        self._coupled_steady_problems.append(next_problem)

pterasoftware/unsteady_ring_vortex_lattice_method.py

@@ -21,6 +21,7 @@
 
 from . import (
     _aerodynamics_functions,
+    _core,
     _functions,
     _logging,
     _panel,
@@ -127,14 +128,20 @@ class UnsteadyRingVortexLatticeMethodSolver:
         "ran",
     )
 
-    def __init__(self, unsteady_problem: problems.UnsteadyProblem) -> None:
+    def __init__(self, unsteady_problem: _core.CoreUnsteadyProblem) -> None:
         """The initialization method.
 
-        :param unsteady_problem: The UnsteadyProblem to be solved.
+        :param unsteady_problem: The CoreUnsteadyProblem to be solved. This can be an
+            UnsteadyProblem (with pre-computed SteadyProblems) or a
+            CoupledUnsteadyProblem (with dynamically generated SteadyProblems). The
+            problem must expose a steady_problems attribute and a movement property.
         :return: None
         """
-        if not isinstance(unsteady_problem, problems.UnsteadyProblem):
-            raise TypeError("unsteady_problem must be an UnsteadyProblem.")
+        if not isinstance(unsteady_problem, _core.CoreUnsteadyProblem):
+            raise TypeError(
+                "unsteady_problem must be a CoreUnsteadyProblem or one of its"
+                " subclasses."
+            )
         self.unsteady_problem = unsteady_problem
 
         self._max_wake_rows = self.unsteady_problem.max_wake_rows
@@ -1449,6 +1456,54 @@ def _calculate_loads(self) -> None:
             + unsteady_forces_GP1
         )
 
+        # Calculate the moments via the hook method. The base implementation
+        # computes moments in the first Airplane's geometry axes, relative to the
+        # first Airplane's CG. Subclasses can override the hook to compute
+        # additional moment representations (e.g., about strip leading edge points).
+        moments_GP1_CgP1 = self._load_calculation_moment_processing_hook(
+            rightLegForces_GP1,
+            frontLegForces_GP1,
+            leftLegForces_GP1,
+            backLegForces_GP1,
+            unsteady_forces_GP1,
+        )
+
+        # TODO: Transform forces_GP1 and moments_GP1_CgP1 to each Airplane's local
+        #  geometry axes before passing to process_solver_loads.
+        _functions.process_solver_loads(self, forces_GP1, moments_GP1_CgP1)
+
+    def _load_calculation_moment_processing_hook(
+        self,
+        rightLegForces_GP1: np.ndarray,
+        frontLegForces_GP1: np.ndarray,
+        leftLegForces_GP1: np.ndarray,
+        backLegForces_GP1: np.ndarray,
+        unsteady_forces_GP1: np.ndarray,
+    ) -> np.ndarray:
+        """Computes moments (in the first Airplane's geometry axes, relative to the
+        first Airplane's CG) from the per-leg and unsteady forces on every Panel.
+
+        Subclasses can override this method to compute additional moment representations
+        (e.g., about strip leading edge points) while still returning the CG-based
+        moments for the standard load processing pipeline.
+
+        :param rightLegForces_GP1: An (N, 3) ndarray of floats representing the forces
+            (in the first Airplane's geometry axes) on each Panel's right LineVortex
+            leg.
+        :param frontLegForces_GP1: An (N, 3) ndarray of floats representing the forces
+            (in the first Airplane's geometry axes) on each Panel's front LineVortex
+            leg.
+        :param leftLegForces_GP1: An (N, 3) ndarray of floats representing the forces
+            (in the first Airplane's geometry axes) on each Panel's left LineVortex leg.
+        :param backLegForces_GP1: An (N, 3) ndarray of floats representing the forces
+            (in the first Airplane's geometry axes) on each Panel's back LineVortex leg.
+        :param unsteady_forces_GP1: An (N, 3) ndarray of floats representing the
+            unsteady component of the force (in the first Airplane's geometry axes) on
+            each Panel, derived from the unsteady Bernoulli equation.
+        :return: An (N, 3) ndarray of floats representing the moments (in the first
+            Airplane's geometry axes, relative to the first Airplane's CG) on every
+            Panel at the current time step.
+        """
         # Find the moments (in the first Airplane's geometry axes, relative to the
         # first Airplane's CG) on the Panels' RingVortex's right LineVortex,
         # front LineVortex, left LineVortex, and back LineVortex.
@@ -1470,23 +1525,21 @@ def _calculate_loads(self) -> None:
         # collocation point, not at the Panel's centroid.
 
         # Find the moments (in the first Airplane's geometry axes, relative to the
-        # first Airplane's CG) due to the unsteady component of the force on each Panel.
+        # first Airplane's CG) due to the unsteady component of the force on each
+        # Panel.
         unsteady_moments_GP1_CgP1 = _functions.numba_1d_explicit_cross(
             self.stackCpp_GP1_CgP1, unsteady_forces_GP1
         )
 
-        moments_GP1_CgP1 = (
+        return cast(
+            np.ndarray,
             rightLegMoments_GP1_CgP1
             + frontLegMoments_GP1_CgP1
             + leftLegMoments_GP1_CgP1
             + backLegMoments_GP1_CgP1
-            + unsteady_moments_GP1_CgP1
+            + unsteady_moments_GP1_CgP1,
         )
 
-        # TODO: Transform forces_GP1 and moments_GP1_CgP1 to each Airplane's local
-        #  geometry axes before passing to process_solver_loads.
-        _functions.process_solver_loads(self, forces_GP1, moments_GP1_CgP1)
-
     def _populate_next_airplanes_wake(self) -> None:
         """Updates the next time step's Airplanes' wakes.
 
@@ -2084,8 +2137,7 @@ def _finalize_loads(self) -> None:
         num_steps_to_average = self.num_steps - self._first_averaging_step
 
         # Determine if this SteadyProblem's geometry is static or variable.
-        this_movement: movements.movement.Movement = self.unsteady_problem.movement
-        static = this_movement.static
+        static = self.unsteady_problem.movement.static
 
         # Initialize ndarrays to hold each Airplane's loads and load coefficients at
         # each of the time steps that calculated the loads.