CoSAppLogo CoSApp examples

Tips & Tricks

Purge all drivers from a system

Often necessary to start off a new, clean simulation workflow.

Drivers are stored in a plain dictionary drivers. Thus, cleaning all drivers from system s is as simple as

s.drivers.clear()

Save and restore system state

You may quickly save and restore the state of a system in an execution sequence using functions get_state and set_state from module cosapp.utils:

from cosapp.utils import get_state, set_state

s = SomeSystem('s')
# ... many design steps, say

# Save state in local object
designed_state = get_state(s)

s.drivers.clear()
s.add_driver(SomeDriver('driver'))

try:
    s.run_drivers()
except:
    # Recover previous state
    set_state(s, designed_state)

Note that the objet returned by get_state contains the inner data of the system (both inputs and outputs), but does not store the description of the system (class names, sub-systems, connections, etc.). Data are transfered by deep (recursive) copy, to avoid all side effects between the original system and the extracted state.

Consistently, set_state does not perform any structure check, and resets system state by deep copy operations.

Set targets on output values

See tutorial on method add_target.

Display the internal off-design problem of a system

The overall off-design problem of a system, gathered from all unknowns and equations defined at setup throughout the sub-system tree, can be extracted from method System.assembled_problem. By default, additional unknowns and equations arising from cyclic loop analysis (if any) are not accounted for in the assembled problem. To include them, one must call method open_loops() before problem assembly.

In the next cell, we show a slight variation of the Sellar case, with off-design constraints on parameters y1, y2 and z:

[1]:
from cosapp.base import System
from cosapp.drivers import NonLinearSolver
import numpy as np


class SellarDiscipline1(System):
    """y1 = F1(x, y2, z)"""
    def setup(self):
        self.add_inward('z', np.zeros(2))
        self.add_inward('x', 1.0)
        self.add_inward('y2', 1.0)

        self.add_outward('y1', 0.0)

    def compute(self):
        self.y1 = self.z[0]**2 + self.z[1] + self.x - 0.2 * self.y2


class SellarDiscipline2(System):
    """y2 = F2(y1, z)"""
    def setup(self):
        self.add_inward('z', np.zeros(2))
        self.add_inward('y1', 1.0)

        self.add_outward('y2', 0.0)

    def compute(self):
        self.y2 = np.sqrt(abs(self.y1)) + self.z[0] + self.z[1]


class ConstrainedSellar(System):
    """System coupling disciplines 1 & 2, with off-design constraints.
    """
    def setup(self):
        d1 = self.add_child(SellarDiscipline1('d1'), pulling=['x', 'z', 'y1'])
        d2 = self.add_child(SellarDiscipline2('d2'), pulling=['z', 'y2'])

        # Couple sub-systems d1 and d2:
        self.connect(d1, d2, ['y1', 'y2'])

        # Add off-design constraints:
        self.add_unknown('z').add_equation(['y1 == 1', 'y2 == 0.5'])


def show_problem(problem, header=None):
    if header:
        print(f"*** {header.upper()} ***")
    print(f"shape = {problem.shape}", problem, "", sep="\n")

# Create and initialize system
s = ConstrainedSellar('s')
s.x = 0.1
s.z = np.zeros(2)
s.run_once()

show_problem(s.assembled_problem(), "Off-design problem")

# Open loops first
s.open_loops()
show_problem(s.assembled_problem(), "Including loop constraints")
# Restore all connectors
s.close_loops()

# Solve problem (solver will eventually open loops)
solver = s.add_driver(NonLinearSolver('solver'))
s.run_drivers()
show_problem(solver.problem, "After resolution")
*** OFF-DESIGN PROBLEM ***
shape = (2, 2)
Unknowns [2]
  z = [0. 0.]
Equations [2]
  y1 == 1 := -0.9
  y2 == 0.5 := -0.18377223398316206

*** INCLUDING LOOP CONSTRAINTS ***
shape = (3, 3)
Unknowns [3]
  z = [0. 0.]
  d1.y2 = 0.0
Equations [3]
  y1 == 1 := -0.9
  y2 == 0.5 := -0.18377223398316206
  d1.y2 == d2.y2 (loop) := -0.31622776601683794

*** AFTER RESOLUTION ***
shape = (3, 3)
Unknowns [3]
  z = [-0.82287566  0.32287566]
  d1.y2 = 0.5
Equations [3]
  y1 == 1 := -1.1102230246251565e-16
  y2 == 0.5 := -1.1102230246251565e-16
  d1.y2 == d2.y2 (loop) := 1.1102230246251565e-16

After loop analysis, it appears an additional unknown, d1.y2, must be computed in order to satisfy condition d1.y2 == d2.y2 imposed by cyclic connections d1.y1 -> d2.y1 and d2.y2 -> d1.y2.

Loop through all elements in a composite system tree

