CoSApp tutorials on multimode systems
Bouncing Ball¶
[1]:
from cosapp.base import System
import numpy as np
class PointDynamics(System):
"""Point mass dynamics"""
def setup(self):
self.add_inward("mass", 1.0)
self.add_inward("acc_ext", np.zeros(3))
self.add_inward("force_ext", np.zeros(3))
self.add_outward("force", np.zeros(3))
self.add_outward("acc", np.zeros(3))
def compute(self):
self.force = self.force_ext + self.mass * self.acc_ext
self.acc = self.force / self.mass
class PointFriction(System):
"""Point mass ~ v^2 friction model"""
def setup(self):
self.add_inward('v', np.zeros(3), desc="Velocity")
self.add_inward('cf', 0.1, desc="Friction coefficient")
self.add_outward("force", np.zeros(3))
def compute(self):
self.force = (-self.cf * np.linalg.norm(self.v)) * self.v
class PointMassFreeFall(System):
"""Point mass free fall model with friction"""
def setup(self):
self.add_child(PointFriction('friction'), pulling=['cf', 'v'])
self.add_child(PointDynamics('dynamics'), pulling={
'mass': 'mass',
'force': 'force',
'acc_ext': 'g',
'acc': 'a',
})
self.connect(self.friction, self.dynamics, {"force": "force_ext"})
self.add_transient('v', der='a')
self.add_transient('x', der='v')
self.g = np.r_[0, 0, -9.81]
self.exec_order = ['friction', 'dynamics']
class BouncingBall(PointMassFreeFall):
"""Extension of `PointMassFreeFall`, containing a rebound event
and a rebound restitution coefficient.
"""
def setup(self):
super().setup()
self.add_event('rebound', trigger="x[2] <= 0")
self.add_inward('cr', 1.0, desc="Rebound restitution coefficient", limits=(0, 1))
self.add_outward_modevar("n_rebounds", init=0, dtype=int)
def transition(self):
if self.rebound.present:
self.n_rebounds += 1
v = self.v
if abs(v[2]) < 1e-6:
v[2] = 0
else:
v[2] *= -self.cr
[2]:
from cosapp.drivers import RungeKutta
from cosapp.recorders import DataFrameRecorder
ball = BouncingBall('ball')
driver = ball.add_driver(RungeKutta(order=3, dt=0.01))
driver.time_interval = (0, 20)
# Add a recorder to capture time evolution in a dataframe
driver.add_recorder(
DataFrameRecorder(includes=['x', 'v', 'a', 'norm(v)']),
period=0.05,
)
# Define a simulation scenario
driver.set_scenario(
init = {
'x': np.zeros(3),
'v': np.array([8, 0, 9.5]),
},
values = {
'mass': 1.5,
'cf': 0.20,
'cr': 0.98,
},
stop = ball.rebound.filter("norm(v) < 1"),
)
ball.run_drivers()
# Retrieve recorded data
data = driver.recorder.export_data()
data = data.drop(['Section', 'Status', 'Error code'], axis=1)
time = np.asarray(data['time'])
traj = np.asarray(data['x'].tolist())
[3]:
# Plot results
import plotly.graph_objs as go
traces = [
go.Scatter(
x = traj[:, 0],
y = traj[:, 2],
mode = 'lines',
name = 'numerical',
line = dict(color='red'),
)
]
layout = go.Layout(
title = "Trajectory",
xaxis = dict(title="x"),
yaxis = dict(
title = "z",
scaleanchor = "x",
scaleratio = 1,
),
hovermode = "x",
)
go.Figure(data=traces, layout=layout)
In the next cell, we filter the recorded data in the vicinity of a rebound. Notice that there are two entries at event time \(t \approx 1.457\), showing the discontinuity in velocity (all other quantities are here continuous).
[4]:
data[data['time'].between(1.4, 1.6)]
[4]:
Reference | a | norm(v) | time | v | x | |
---|---|---|---|---|---|---|
28 | t=1.4 | [-2.1942952103492357, 0.0, -5.1152226632109015] | 6.234339 | 1.400000 | [2.639768837352256, 0.0, -5.6478849580093815] | [6.268834123218377, 0.0, 0.33123965546824685] |
29 | t=1.45 | [-2.164715456683632, 0.0, -4.7677199187252395] | 6.415203 | 1.450000 | [2.530764095587121, 0.0, -5.894918590886763] | [6.3980912779875725, 0.0, 0.0425971797247012] |
30 | t=1.4572051635645 | [-2.159879336261959, 0.0, -4.718474119601482] | 6.440520 | 1.457205 | [2.515184306625333, 0.0, -5.929093248939675] | [6.416269698167304, 0.0, -1.0894063429134349e-15] |
31 | ball.rebound | [-2.123326599905211, 0.0, -14.715252210789187] | 6.331524 | 1.457205 | [2.515184306625333, 0.0, 5.810511383960882] | [6.416269698167304, 0.0, -1.0894063429134349e-15] |
32 | t=1.5 | [-1.860088399632333, 0.0, -13.791159219544971] | 5.740817 | 1.500000 | [2.4300831558767904, 0.0, 5.201122679003914] | [6.522045524499052, 0.0, 0.2354796735062687] |
33 | t=1.55 | [-1.595276350588942, 0.0, -12.896379761130461] | 5.104618 | 1.550000 | [2.3438721801853557, 0.0, 4.534687458339203] | [6.64133929531685, 0.0, 0.4786887292884219] |
34 | t=1.6 | [-1.3680298317851864, 0.0, -12.16570737139476] | 4.520060 | 1.600000 | [2.269930699223769, 0.0, 3.9087542950279386] | [6.7566370670667295, 0.0, 0.6896227139958433] |
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.
[5]:
driver.recorded_events
[5]:
[EventRecord(time=1.4572051635645258, events=[<Event ball.rebound>]),
EventRecord(time=2.5179323285203927, events=[<Event ball.rebound>]),
EventRecord(time=3.3959818563126203, events=[<Event ball.rebound>]),
EventRecord(time=4.160288296667009, events=[<Event ball.rebound>]),
EventRecord(time=4.844135950925513, events=[<Event ball.rebound>]),
EventRecord(time=5.466803595620114, events=[<Event ball.rebound>]),
EventRecord(time=6.040711897170552, events=[<Event ball.rebound>]),
EventRecord(time=6.574449489921965, events=[<Event ball.rebound>]),
EventRecord(time=7.074264201047492, events=[<Event ball.rebound>]),
EventRecord(time=7.544878421625021, events=[<Event ball.rebound>]),
EventRecord(time=7.989970249199138, events=[<Event ball.rebound>]),
EventRecord(time=8.412474526810431, events=[<Event ball.rebound>]),
EventRecord(time=8.814780162881563, events=[<Event ball.rebound>]),
EventRecord(time=9.198864468284205, events=[<Event ball.rebound>]),
EventRecord(time=9.56638754664407, events=[<Event ball.rebound>]),
EventRecord(time=9.918760402315899, events=[<Event ball.rebound>]),
EventRecord(time=10.2571952264516, events=[<Event ball.rebound>]),
EventRecord(time=10.582743266420342, events=[<Event ball.rebound>]),
EventRecord(time=10.896323842616605, events=[<Event ball.rebound>]),
EventRecord(time=11.198746929216306, events=[<Event ball.rebound>]),
EventRecord(time=11.490730966804389, events=[<Event ball.rebound>]),
EventRecord(time=11.772917081967268, events=[<Event ball.rebound>]),
EventRecord(time=12.045880568609277, events=[<Event ball.rebound>]),
EventRecord(time=12.310140255200968, events=[<Event ball.rebound>]),
EventRecord(time=12.566166211079054, events=[<Event ball.rebound>]),
EventRecord(time=12.814386137028285, events=[<Event ball.rebound>]),
EventRecord(time=13.05519071243782, events=[<Event ball.rebound>]),
EventRecord(time=13.288938099032451, events=[<Event ball.rebound>]),
EventRecord(time=13.515957756114215, events=[<Event ball.rebound>]),
EventRecord(time=13.736553690847542, events=[<Event ball.rebound>]),
EventRecord(time=13.951007249193816, events=[<Event ball.rebound>]),
EventRecord(time=14.159579520456148, events=[<Event ball.rebound>]),
EventRecord(time=14.362513417912503, events=[<Event ball.rebound>]),
EventRecord(time=14.56003548901974, events=[<Event ball.rebound>, <Event ball.stop>])]
Time driver property event_data
contains the same fields as the recorder, but only at event times:
[6]:
driver.event_data.drop(['Section', 'Status', 'Error code'], axis=1)
[6]:
Reference | a | norm(v) | time | v | x | |
---|---|---|---|---|---|---|
0 | t=1.4572051635645 | [-2.159879336261959, 0.0, -4.718474119601482] | 6.440520 | 1.457205 | [2.515184306625333, 0.0, -5.929093248939675] | [6.416269698167304, 0.0, -1.0894063429134349e-15] |
1 | ball.rebound | [-2.159879336261959, 0.0, -4.718474119601482] | 6.331524 | 1.457205 | [2.515184306625333, 0.0, 5.810511383960882] | [6.416269698167304, 0.0, -1.0894063429134349e-15] |
2 | t=2.5179323285204 | [-1.026829786467798, 0.0, -6.68625007328442] | 4.966027 | 2.517932 | [1.5507816107708083, 0.0, -4.717679606530485] | [8.508373346715413, 0.0, -2.7824964554667986e-15] |
3 | ball.rebound | [-1.026829786467798, 0.0, -6.68625007328442] | 4.876481 | 2.517932 | [1.5507816107708083, 0.0, 4.623326014399876] | [8.508373346715413, 0.0, -2.7824964554667986e-15] |
4 | t=3.3959818563126 | [-0.6356294918622605, 0.0, -7.5593022006804995] | 4.188136 | 3.395982 | [1.1382678785776472, 0.0, -4.030487952730085] | [9.674569363224247, 0.0, -6.45108966246255e-14] |
... | ... | ... | ... | ... | ... | ... |
63 | ball.rebound | [-0.028778865925716744, 0.0, -9.668635608755256] | 1.020215 | 14.159580 | [0.20750341921211027, 0.0, 0.9988900438801579] | [14.53024338310035, 0.0, -1.5371384720630488e-14] |
64 | t=14.362513417913 | [-0.027601026244375052, 0.0, -9.676058041450656] | 1.012755 | 14.362513 | [0.2044006679418881, 0.0, -0.9919133278064299] | [14.572036847736637, 0.0, -2.6281060661048627e... |
65 | ball.rebound | [-0.027601026244375052, 0.0, -9.676058041450656] | 0.993333 | 14.362513 | [0.2044006679418881, 0.0, 0.9720750612503013] | [14.572036847736637, 0.0, -2.6281060661048627e... |
66 | t=14.56003548902 | [-0.026501918384140885, 0.0, -9.682995123964925] | 0.986435 | 14.560035 | [0.20149775240327486, 0.0, -0.9656356454798461] | [14.61212277116154, 0.0, 1.3764624206061282e-16] |
67 | ball.rebound | [-0.026501918384140885, 0.0, -9.682995123964925] | 0.967537 | 14.560035 | [0.20149775240327486, 0.0, 0.9463229325702492] | [14.61212277116154, 0.0, 1.3764624206061282e-16] |
68 rows × 6 columns