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.194295210349236, 0.0, -5.1152226632109015]	6.234339	1.400000	[2.6397688373522565, 0.0, -5.6478849580093815]	[6.268834123218378, 0.0, 0.3312396554682464]
29	t=1.45	[-2.1647154566836324, 0.0, -4.7677199187252395]	6.415203	1.450000	[2.5307640955871213, 0.0, -5.894918590886763]	[6.398091277987573, 0.0, 0.04259717972470076]
30	t=1.4572051635645	[-2.1598793362619593, 0.0, -4.718474119601484]	6.440520	1.457205	[2.515184306625334, 0.0, -5.929093248939673]	[6.416269698167304, 0.0, -2.1510571102112408e-16]
31	ball.rebound	[-2.1233265999052118, 0.0, -14.715252210789183]	6.331524	1.457205	[2.515184306625334, 0.0, 5.81051138396088]	[6.416269698167304, 0.0, -2.1510571102112408e-16]
32	t=1.5	[-1.8600883996323319, 0.0, -13.791159219544964]	5.740817	1.500000	[2.430083155876791, 0.0, 5.201122679003909]	[6.522045524499052, 0.0, 0.2354796735062707]
33	t=1.55	[-1.595276350588941, 0.0, -12.896379761130454]	5.104618	1.550000	[2.343872180185356, 0.0, 4.5346874583391985]	[6.64133929531685, 0.0, 0.47868872928842365]
34	t=1.6	[-1.3680298317851862, 0.0, -12.165707371394758]	4.520060	1.600000	[2.2699306992237696, 0.0, 3.908754295027936]	[6.7566370670667295, 0.0, 0.6896227139958451]

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.4572051635645256, events=[<Event ball.rebound>]),
 EventRecord(time=2.5179323285203927, events=[<Event ball.rebound>]),
 EventRecord(time=3.3959818563126207, 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.040711897170551, events=[<Event ball.rebound>]),
 EventRecord(time=6.574449489921964, events=[<Event ball.rebound>]),
 EventRecord(time=7.074264201047491, events=[<Event ball.rebound>]),
 EventRecord(time=7.544878421625021, events=[<Event ball.rebound>]),
 EventRecord(time=7.989970249199137, 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.566387546644071, events=[<Event ball.rebound>]),
 EventRecord(time=9.9187604023159, events=[<Event ball.rebound>]),
 EventRecord(time=10.2571952264516, events=[<Event ball.rebound>]),
 EventRecord(time=10.58274326642034, events=[<Event ball.rebound>]),
 EventRecord(time=10.896323842616601, events=[<Event ball.rebound>]),
 EventRecord(time=11.1987469292163, events=[<Event ball.rebound>]),
 EventRecord(time=11.490730966804382, events=[<Event ball.rebound>]),
 EventRecord(time=11.77291708196726, events=[<Event ball.rebound>]),
 EventRecord(time=12.045880568609268, events=[<Event ball.rebound>]),
 EventRecord(time=12.31014025520096, events=[<Event ball.rebound>]),
 EventRecord(time=12.566166211079043, events=[<Event ball.rebound>]),
 EventRecord(time=12.814386137028272, events=[<Event ball.rebound>]),
 EventRecord(time=13.055190712437806, events=[<Event ball.rebound>]),
 EventRecord(time=13.288938099032437, events=[<Event ball.rebound>]),
 EventRecord(time=13.515957756114199, events=[<Event ball.rebound>]),
 EventRecord(time=13.736553690847526, events=[<Event ball.rebound>]),
 EventRecord(time=13.951007249193802, events=[<Event ball.rebound>]),
 EventRecord(time=14.159579520456134, events=[<Event ball.rebound>]),
 EventRecord(time=14.362513417912487, events=[<Event ball.rebound>]),
 EventRecord(time=14.560035489019723, 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.1598793362619593, 0.0, -4.718474119601484]	6.440520	1.457205	[2.515184306625334, 0.0, -5.929093248939673]	[6.416269698167304, 0.0, -2.1510571102112408e-16]
1	ball.rebound	[-2.1233265999052118, 0.0, -14.715252210789183]	6.331524	1.457205	[2.515184306625334, 0.0, 5.81051138396088]	[6.416269698167304, 0.0, -2.1510571102112408e-16]
2	t=2.5179323285204	[-1.0268297864677984, 0.0, -6.686250073284417]	4.966027	2.517932	[1.5507816107708086, 0.0, -4.717679606530486]	[8.508373346715413, 0.0, -3.434752482434078e-15]
3	ball.rebound	[-1.0083142767195679, 0.0, -12.816074868228071]	4.876481	2.517932	[1.5507816107708086, 0.0, 4.623326014399876]	[8.508373346715413, 0.0, -3.434752482434078e-15]
4	t=3.3959818563126	[-0.6356294918622606, 0.0, -7.5593022006804995]	4.188136	3.395982	[1.1382678785776474, 0.0, -4.030487952730085]	[9.674569363224247, 0.0, -6.495498583447556e-14]
...	...	...	...	...	...	...
63	ball.rebound	[-0.028226418088618272, 0.0, -9.945877703173165]	1.020215	14.159580	[0.20750341921211013, 0.0, 0.9988900438801666]	[14.530243383100348, 0.0, -1.6014967130217883e...
64	t=14.362513417912	[-0.02760102624437481, 0.0, -9.676058041450657]	1.012755	14.362513	[0.20440066794188802, 0.0, -0.9919133278064214]	[14.572036847736637, 0.0, 8.374377580278036e-16]
65	ball.rebound	[-0.027071711600498465, 0.0, -9.938745840105016]	0.993333	14.362513	[0.20440066794188802, 0.0, 0.972075061250293]	[14.572036847736637, 0.0, 8.374377580278036e-16]
66	t=14.56003548902	[-0.026501918384141097, 0.0, -9.682995123964924]	0.986435	14.560035	[0.20149775240327483, 0.0, -0.9656356454798545]	[14.612122771161538, 0.0, -4.0365524507993134e...
67	ball.rebound	[-0.02599421231491442, 0.0, -9.932080365335649]	0.967537	14.560035	[0.20149775240327483, 0.0, 0.9463229325702573]	[14.612122771161538, 0.0, -4.0365524507993134e...

68 rows × 6 columns