Browsing through composite trees is usually tackled with recursive functions. However, writing such functions can sometimes be tricky.

As of version 0.11.7, a recursive iterator makes it very easy to loop through system or driver trees.

system = SomeComplexSystem('system')

for elem in system.tree():
    # do whatever

This functionality might come handy if one wishes to collect data from all sub-systems, say.

[2]:
from cosapp.base import System

a = System('a')
# First layer
a.add_child(System('aa'))
a.add_child(System('ab'))
# Second layer
a.ab.add_child(System('aba'))
a.ab.add_child(System('abb'))
a.ab.add_child(System('abc'))

print([elem.name for elem in a.tree()])
['aa', 'aba', 'abb', 'abc', 'ab', 'a']

System.tree() accepts optional argument downwards, which determines whether the iterator should yield elements from top to bottom (downwards=True), or from bottom to top (downwards=False - default behaviour).

[3]:
print([elem.name for elem in a.tree(downwards=True)])
['a', 'aa', 'ab', 'aba', 'abb', 'abc']

Note that the iterator follows the execution order of all elements in the tree:

[4]:
print([elem.name for elem in a.tree()])

a.exec_order = ['ab', 'aa']
a.ab.exec_order = ['abc', 'aba', 'abb']
print([elem.name for elem in a.tree()])
['aa', 'aba', 'abb', 'abc', 'ab', 'a']
['abc', 'aba', 'abb', 'ab', 'aa', 'a']

Synchronize sub-system inputs

Inputs/inwards of a system specify the necessary information for local execution of the system. Very often, though, inputs of two or more sub-systems must take the same value, in the context of their parent assembly.

For instance, consider a system PlaneWing with inward length. When gathered in a higher-level system Airplane as left_wing and right_wing, it may be desirable to force both wings to have the same length, by construction.

Such input synchronization is natively handled by CoSApp through the pulling mechanism:

class Airplane(System):
    def setup(self):
        self.add_child(PlaneWing('left_wing', pulling={'length': 'wing_length'}))
        self.add_child(PlaneWing('right_wing', pulling={'length': 'wing_length'}))

Explaination:

Pulling an input/inward creates a parent-to-child downward connector. Thus, the first add_child statement will create new inward wing_length on the fly, as well as a connector wing_length -> left_wing.length. When the right wing sub-system is created, a second connector wing_length -> right_wing.length is generated, using existing inward wing_length.

plane = Airplane('plane')

plane.wing_length = 12.3
plane.run_once()  # will propagate `wing_length` to both wings

Notes

  • Connectors only transfer information when run_once() or run_drivers() is invoked.

  • Any value directly assigned to plane.left_wing.length, say, will be eventually superseded by plane.wing_length at each model execution.

  • Renaming pulled inputs/inwards at parent level is not mandatory. Passing a list of attribute names rather than a dictionary as pulling argument will trigger the same mechanisms, with no name changes. In the example above, it is just more explicit to expose inward wing_length rather than just length in the context of the airplane.

Create collections of sub-systems and/or ports in systems

Consider an arbitrary number of resistors in series, modeled as an assembly of subsystems R1, R2, etc.

In this case, the number of resistors n is provided at setup, and the n resistors are created with a loop of the kind:

def setup(self, n=2):
    for i in range(n)
        self.add_child(Resistor(f"R{i}"))

While this is creation loop is unavoidable, it can be interesting to get an attribute storing all resistors as an iterable collection.

One can pull this trick using method add_property, as illustrated below:

[5]:
from cosapp.base import Port, System


class ElecPort(Port):
    def setup(self):
        self.add_variable("I", 1.0, unit="A", desc="Current")
        self.add_variable("V", 0.0, unit="V", desc="Voltage")


class Resistor(System):
    def setup(self):
        self.add_input(ElecPort, 'elec_in')
        self.add_output(ElecPort, 'elec_out')
        self.add_inward("R", 1e2, unit="ohm", desc="Resistance")

    def compute(self):
        self.elec_out.I = I = self.elec_in.I
        self.elec_out.V = self.elec_in.V - self.R * I


class ResistorSeries(System):
    def setup(self, n=2):
        self.add_property('n', max(int(n), 2))

        # Create a tuple of children
        resistors = tuple(
            self.add_child(Resistor(f"R{i}"))
            for i in range(self.n)
        )
        # Store collection as a read-only property
        self.add_property('resistors', resistors)

        # Connect resistors in series
        for previous, current in zip(resistors, resistors[1:]):
            self.connect(current.elec_in, previous.elec_out)

        # Pull first `elec_in` and last `elec_out`
        self.add_input(ElecPort, 'elec_in')
        self.add_output(ElecPort, 'elec_out')

        self.connect(self.elec_in, resistors[0].elec_in)
        self.connect(self.elec_out, resistors[-1].elec_out)

The same trick can be used to create collections of ports, if needed.

