Transport System

• 

• (pyinduct.trajectory.SignalGenerator):

  

• :

  

• source code:

import pyinduct.core as cr
import pyinduct.placeholder as ph
import pyinduct.registry as reg
import pyinduct.shapefunctions as sh
import pyinduct.simulation as sim
import pyinduct.trajectory as tr
import pyinduct.visualization as vis
import numpy as np
import pyqtgraph as pg
import matplotlib.pyplot as plt

sys_name = 'transport system'
v = 10
l = 5
T = 5
spat_domain = sim.Domain(bounds=(0, l), num=51)
temp_domain = sim.Domain(bounds=(0, T), num=1e2)

init_x = cr.Function(lambda z: 0)

nodes, init_funcs = sh.cure_interval(sh.LagrangeSecondOrder, spat_domain.bounds, node_count=len(spat_domain))
func_label = 'init_funcs'
reg.register_base(func_label, init_funcs)

u = sim.SimulationInputSum([
    tr.SignalGenerator('square', np.array(temp_domain), frequency=0.3, scale=2, offset=4, phase_shift=1),
    tr.SignalGenerator('gausspulse', np.array(temp_domain), phase_shift=temp_domain[15]),
    tr.SignalGenerator('gausspulse', np.array(temp_domain), phase_shift=temp_domain[25], scale=-4),
    tr.SignalGenerator('gausspulse', np.array(temp_domain), phase_shift=temp_domain[35]),
    tr.SignalGenerator('gausspulse', np.array(temp_domain), phase_shift=temp_domain[60], scale=-2),
])

x = ph.FieldVariable(func_label)
psi = ph.TestFunction(func_label)
weak_form = sim.WeakFormulation(
    [
        ph.IntegralTerm(ph.Product(x.derive(temp_order=1), psi), spat_domain.bounds),
        ph.IntegralTerm(ph.Product(x, psi.derive(1)), spat_domain.bounds, scale=-v),
        ph.ScalarTerm(ph.Product(x(l), psi(l)), scale=v),
        ph.ScalarTerm(ph.Product(ph.Input(u), psi(0)), scale=-v)
    ], name=sys_name)

eval_data = sim.simulate_system(weak_form, init_x, temp_domain, spat_domain)

win0 = pg.plot(np.array(eval_data[0].input_data[0]), u.get_results(eval_data[0].input_data[0]),
               labels=dict(left='u(t)', bottom='t'), pen='b')
win0.showGrid(x=False, y=True, alpha=0.5)
vis.save_2d_pg_plot(win0, 'transport_system')
win1 = vis.PgAnimatedPlot(eval_data, title=eval_data[0].name,
                          save_pics=False, labels=dict(left='x(z,t)', bottom='z'))
win2 = vis.MplSlicePlot(eval_data, spatial_point=0, ylabel="x(0,t)")
plt.show()
pg.QtGui.QApplication.instance().exec_()