CoSApp tutorials on multimode systems
Simple Newton PendulumΒΆ
[1]:
from cosapp.base import System
import numpy as np
class Pendulum(System):
"""Point mass planar pendulum model"""
def setup(self):
self.add_inward('L', 1.00, desc='Length of the rod')
self.add_inward('g', 9.81, desc='Gravitational acceleration')
self.add_outward('acc', 0., desc="Angular acceleration of the mass")
self.add_transient('omega', der='acc', desc="Angular velocity")
self.add_transient('theta', der='omega', desc="Positional angle")
def compute(self):
self.acc = -(self.g / self.L) * np.sin(self.theta)
class NewtonPendulum(System):
def setup(self):
self.add_child(Pendulum('p1'), pulling=['L', 'g'])
self.add_child(Pendulum('p2'), pulling=['L', 'g'])
contact = self.add_event('contact', trigger="p1.theta == p2.theta")
self.add_event('collision', trigger=contact.filter("p1.omega > p2.omega"))
def transition(self):
if self.collision.present:
p1, p2 = self.p1, self.p2
# swap angular velocities
p1.omega, p2.omega = p2.omega, p1.omega
[2]:
from cosapp.drivers import RungeKutta
from cosapp.recorders import DataFrameRecorder
pend = NewtonPendulum("pend")
driver = pend.add_driver(
RungeKutta(order=3, dt=0.01, time_interval=[0, 10])
)
driver.set_scenario(
init = {
"p1.theta" : np.radians(-60),
"p1.omega" : 2,
"p2.theta" : np.radians(22),
"p2.omega" : 0.5,
},
values = {
"L" : 0.5,
"g" : 9.81,
},
# stop = pend.collision.filter('p1.theta < 0'), # try it!
)
# Recorder
driver.add_recorder(
DataFrameRecorder(includes=[
'p?.theta', 'p?.omega',
]),
period = 0.05,
)
# Run simulation & retrieve recorded data
pend.run_drivers()
data = driver.recorder.export_data()
data = data.drop(['Section', 'Status', 'Error code'], axis=1)
[3]:
# Plot results
import plotly.graph_objs as go
traces = [
go.Scatter(
x = data['time'],
y = data['p1.theta'],
mode='lines', name='Pendulum 1', line=dict(color='red'),
),
go.Scatter(
x = data['time'],
y = data['p2.theta'],
mode='lines', name='Pendulum 2', line=dict(color='blue'),
),
]
layout = go.Layout(
title = "Angular positions",
xaxis = dict(title="Time"),
height = 450,
hovermode = "x",
)
go.Figure(data=traces, layout=layout)
Time driver property recorded_events
contains a list of all recorded events. Each element of the list is a named tuple EventRecord(time, events)
, containing the occurrence time, and the list of cascading events at that time, starting by the primitive event.
[4]:
driver.recorded_events
[4]:
[EventRecord(time=0.3129308816366094, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=1.0724999600884106, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=1.8314905848560308, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=2.5902682511923247, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=3.349180640271147, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=4.108586734528638, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=4.868887485626274, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=5.6305622357312926, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=6.394215419574377, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=7.1606375003784, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=7.930878893467985, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=8.706314427786635, events=[<Event pend.contact>, <Event pend.collision>]),
EventRecord(time=9.488609014935344, events=[<Event pend.contact>, <Event pend.collision>])]
Time driver property event_data
contains the same fields as the recorder, but only at event times:
[5]:
driver.event_data.drop(['Section', 'Status', 'Error code'], axis=1)
[5]:
Reference | p1.omega | p1.theta | p2.omega | p2.theta | time | |
---|---|---|---|---|---|---|
0 | t=0.31293088163661 | 4.787848 | 0.188669 | -1.552723 | 0.188669 | 0.312931 |
1 | pend.contact | -1.552723 | 0.188669 | 4.787848 | 0.188669 | 0.312931 |
2 | t=1.0724999600884 | 1.682617 | -0.118677 | -4.831507 | -0.118677 | 1.072500 |
3 | pend.contact | -4.831507 | -0.118677 | 1.682617 | -0.118677 | 1.072500 |
4 | t=1.831490584856 | 4.855791 | 0.045391 | -1.751205 | 0.045391 | 1.831491 |
5 | pend.contact | -1.751205 | 0.045391 | 4.855791 | 0.045391 | 1.831491 |
6 | t=2.5902682511923 | 1.757959 | 0.029126 | -4.858203 | 0.029126 | 2.590268 |
7 | pend.contact | -4.858203 | 0.029126 | 1.757959 | 0.029126 | 2.590268 |
8 | t=3.3491806402711 | 4.838499 | -0.102850 | -1.702831 | -0.102850 | 3.349181 |
9 | pend.contact | -1.702831 | -0.102850 | 4.838499 | -0.102850 | 3.349181 |
10 | t=4.1085867345286 | 1.586258 | 0.173733 | -4.798684 | 0.173733 | 4.108587 |
11 | pend.contact | -4.798684 | 0.173733 | 1.586258 | 0.173733 | 4.108587 |
12 | t=4.8688874856263 | 4.743040 | -0.239635 | -1.409096 | -0.239635 | 4.868887 |
13 | pend.contact | -1.409096 | -0.239635 | 4.743040 | -0.239635 | 4.868887 |
14 | t=5.6305622357313 | 1.172499 | 0.298221 | -4.678179 | 0.298221 | 5.630562 |
15 | pend.contact | -4.678179 | 0.298221 | 1.172499 | 0.298221 | 5.630562 |
16 | t=6.3942154195744 | 4.613127 | -0.346796 | -0.877849 | -0.346796 | 6.394215 |
17 | pend.contact | -0.877849 | -0.346796 | 4.613127 | -0.346796 | 6.394215 |
18 | t=7.1606375003784 | 0.526950 | 0.382039 | -4.559355 | 0.382039 | 7.160638 |
19 | pend.contact | -4.559355 | 0.382039 | 0.526950 | 0.382039 | 7.160638 |
20 | t=7.930878893468 | 4.530440 | -0.399599 | -0.122997 | -0.399599 | 7.930879 |
21 | pend.contact | -0.122997 | -0.399599 | 4.530440 | -0.399599 | 7.930879 |
22 | t=8.7063144277866 | -0.326529 | 0.393557 | -4.540494 | 0.393557 | 8.706314 |
23 | pend.contact | -4.540494 | 0.393557 | -0.326529 | 0.393557 | 8.706314 |
24 | t=9.4886090149353 | 4.599304 | -0.356116 | 0.802758 | -0.356116 | 9.488609 |
25 | pend.contact | 0.802758 | -0.356116 | 4.599304 | -0.356116 | 9.488609 |