Use of property resistors makes the code clearer, and allows one to loop through the child/port collection without resorting to syntaxes of the kind

s[f"R{i}"] for i in range(s.n)

In the next cell, for example, we show three identical ways of initializing the resistances:

[6]:
s = ResistorSeries('s', n=5)

for res in s.resistors:
    res.R = 1.25e3

# instead of:
for i in range(s.n):
    s[f"R{i}"].R = 1.25e3

# or:
for i in range(s.n):
    setattr(s, f"R{i}.R", 1.25e3)

Note in particular that the first expression, apart from being much simpler, is also more robust to name changes in ResistorSeries, as we do not need to know the naming convention of sub-systems (R0, R1, etc.), nor their total number n.

System properties are immutable by nature, so resistors cannot be redefined. Furthermore, use of tuple garanties that each individual element of the collection is also immutable, which is exactly what we wish to achieve here.

[7]:
import logging

# Check that resistors cannot be individually reassigned:
try:
    s.resistors[0] = None

except Exception as error:
    logging.error(f"Caught {type(error).__name__}: {error}")

ERROR:root:Caught TypeError: 'tuple' object does not support item assignment

Write alternative ways to create a system

Consider a system class containing a 2D numpy array, which we would like to define from local data or from a file. In the next cell, we show a first, awkward implementation, where all parameters are provided to the setup method, and different courses of action are decided, depending on data type.

[8]:
from cosapp.base import System
import numpy

class MySystem(System):
    def setup(self, n: int, **options):
        self.add_property('n', int(n))

        self.add_inward('x', numpy.zeros((self.n, 3)))

        # Get data from array
        try:
            x = options['x']
        except KeyError:
            pass
        else:
            if numpy.shape(x) == self.x.shape:
                self.x = numpy.copy(x)
        # Get data from file
        try:
            filename = options['filename']
        except KeyError:
            pass
        else:
            x = numpy.load(filename)
            if x.shape == self.x.shape:
                self.x = x
        # And so on...
[9]:
s = MySystem('s', n=2, x=[[0, 0.1, 0.2], [-0.5, 0.9, 0.4]])

print(f"s.x = \n{s.x!s}")
s.x =
[[ 0.   0.1  0.2]
 [-0.5  0.9  0.4]]

The complexity of parsing through options at setup can be simplified using ad-hoc functions referred to as factories, declared as class methods by decorator @classmethod. Factories are special class methods whose job is to create and return a new class instance (an object) from a specific set of arguments.

Class methods are invoked with syntax ClassName.method_name(...).

In the next cell, we implement factory methods from_data and from_file, each with its own signature. Note that the first argument of a class method is always cls, denoting the class itself, as opposed to self, used in traditional object-bound methods.

[10]:
from __future__ import annotations
from cosapp.base import System
import numpy


class MySystem(System):
    def setup(self, n: int):
        self.add_property('n', int(n))
        self.add_inward('x', numpy.zeros((self.n, 3)))

    @classmethod
    def from_data(cls, name: str, data: numpy.ndarray) -> MySystem:
        """Factory creating a new system from a numpy array."""
        x = cls.reshape(data)
        s = cls(name, n=x.shape[0])  # newly created system
        s.x = numpy.asarray(data, dtype=float)
        return s

    @classmethod
    def from_file(cls, name: str, filename) -> MySystem:
        """Factory creating a new system from data stored on file."""
        data = numpy.load(filename)
        return cls.from_data(name, data)

    @staticmethod
    def reshape(data) -> numpy.ndarray:
        """Checks that `data` can be reshaped into a (N x 3) array.
        If so, returns the reshaped array; otherwise, raises `ValueError`.
        """
        try:
            return numpy.reshape(data, (-1, 3))
        except ValueError:
            raise ValueError("data must be interpretable as a (N x 3) array")

[11]:
s1 = MySystem('s1', n=4)  # usual way

print(f"{s1.n = }", f"s1.x =\n{s1.x}", sep="\n")

# Create a system from existing data with `from_data`
s2 = MySystem.from_data('s2', [0, 0.1, 0.2, -0.5, 0.9, 0.4])

print("", f"{s2.n = }", f"s2.x =\n{s2.x}", sep="\n")
s1.n = 4
s1.x =
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]

s2.n = 2
s2.x =
[ 0.   0.1  0.2 -0.5  0.9  0.4]
[12]:
# Create a system from existing file with `from_file`:
try:
    s3 = MySystem.from_file('s3', 'data.pickle')

except Exception as error:
    print(
        f"{type(error).__name__}: {error}",
        f"Does not work for s3 since the file does not exist, but you get the idea!",
        sep="\n",
    )
else:
    print(f"{s3.n = }", f"s3.x =\n{s3.x}", sep="\n")

FileNotFoundError: [Errno 2] No such file or directory: 'data.pickle'
Does not work for s3 since the file does not exist, but you get the idea!