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