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