1"""Diode Bridge simulation. See examples manual.
2Test purpose --> compare with reference results.
3"""
4import os
5from siconos.tests_setup import working_dir
6from siconos.kernel import FirstOrderLinearDS, FirstOrderLinearTIR, \
7    ComplementarityConditionNSL, Interaction, NonSmoothDynamicalSystem, EulerMoreauOSI, \
8    TimeDiscretisation, LCP, TimeStepping
9from numpy import empty
10from siconos.kernel import SimpleMatrix, getMatrix
11from numpy.linalg import norm
12
13
14# Common data
15t0 = 0.0
16total_time = 5.0e-3
17time_step = 1.0e-6
18inductance = 1e-2
19capacitance = 1e-6
20resistance = 1e3
21initial_voltage = 10.0
22model_title = "DiodeBridge"
23init_state = [initial_voltage, 0]
24
25A = [[0, -1.0 / capacitance], [1.0 / inductance, 0]]
26
27C = [[0., 0.], [0, 0.], [-1., 0.], [1., 0.]]
28
29D = [[1. / resistance, 1. / resistance, -1., 0.],
30     [1. / resistance, 1. / resistance, 0., -1.],
31     [1., 0., 0., 0.],
32     [0., 1., 0., 0.]]
33
34B = [[0., 0., -1. / capacitance, 1. / capacitance],
35     [0., 0., 0., 0.]]
36
37
38def test_diode_bridge():
39    """Build diode bridge model"""
40    # dynamical system
41    bridge_ds = FirstOrderLinearDS(init_state, A)
42    # interaction
43    diode_bridge_relation = FirstOrderLinearTIR(C, B)
44    diode_bridge_relation.setDPtr(D)
45
46    nslaw = ComplementarityConditionNSL(4)
47    bridge_interaction = Interaction(nslaw, diode_bridge_relation)
48
49    # Model
50    diode_bridge = NonSmoothDynamicalSystem(t0, total_time)
51    diode_bridge.setTitle(model_title)
52    #  add the dynamical system in the non smooth dynamical system
53    diode_bridge.insertDynamicalSystem(bridge_ds)
54
55    #   link the interaction and the dynamical system
56    diode_bridge.link(bridge_interaction, bridge_ds)
57
58    # Simulation
59
60    # (1) OneStepIntegrators
61    theta = 0.5
62    integrator = EulerMoreauOSI(theta)
63    # (2) Time discretisation
64    time_discretisation = TimeDiscretisation(t0, time_step)
65
66    # (3) Non smooth problem
67    non_smooth_problem = LCP()
68
69    # (4) Simulation setup with (1) (2) (3)
70    bridge_simulation = TimeStepping(diode_bridge,time_discretisation,
71                                     integrator, non_smooth_problem)
72
73
74    k = 0
75    h = bridge_simulation.timeStep()
76    # Number of time steps
77    N = int((total_time - t0) / h)
78
79    # Get the values to be plotted
80    # ->saved in a matrix dataPlot
81    data_plot = empty([N, 8])
82
83    x = bridge_ds.x()
84    print("Initial state : ", x)
85    y = bridge_interaction.y(0)
86    print("First y : ", y)
87    lambda_ = bridge_interaction.lambda_(0)
88
89    # For the initial time step:
90    # time
91    data_plot[k, 0] = t0
92
93    #  inductor voltage
94    data_plot[k, 1] = x[0]
95
96    # inductor current
97    data_plot[k, 2] = x[1]
98
99    # diode R1 current
100    data_plot[k, 3] = y[0]
101
102    # diode R1 voltage
103    data_plot[k, 4] = - lambda_[0]
104
105    # diode F2 voltage
106    data_plot[k, 5] = - lambda_[1]
107
108    # diode F1 current
109    data_plot[k, 6] = lambda_[2]
110
111    # resistor current
112    data_plot[k, 7] = y[0] + lambda_[2]
113
114    k += 1
115    while k < N:
116        bridge_simulation.computeOneStep()
117        #non_smooth_problem.display()
118        data_plot[k, 0] = bridge_simulation.nextTime()
119        #  inductor voltage
120        data_plot[k, 1] = x[0]
121        # inductor current
122        data_plot[k, 2] = x[1]
123        # diode R1 current
124        data_plot[k, 3] = y[0]
125        # diode R1 voltage
126        data_plot[k, 4] = - lambda_[0]
127        # diode F2 voltage
128        data_plot[k, 5] = - lambda_[1]
129        # diode F1 current
130        data_plot[k, 6] = lambda_[2]
131        # resistor current
132        data_plot[k, 7] = y[0] + lambda_[2]
133        k += 1
134        bridge_simulation.nextStep()
135
136    #
137    # comparison with the reference file
138    #
139    ref = getMatrix(SimpleMatrix(os.path.join(working_dir,
140                                              "data/diode_bridge.ref")))
141    assert norm(data_plot - ref) < 1e-12
142    return ref, data_plot
143