diff --git a/pySDC/projects/DAE/problems/DiscontinuousTestDAE.py b/pySDC/projects/DAE/problems/DiscontinuousTestDAE.py new file mode 100644 index 0000000000..8df106a828 --- /dev/null +++ b/pySDC/projects/DAE/problems/DiscontinuousTestDAE.py @@ -0,0 +1,179 @@ +import numpy as np + +from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae + + +class DiscontinuousTestDAE(ptype_dae): + r""" + This class implements a scalar test discontinuous differential-algebraic equation (DAE) similar to [1]_. The event function + is defined as :math:`h(y):= 2y - 100`. Then, the discontinuous DAE model reads: + + - if :math:`h(y) \leq 0`: + + .. math:: + \dfrac{d y(t)}{dt} = z(t), + + .. math:: + y(t)^2 - z(t)^2 - 1 = 0, + + else: + + .. math:: + \dfrac{d y(t)}{dt} = 0, + + .. math:: + y(t)^2 - z(t)^2 - 1 = 0, + + for :math:`t \geq 1`. If :math:`h(y) < 0`, the solution is given by + + .. math:: + (y(t), z(t)) = (cosh(t), sinh(t)), + + and the event takes place at :math:`t^* = 0.5 * arccosh(50) = 4.60507` and event point :math:`(cosh(t^*), sinh(t^*))`. + As soon as :math:`t \geq t^*`, i.e., for :math:`h(y) \geq 0`, the solution is given by the constant value + :math:`(cosh(t^*), sinh(t^*))`. + + Parameters + ---------- + nvars : int + Number of unknowns of the system of DAEs. + newton_tol : float + Tolerance for Newton solver. + + Attributes + ---------- + t_switch_exact: float + Exact time of the event :math:`t^*`. + t_switch: float + Time point of the event found by switch estimation. + nswitches: int + Number of switches found by switch estimation. + + References + ---------- + .. [1] L. Lopez, S. Maset. Numerical event location techniques in discontinuous differential algebraic equations. + Appl. Numer. Math. 178, 98-122 (2022). + """ + + def __init__(self, newton_tol=1e-12): + """Initialization routine""" + nvars = 2 + super().__init__(nvars, newton_tol) + self._makeAttributeAndRegister('nvars', localVars=locals(), readOnly=True) + self._makeAttributeAndRegister('newton_tol', localVars=locals()) + + self.t_switch_exact = np.arccosh(50) + self.t_switch = None + self.nswitches = 0 + + def eval_f(self, u, du, t): + r""" + Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. + + Parameters + ---------- + u : dtype_u + Current values of the numerical solution at time t. + du : dtype_u + Current values of the derivative of the numerical solution at time t. + t : float + Current time of the numerical solution. + + Returns + ------- + f : dtype_f + The right-hand side of f (contains two components). + """ + + y, z = u[0], u[1] + dy = du[0] + + t_switch = np.inf if self.t_switch is None else self.t_switch + + h = 2 * y - 100 + f = self.dtype_f(self.init) + + if h >= 0 or t >= t_switch: + f[:] = ( + dy, + y**2 - z**2 - 1, + ) + else: + f[:] = ( + dy - z, + y**2 - z**2 - 1, + ) + return f + + def u_exact(self, t, **kwargs): + r""" + Routine for the exact solution at time :math:`t \leq 1`. For this problem, the exact + solution is piecewise. + + Parameters + ---------- + t : float + Time of the exact solution. + + Returns + ------- + me : dtype_u + Exact solution. + """ + + assert t >= 1, 'ERROR: u_exact only available for t>=1' + + me = self.dtype_u(self.init) + if t <= self.t_switch_exact: + me[:] = (np.cosh(t), np.sinh(t)) + else: + me[:] = (np.cosh(self.t_switch_exact), np.sinh(self.t_switch_exact)) + return me + + def get_switching_info(self, u, t): + r""" + Provides information about the state function of the problem. A change in sign of the state function + indicates an event. If a sign change is detected, a root can be found within the step according to the + intermediate value theorem. + + The state function for this problem is given by + + .. math:: + h(y(t)) = 2 y(t) - 100. + + Parameters + ---------- + u : dtype_u + Current values of the numerical solution at time :math:`t`. + t : float + Current time of the numerical solution. + + Returns + ------- + switch_detected : bool + Indicates whether a discrete event is found or not. + m_guess : int + The index before the sign changes. + state_function : list + Defines the values of the state function at collocation nodes where it changes the sign. + """ + + switch_detected = False + m_guess = -100 + + for m in range(1, len(u)): + h_prev_node = 2 * u[m - 1][0] - 100 + h_curr_node = 2 * u[m][0] - 100 + if h_prev_node < 0 and h_curr_node >= 0: + switch_detected = True + m_guess = m - 1 + break + + state_function = [2 * u[m][0] - 100 for m in range(len(u))] + return switch_detected, m_guess, state_function + + def count_switches(self): + """ + Setter to update the number of switches if one is found. + """ + self.nswitches += 1 diff --git a/pySDC/projects/DAE/problems/WSCC9BusSystem.py b/pySDC/projects/DAE/problems/WSCC9BusSystem.py new file mode 100644 index 0000000000..ad472dc636 --- /dev/null +++ b/pySDC/projects/DAE/problems/WSCC9BusSystem.py @@ -0,0 +1,1225 @@ +import numpy as np +from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae +from pySDC.implementations.datatype_classes.mesh import mesh +from pySDC.core.Errors import ParameterError + + +def WSCC9Bus(): + """ + Returns the Ybus for the power system. + + Returns + ------- + ppc_res : dict + The data with buses, branches, generators (with power flow result) and the Ybus to define the power system. + """ + ppc_res = {} + + ppc_res['baseMVA'] = 100 + ppc_res['Ybus'] = get_initial_Ybus() + ppc_res['bus'] = np.array( + [ + [ + 1.0, + 3.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 1.0, + 0.0, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 2.0, + 2.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 9.999999999999998890e-01, + 9.668741126628123794e00, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 3.0, + 2.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 9.999999999999998890e-01, + 4.771073237177319015e00, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 4.0, + 1.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 9.870068523919054426e-01, + -2.406643919519410257e00, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 5.0, + 1.0, + 90.0, + 30.0, + 0.0, + 0.0, + 1.0, + 9.754721770850530715e-01, + -4.017264326707549849e00, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 6.0, + 1.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 1.003375436452800251e00, + 1.925601686828564363e00, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 7.0, + 1.0, + 100.0, + 35.0, + 0.0, + 0.0, + 1.0, + 9.856448817249467975e-01, + 6.215445553889322738e-01, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 8.0, + 1.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 9.961852458090698637e-01, + 3.799120192692319264e00, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + [ + 9.0, + 1.0, + 125.0, + 50.0, + 0.0, + 0.0, + 1.0, + 9.576210404299042578e-01, + -4.349933576561006987e00, + 345.0, + 1.0, + 1.100000000000000089e00, + 9.000000000000000222e-01, + ], + ] + ) + ppc_res['branch'] = np.array( + [ + [ + 1.0, + 4.0, + 0.0, + 5.759999999999999842e-02, + 0.0, + 250.0, + 250.0, + 250.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + 7.195470158922189796e01, + 2.406895777275899206e01, + -7.195470158922189796e01, + -2.075304453873995314e01, + ], + [ + 4.0, + 5.0, + 1.700000000000000122e-02, + 9.199999999999999845e-02, + 1.580000000000000016e-01, + 250.0, + 250.0, + 250.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + 3.072828027973678999e01, + -5.858508226424568033e-01, + -3.055468555805444453e01, + -1.368795049942141873e01, + ], + [ + 5.0, + 6.0, + 3.899999999999999994e-02, + 1.700000000000000122e-01, + 3.579999999999999849e-01, + 150.0, + 150.0, + 150.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + -5.944531444194475966e01, + -1.631204950057851022e01, + 6.089386583276659337e01, + -1.242746953108854591e01, + ], + [ + 3.0, + 6.0, + 0.0, + 5.859999999999999931e-02, + 0.0, + 300.0, + 300.0, + 300.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + 8.499999999999997158e01, + -3.649025534209526800e00, + -8.499999999999997158e01, + 7.890678351196221740e00, + ], + [ + 6.0, + 7.0, + 1.190000000000000085e-02, + 1.008000000000000007e-01, + 2.089999999999999913e-01, + 150.0, + 150.0, + 150.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + 2.410613416723325741e01, + 4.536791179891427994e00, + -2.401064777894146474e01, + -2.440076244077697609e01, + ], + [ + 7.0, + 8.0, + 8.500000000000000611e-03, + 7.199999999999999456e-02, + 1.489999999999999936e-01, + 250.0, + 250.0, + 250.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + -7.598935222105758669e01, + -1.059923755922268285e01, + 7.649556434279409700e01, + 2.562394697223899231e-01, + ], + [ + 8.0, + 2.0, + 0.0, + 6.250000000000000000e-02, + 0.0, + 250.0, + 250.0, + 250.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + -1.629999999999997158e02, + 2.276189879408803574e00, + 1.629999999999997442e02, + 1.446011953112515869e01, + ], + [ + 8.0, + 9.0, + 3.200000000000000067e-02, + 1.610000000000000042e-01, + 3.059999999999999942e-01, + 250.0, + 250.0, + 250.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + 8.650443565720313188e01, + -2.532429349130207452e00, + -8.403988686535042518e01, + -1.428198298779915731e01, + ], + [ + 9.0, + 4.0, + 1.000000000000000021e-02, + 8.500000000000000611e-02, + 1.759999999999999898e-01, + 250.0, + 250.0, + 250.0, + 0.0, + 0.0, + 1.0, + -360.0, + 360.0, + -4.096011313464404680e01, + -3.571801701219811775e01, + 4.122642130948177197e01, + 2.133889536138378062e01, + ], + ] + ) + ppc_res['gen'] = np.array( + [ + [ + 1.0, + 71.0, + 24.0, + 300.0, + -300.0, + 1.0, + 100.0, + 1.0, + 250.0, + 10.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + ], + [ + 2.0, + 163.0, + 14.0, + 300.0, + -300.0, + 1.0, + 100.0, + 1.0, + 300.0, + 10.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + ], + [ + 3.0, + 85.0, + -3.0, + 300.0, + -300.0, + 1.0, + 100.0, + 1.0, + 270.0, + 10.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + ], + ] + ) + + return ppc_res + + +def get_initial_Ybus(): + ybus = np.array( + [ + [0 - 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j], + [0 + 0j, 0 - 16j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 16j, 0 + 0j], + [0 + 0j, 0 + 0j, 0 - 17.06484641638225j, 0 + 0j, 0 + 0j, 0 + 17.06484641638225j, 0 + 0j, 0 + 0j, 0 + 0j], + [ + 0 + 17.36111111111111j, + 0 + 0j, + 0 + 0j, + 3.307378962025306 - 39.30888872611897j, + -1.942191248714727 + 10.51068205186793j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.36518771331058 + 11.60409556313993j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.942191248714727 + 10.51068205186793j, + 3.224200387138842 - 15.84092701422946j, + -1.282009138424115 + 5.588244962361526j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 17.06484641638225j, + 0 + 0j, + -1.282009138424115 + 5.588244962361526j, + 2.437096619314212 - 32.15386180510696j, + -1.155087480890097 + 9.784270426363173j, + 0 + 0j, + 0 + 0j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.155087480890097 + 9.784270426363173j, + 2.772209954136233 - 23.30324902327162j, + -1.617122473246136 + 13.69797859690844j, + 0 + 0j, + ], + [ + 0 + 0j, + 0 + 16j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.617122473246136 + 13.69797859690844j, + 2.804726852537284 - 35.44561313021703j, + -1.187604379291148 + 5.975134533308591j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.36518771331058 + 11.60409556313993j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.187604379291148 + 5.975134533308591j, + 2.552792092601728 - 17.33823009644852j, + ], + ], + dtype=complex, + ) + + return ybus + + +def get_event_Ybus(): + ybus = np.array( + [ + [0 - 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 17.36111111111111j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j], + [0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j], + [0 + 0j, 0 + 0j, 0 - 17.06484641638225j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 0j, 0 + 17.06484641638225j], + [ + 0 + 17.36111111111111j, + 0 + 0j, + 0 + 0j, + 3.307378962025306 - 39.30888872611897j, + -1.36518771331058 + 11.60409556313993j, + -1.942191248714727 + 10.51068205186793j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.36518771331058 + 11.60409556313993j, + 2.552792092601728 - 17.33823009644852j, + 0 + 0j, + -1.187604379291148 + 5.975134533308591j, + 0 + 0j, + 0 + 0j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.942191248714727 + 10.51068205186793j, + 0 + 0j, + 3.224200387138842 - 15.84092701422946j, + 0 + 0j, + 0 + 0j, + -1.282009138424115 + 5.588244962361526j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.187604379291148 + 5.975134533308591j, + 0 + 0j, + 2.804726852537284 - 19.44561313021703j, + -1.617122473246136 + 13.69797859690844j, + 0 + 0j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + 0 + 0j, + -1.617122473246136 + 13.69797859690844j, + 2.772209954136233 - 23.30324902327162j, + -1.155087480890097 + 9.784270426363173j, + ], + [ + 0 + 0j, + 0 + 0j, + 0 + 17.06484641638225j, + 0 + 0j, + 0 + 0j, + -1.282009138424115 + 5.588244962361526j, + 0 + 0j, + -1.155087480890097 + 9.784270426363173j, + 2.437096619314212 - 32.15386180510696j, + ], + ], + dtype=complex, + ) + + return ybus + + +# def get_YBus(ppc): + +# ppci = ext2int(ppc) +# Ybus, yf, yt = makeYbus(ppci['baseMVA'],ppci['bus'],ppci['branch']) + +# return ppc['Ybus']() + + +class WSCC9BusSystem(ptype_dae): + r""" + Example implementing the WSCC 9 Bus system [1]_. For this complex model, the equations can be found in [2]_, and + sub-transient and turbine parameters are taken from [3]_. The network data of the system are taken from MATPOWER [4]_. + + Parameters + ---------- + nvars : int + Number of unknowns of the system of DAEs (not used here, since it is set up inside this class). + newton_tol : float + Tolerance for Newton solver. + + Attributes + ---------- + mpc : dict + Contains the data for the buses, branches, generators, and the Ybus. + m : int + Number of machines used in the network. + n : int + Number of buses used in the network. + baseMVA : float + Base value of the apparent power. + ws : float + Generator synchronous speed in rad per second. + ws_vector : np.1darray + Vector containing ``ws``. + MD : np.2darray + Machine data. + ED : np.2darray + Excitation data. + TD : np.2darray + Turbine data. + bus : np.2darray + Data for the buses. + branch : np.2darray + Data for the branches in the power system. + gen : np.2darray + Data for generators in the system. + Ybus : np.2darray + Ybus. + YBus_line6_8_outage : np.2darray + Contains the data for the line outage in the power system, where line at bus6 is outaged. + psv_max : float + Maximum valve position. + IC1 : list + Contains the :math:`8`-th row of the ``bus`` matrix. + IC2 : list + Contains the :math:`9`-th row of the ``bus`` matrix. + IC3 : list + Generator values divided by ``baseMVA``. + IC4 : list + Generator values divided by ``baseMVA``. + IC5 : list + Loads divided by ``baseMVA``. + IC6 : list + Loads divided by ``baseMVA``. + genP : list + Contains data for one generator from the ``mpc``. + IC : list + Contains all the values of `IC1, IC2, IC3, IC4, IC5, IC6`. + PL : list + Contains the :math:`5`-th column of ``IC``. + QL : list + Contains the :math:`6`-th column of ``IC``. + PG : np.1darray + Contains the :math:`3`-rd column of ``IC``. + QG : np.1darray + Contains the :math:`4`-th column of ``IC``. + TH0 : np.1darray + Initial condition for angle of bus voltage in rad. + V0 : np.1darray + Contains the :math:`1`-st column of ``IC``, initial condition for magnitude of bus voltage in per unit. + VG0 : np.1darray + Initial condition for complex voltage phasor. + THG0 : np.1darray + Initial condition for angle of the bus voltage in rad. + H : np.1darray + Shaft inertia constant in second. + Xd : np.1darray + d-axis reactance in per unit. + Xdp : np.1darray + Transient d-axis reactance in per unit. + Xdpp : np.1darray + Sub-transient d-axis reactance in per unit. + Xq : np.1darray + q-axis reactance in per unit. + Xqp : np.1darray + Transient q-axis reactance in per unit. + Xqpp : np.1darray + Sub-transient q-axis reactance in per unit. + Td0p : np.1darray + d-axis time constant associated with :math:`E_q'` in second. + Td0pp : np.1darray + d-axis time constant associated with :math:`\psi_{1d}` in second. + Tq0p : np.1darray + q-axis time constant associated with :math:`E_d'` in second. + Tq0pp : np.1darray + q-axis time constant associated with :math:`\psi_{2q}` in second. + Rs : np.1darray + Stator resistance in per unit. + Xls : np.1darray + Parameter :math:`X_{\ell s}`. + Dm : np.1darray + Rotor angle in rad. + KA : np.1darray + Amplifier gain. + TA : np.1darray + Amplifier time constant in second. + KE : np.1darray + Separate or self-excited constant. + TE : np.1darray + Parameter :math:`T_E`. + KF : np.1darray + Parameter _math:`K_F`. + TF : np.1darray + Parameter :math:`T_F`. + Ax : np.1darray + Constant :math:`A_x` of the saturation function :math:`S_{E_i}`. + Bx : np.1darray + Constant :math:`B_x` of the saturation function :math:`S_{E_i}`. + TCH : np.1darray + Incremental steam chest time constant in second. + TSV : np.1darray + Steam valve time constant in second. + RD : np.1darray + Speed regulation quantity in Hz/per unit. + MH : float + Factor :math:`\frac{2 H_i}{w_s}`. + QG : np.1darray + Used to compute :math:`I_{phasor}`. + Vphasor : np.1darray + Complex voltage phasor. + Iphasor : np.1darray + Complex current phasor. + E0 : np.1darray + Initial internal voltage of the synchronous generator. + Em : np.1darray + Absolute values of ``E0``. + D0 : np.1darray + Initial condition for rotor angle in rad. + Id0 : np.1darray + Initial condition for d-axis current in per unit. + Iq0 : np.1darray + Initial condition for q-axis current in per unit. + Edp0 : np.1darray + Initial condition for d-axis transient internal voltages in per unit. + Si2q0 : np.1darray + Initial condition for damper winding 2q flux linkages in per unit. + Eqp0 : np.1darray + Initial condition for q-axis transient internal voltages in per unit. + Si1d0 : np.1darray + Initial condition for damper winding 1d flux linkages in per unit. + Efd0 : np.1darray + Initial condition for field winding fd flux linkages in per unit. + TM0 : np.1darray + Initial condition for mechanical input torque in per unit. + VR0 : np.1darray + Initial condition for exciter input in per unit. + RF0 : np.1darray + Initial condition for exciter feedback in per unit. + Vref : np.1darray + Reference voltage input in per unit. + PSV0 : np.1darray + Initial condition for steam valve position in per unit. + PC : np.1darray + Initial condition for control power input in per unit. + alpha : int + Active load parameter. + beta : int + Reactive load parameter. + bb1, aa1 : list of ndarrays + Used to access on specific values of ``TH``. + bb2, aa2 : list of ndarrays + Used to access on specific values of ``TH``. + t_switch : float + Time the event found by detection. + nswitches : int + Number of events found by detection. + + References + ---------- + .. [1] WSCC 9-Bus System - Illinois Center for a Smarter Electric Grid. https://icseg.iti.illinois.edu/wscc-9-bus-system/ + .. [2] P. W. Sauer, M. A. Pai. Power System Dynamics and Analysis. John Wiley & Sons (2008). + .. [3] I. Abdulrahman. MATLAB-Based Programs for Power System Dynamics Analysis. IEEE Open Access Journal of Power and Energy. + Vol. 7, No. 1, pp. 59–69 (2020). + .. [4] R. D. Zimmerman, C. E. Murillo-Sánchez, R. J. Thomas. MATPOWER: Steady-State Operations, Planning, and Analysis Tools + for Power Systems Research and Education. IEEE Transactions on Power Systems. Vol. 26, No. 1, pp. 12–19 (2011). + """ + + dtype_u = mesh + dtype_f = mesh + + def __init__(self, nvars=None, newton_tol=1e-10, m=3, n=9): + """Initialization routine""" + + nvars = 11 * m + 2 * m + 2 * n + # invoke super init, passing number of dofs + super().__init__(nvars, newton_tol) + self._makeAttributeAndRegister('nvars', 'newton_tol', localVars=locals(), readOnly=True) + self._makeAttributeAndRegister('m', 'n', localVars=locals()) + self.mpc = WSCC9Bus() + + self.baseMVA = self.mpc['baseMVA'] + self.ws = 2 * np.pi * 60 + self.ws_vector = self.ws * np.ones(self.m) + + # Machine data (MD) as a 2D NumPy array + self.MD = np.array( + [ + [23.640, 6.4000, 3.0100], # 1 - H + [0.1460, 0.8958, 1.3125], # 2 - Xd + [0.0608, 0.1198, 0.1813], # 3 - Xdp + [0.0489, 0.0881, 0.1133], # 4 - Xdpp + [0.0969, 0.8645, 1.2578], # 5 - Xq + [0.0969, 0.1969, 0.2500], # 6 - Xqp + [0.0396, 0.0887, 0.0833], # 7 - Xqpp + [8.960000000000001, 6.0000, 5.8900], # 8 - Tdop + [0.1150, 0.0337, 0.0420], # 9 - Td0pp + [0.3100, 0.5350, 0.6000], # 10 - Tqop + [0.0330, 0.0780, 0.1875], # 11 - Tq0pp + [0.0041, 0.0026, 0.0035], # 12 - RS + [0.1200, 0.1020, 0.0750], # 13 - Xls + [ + 0.1 * (2 * 23.64) / self.ws, + 0.2 * (2 * 6.4) / self.ws, + 0.3 * (2 * 3.01) / self.ws, + ], # 14 - Dm (ws should be defined) + ] + ) + + # Excitation data (ED) as a 2D NumPy array + self.ED = np.array( + [ + 20.000 * np.ones(self.m), # 1 - KA + 0.2000 * np.ones(self.m), # 2 - TA + 1.0000 * np.ones(self.m), # 3 - KE + 0.3140 * np.ones(self.m), # 4 - TE + 0.0630 * np.ones(self.m), # 5 - KF + 0.3500 * np.ones(self.m), # 6 - TF + 0.0039 * np.ones(self.m), # 7 - Ax + 1.5550 * np.ones(self.m), # 8 - Bx + ] + ) + + # Turbine data (TD) as a 2D NumPy array + self.TD = np.array( + [ + 0.10 * np.ones(self.m), # 1 - TCH + 0.05 * np.ones(self.m), # 2 - TSV + 0.05 * np.ones(self.m), # 3 - RD + ] + ) + + self.bus = self.mpc['bus'] + self.branch = self.mpc['branch'] + self.gen = self.mpc['gen'] + self.YBus = get_initial_Ybus() + + temp_mpc = self.mpc + temp_mpc['branch'] = np.delete( + temp_mpc['branch'], 6, 0 + ) # note that this is correct but not necessary, because we have the event Ybus already + self.YBus_line6_8_outage = get_event_Ybus() + + # excitation limiter vmax + # self.vmax = 2.1 + self.psv_max = 1.0 + + self.IC1 = [row[7] for row in self.bus] # Column 8 in MATLAB is indexed as 7 in Python (0-based index) + self.IC2 = [row[8] for row in self.bus] # Column 9 in MATLAB is indexed as 8 in Python + + n_prev, m_prev = self.n, self.m + self.n = len(self.bus) # Number of rows in 'bus' list; self.n already defined above?! + self.m = len(self.gen) # Number of rows in 'gen' list; self.m already defined above?! + if n_prev != 9 or m_prev != 3: + raise ParameterError("Number of rows in bus or gen not equal to initialised n or m!") + + gen0 = [0] * self.n + for i in range(self.m): + gen0[i] = self.gen[i][1] + self.genP = gen0 + self.IC3 = [val / self.baseMVA for val in self.genP] + + gen0 = [0] * self.n + for i in range(self.m): + gen0[i] = self.gen[i][2] + genQ = gen0 + for i in range(self.n): + genQ[i] += self.bus[i][5] # Column 6 in MATLAB is indexed as 5 in Python + self.IC4 = [val / self.baseMVA for val in genQ] + + self.IC5 = [row[2] for row in self.bus] # Column 3 in MATLAB is indexed as 2 in Python + self.IC5 = [val / self.baseMVA for val in self.IC5] + + self.IC6 = [row[3] for row in self.bus] # Column 4 in MATLAB is indexed as 3 in Python + self.IC6 = [val / self.baseMVA for val in self.IC6] + + self.IC = list(zip(self.IC1, self.IC2, self.IC3, self.IC4, self.IC5, self.IC6)) + + self.PL = [row[4] for row in self.IC] # Column 5 in MATLAB is indexed as 4 in Python + self.QL = [row[5] for row in self.IC] # Column 6 in MATLAB is indexed as 5 in Python + + self.PG = np.array([row[2] for row in self.IC]) # Column 3 in MATLAB is indexed as 2 in Python + self.QG = np.array([row[3] for row in self.IC]) # Column 4 in MATLAB is indexed as 3 in Python + + self.TH0 = np.array([row[1] * np.pi / 180 for row in self.IC]) + self.V0 = np.array([row[0] for row in self.IC]) + self.VG0 = self.V0[: self.m] + self.THG0 = self.TH0[: self.m] + + # Extracting values from the MD array + self.H = self.MD[0, :] + self.Xd = self.MD[1, :] + self.Xdp = self.MD[2, :] + self.Xdpp = self.MD[3, :] + self.Xq = self.MD[4, :] + self.Xqp = self.MD[5, :] + self.Xqpp = self.MD[6, :] + self.Td0p = self.MD[7, :] + self.Td0pp = self.MD[8, :] + self.Tq0p = self.MD[9, :] + self.Tq0pp = self.MD[10, :] + self.Rs = self.MD[11, :] + self.Xls = self.MD[12, :] + self.Dm = self.MD[13, :] + + # Extracting values from the ED array + self.KA = self.ED[0, :] + self.TA = self.ED[1, :] + self.KE = self.ED[2, :] + self.TE = self.ED[3, :] + self.KF = self.ED[4, :] + self.TF = self.ED[5, :] + self.Ax = self.ED[6, :] + self.Bx = self.ED[7, :] + + # Extracting values from the TD array + self.TCH = self.TD[0, :] + self.TSV = self.TD[1, :] + self.RD = self.TD[2, :] + + # Calculate MH + self.MH = 2 * self.H / self.ws + + # Represent QG as complex numbers + self.QG = self.QG.astype(complex) + + # Convert VG0 and THG0 to complex phasors + self.Vphasor = self.VG0 * np.exp(1j * self.THG0) + + # Calculate Iphasor + self.Iphasor = np.conj(np.divide(self.PG[:m] + 1j * self.QG[:m], self.Vphasor)) + + # Calculate E0 + self.E0 = self.Vphasor + (self.Rs + 1j * self.Xq) * self.Iphasor + + # Calculate Em, D0, Id0, and Iq0 + self.Em = np.abs(self.E0) + self.D0 = np.angle(self.E0) + self.Id0 = np.real(self.Iphasor * np.exp(-1j * (self.D0 - np.pi / 2))) + self.Iq0 = np.imag(self.Iphasor * np.exp(-1j * (self.D0 - np.pi / 2))) + + # Calculate Edp0, Si2q0, Eqp0, and Si1d0 + self.Edp0 = (self.Xq - self.Xqp) * self.Iq0 + self.Si2q0 = (self.Xls - self.Xq) * self.Iq0 + self.Eqp0 = self.Rs * self.Iq0 + self.Xdp * self.Id0 + self.V0[: self.m] * np.cos(self.D0 - self.TH0[: self.m]) + self.Si1d0 = self.Eqp0 - (self.Xdp - self.Xls) * self.Id0 + + # Calculate Efd0 and TM0 + self.Efd0 = self.Eqp0 + (self.Xd - self.Xdp) * self.Id0 + self.TM0 = ( + ((self.Xdpp - self.Xls) / (self.Xdp - self.Xls)) * self.Eqp0 * self.Iq0 + + ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls)) * self.Si1d0 * self.Iq0 + + ((self.Xqpp - self.Xls) / (self.Xqp - self.Xls)) * self.Edp0 * self.Id0 + - ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls)) * self.Si2q0 * self.Id0 + + (self.Xqpp - self.Xdpp) * self.Id0 * self.Iq0 + ) + + # Calculate VR0 and RF0 + self.VR0 = (self.KE + self.Ax * np.exp(self.Bx * self.Efd0)) * self.Efd0 + self.RF0 = (self.KF / self.TF) * self.Efd0 + + # Calculate Vref and PSV0 + self.Vref = self.V0[: self.m] + self.VR0 / self.KA + self.PSV0 = self.TM0 + self.PC = self.PSV0 + + self.alpha = 2 + self.beta = 2 + + self.bb1, self.aa1 = np.meshgrid(np.arange(0, self.m), np.arange(0, self.n)) + self.bb1, self.aa1 = self.bb1.astype(int), self.aa1.astype(int) + + # Create matrix grid to eliminate for-loops (load buses) + self.bb2, self.aa2 = np.meshgrid(np.arange(self.m, self.n), np.arange(0, self.n)) + self.bb2, self.aa2 = self.bb2.astype(int), self.aa2.astype(int) + + self.t_switch = None + self.nswitches = 0 + + def eval_f(self, u, du, t): + r""" + Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. + + Parameters + ---------- + u : dtype_u + Current values of the numerical solution at time t. + du : dtype_u + Current values of the derivative of the numerical solution at time t. + t : float + Current time of the numerical solution. + + Returns + ------- + f : dtype_f + The right-hand side of f (contains two components). + """ + + dEqp, dSi1d, dEdp = du[0 : self.m], du[self.m : 2 * self.m], du[2 * self.m : 3 * self.m] + dSi2q, dDelta = du[3 * self.m : 4 * self.m], du[4 * self.m : 5 * self.m] + dw, dEfd, dRF = du[5 * self.m : 6 * self.m], du[6 * self.m : 7 * self.m], du[7 * self.m : 8 * self.m] + dVR, dTM, dPSV = du[8 * self.m : 9 * self.m], du[9 * self.m : 10 * self.m], du[10 * self.m : 11 * self.m] + + Eqp, Si1d, Edp = u[0 : self.m], u[self.m : 2 * self.m], u[2 * self.m : 3 * self.m] + Si2q, Delta = u[3 * self.m : 4 * self.m], u[4 * self.m : 5 * self.m] + w, Efd, RF = u[5 * self.m : 6 * self.m], u[6 * self.m : 7 * self.m], u[7 * self.m : 8 * self.m] + VR, TM, PSV = u[8 * self.m : 9 * self.m], u[9 * self.m : 10 * self.m], u[10 * self.m : 11 * self.m] + + Id, Iq = u[11 * self.m : 11 * self.m + self.m], u[11 * self.m + self.m : 11 * self.m + 2 * self.m] + V = u[11 * self.m + 2 * self.m : 11 * self.m + 2 * self.m + self.n] + TH = u[11 * self.m + 2 * self.m + self.n : 11 * self.m + 2 * self.m + 2 * self.n] + + # line outage disturbance: + if t >= 0.05: + self.YBus = self.YBus_line6_8_outage + + self.Yang = np.angle(self.YBus) + self.Yabs = np.abs(self.YBus) + + COI = np.sum(w * self.MH) / np.sum(self.MH) + + # Voltage-dependent active loads PL2, and voltage-dependent reactive loads QL2 + PL2 = np.array(self.PL) + QL2 = np.array(self.QL) + + V = V.T + + # Vectorized calculations + Vectorized_angle1 = ( + np.array([TH.take(indices) for indices in self.bb1.T]) + - np.array([TH.take(indices) for indices in self.aa1.T]) + - self.Yang[: self.m, : self.n] + ) + Vectorized_mag1 = (V[: self.m] * V[: self.n].reshape(-1, 1)).T * self.Yabs[: self.m, : self.n] + + sum1 = np.sum(Vectorized_mag1 * np.cos(Vectorized_angle1), axis=1) + sum2 = np.sum(Vectorized_mag1 * np.sin(Vectorized_angle1), axis=1) + + VG = V[: self.m] + THG = TH[: self.m] + Angle_diff = Delta - THG + + Vectorized_angle2 = ( + np.array([TH.take(indices) for indices in self.bb2.T]) + - np.array([TH.take(indices) for indices in self.aa2.T]) + - self.Yang[self.m : self.n, : self.n] + ) + Vectorized_mag2 = (V[self.m : self.n] * V[: self.n].reshape(-1, 1)).T * self.Yabs[self.m : self.n, : self.n] + + sum3 = np.sum(Vectorized_mag2 * np.cos(Vectorized_angle2), axis=1) + sum4 = np.sum(Vectorized_mag2 * np.sin(Vectorized_angle2), axis=1) + + # Initialise f + f = self.dtype_f(self.init) + + t_switch = np.inf if self.t_switch is None else self.t_switch + + # Equations as list + eqs = [] + eqs.append( + (1.0 / self.Td0p) + * ( + -Eqp + - (self.Xd - self.Xdp) + * ( + Id + - ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls) ** 2) * (Si1d + (self.Xdp - self.Xls) * Id - Eqp) + ) + + Efd + ) + - dEqp + ) # (1) + eqs.append((1.0 / self.Td0pp) * (-Si1d + Eqp - (self.Xdp - self.Xls) * Id) - dSi1d) # (2) + eqs.append( + (1.0 / self.Tq0p) + * ( + -Edp + + (self.Xq - self.Xqp) + * ( + Iq + - ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls) ** 2) * (Si2q + (self.Xqp - self.Xls) * Iq + Edp) + ) + ) + - dEdp + ) # (3) + eqs.append((1.0 / self.Tq0pp) * (-Si2q - Edp - (self.Xqp - self.Xls) * Iq) - dSi2q) # (4) + eqs.append(w - COI - dDelta) # (5) + eqs.append( + (self.ws / (2.0 * self.H)) + * ( + TM + - ((self.Xdpp - self.Xls) / (self.Xdp - self.Xls)) * Eqp * Iq + - ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls)) * Si1d * Iq + - ((self.Xqpp - self.Xls) / (self.Xqp - self.Xls)) * Edp * Id + + ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls)) * Si2q * Id + - (self.Xqpp - self.Xdpp) * Id * Iq + - self.Dm * (w - self.ws) + ) + - dw + ) # (6) + eqs.append((1.0 / self.TE) * ((-(self.KE + self.Ax * np.exp(self.Bx * Efd))) * Efd + VR) - dEfd) # (7) + eqs.append((1.0 / self.TF) * (-RF + (self.KF / self.TF) * Efd) - dRF) # (8) + eqs.append( + (1.0 / self.TA) + * (-VR + self.KA * RF - ((self.KA * self.KF) / self.TF) * Efd + self.KA * (self.Vref - V[: self.m])) + - dVR + ) # (9) + + # Limitation of valve position Psv with limiter start + if PSV[0] >= self.psv_max or t >= t_switch: + _temp_dPSV_g1 = (1.0 / self.TSV[1]) * ( + -PSV[1] + self.PSV0[1] - (1.0 / self.RD[1]) * (w[1] / self.ws - 1) + ) - dPSV[1] + _temp_dPSV_g2 = (1.0 / self.TSV[2]) * ( + -PSV[2] + self.PSV0[2] - (1.0 / self.RD[2]) * (w[2] / self.ws - 1) + ) - dPSV[2] + eqs.append(np.array([dPSV[0], _temp_dPSV_g1, _temp_dPSV_g2])) + else: + eqs.append((1.0 / self.TSV) * (-PSV + self.PSV0 - (1.0 / self.RD) * (w / self.ws - 1)) - dPSV) + # Limitation of valve position Psv with limiter end + + eqs.append((1.0 / self.TCH) * (-TM + PSV) - dTM) # (10) + eqs.append( + self.Rs * Id + - self.Xqpp * Iq + - ((self.Xqpp - self.Xls) / (self.Xqp - self.Xls)) * Edp + + ((self.Xqp - self.Xqpp) / (self.Xqp - self.Xls)) * Si2q + + VG * np.sin(Angle_diff) + ) # (12) + eqs.append( + self.Rs * Iq + + self.Xdpp * Id + - ((self.Xdpp - self.Xls) / (self.Xdp - self.Xls)) * Eqp + - ((self.Xdp - self.Xdpp) / (self.Xdp - self.Xls)) * Si1d + + VG * np.cos(Angle_diff) + ) # (13) + eqs.append((Id * VG.T * np.sin(Angle_diff) + Iq * VG.T * np.cos(Angle_diff)) - PL2[0 : self.m] - sum1) # (14) + eqs.append((Id * VG.T * np.cos(Angle_diff) - Iq * VG.T * np.sin(Angle_diff)) - QL2[0 : self.m] - sum2) # (15) + eqs.append(-PL2[self.m : self.n] - sum3) # (16) + eqs.append(-QL2[self.m : self.n] - sum4) # (17) + eqs_flatten = [item for sublist in eqs for item in sublist] + + f[:] = eqs_flatten + return f + + def u_exact(self, t): + r""" + Returns the initial conditions at time :math:`t=0`. + + Parameters + ---------- + t : float + Time of the initial conditions. + + Returns + ------- + me : dtype_u + Initial conditions. + """ + assert t == 0, 'ERROR: u_exact only valid for t=0' + + me = self.dtype_u(self.init) + me[0 : self.m] = self.Eqp0 + me[self.m : 2 * self.m] = self.Si1d0 + me[2 * self.m : 3 * self.m] = self.Edp0 + me[3 * self.m : 4 * self.m] = self.Si2q0 + me[4 * self.m : 5 * self.m] = self.D0 + me[5 * self.m : 6 * self.m] = self.ws_vector + me[6 * self.m : 7 * self.m] = self.Efd0 + me[7 * self.m : 8 * self.m] = self.RF0 + me[8 * self.m : 9 * self.m] = self.VR0 + me[9 * self.m : 10 * self.m] = self.TM0 + me[10 * self.m : 11 * self.m] = self.PSV0 + me[11 * self.m : 11 * self.m + self.m] = self.Id0 + me[11 * self.m + self.m : 11 * self.m + 2 * self.m] = self.Iq0 + me[11 * self.m + 2 * self.m : 11 * self.m + 2 * self.m + self.n] = self.V0 + me[11 * self.m + 2 * self.m + self.n : 11 * self.m + 2 * self.m + 2 * self.n] = self.TH0 + return me + + def get_switching_info(self, u, t, du=None): + r""" + Provides information about the state function of the problem. When the state function changes its sign, + typically an event occurs. So the check for an event should be done in the way that the state function + is checked for a sign change. If this is the case, the intermediate value theorem states a root in this + step. + + The state function for this problem is given by + + .. math:: + h(P_{SV,1}(t)) = P_{SV,1}(t) - P_{SV,1, max} + + for :math:`P_{SV,1,max}=1.0`. + + Parameters + ---------- + u : dtype_u + Current values of the numerical solution at time :math:`t`. + t : float + Current time of the numerical solution. + + Returns + ------- + switch_detected : bool + Indicates whether a discrete event is found or not. + m_guess : int + The index before the sign changes. + state_function : list + Defines the values of the state function at collocation nodes where it changes the sign. + """ + + switch_detected = False + m_guess = -100 + for m in range(1, len(u)): + h_prev_node = u[m - 1][10 * self.m] - self.psv_max + h_curr_node = u[m][10 * self.m] - self.psv_max + if h_prev_node < 0 and h_curr_node >= 0: + switch_detected = True + m_guess = m - 1 + break + + state_function = [u[m][10 * self.m] - self.psv_max for m in range(len(u))] + return switch_detected, m_guess, state_function + + def count_switches(self): + """ + Setter to update the number of switches if one is found. + """ + self.nswitches += 1 diff --git a/pySDC/projects/PinTSimE/paper_PSCC2024/README.md b/pySDC/projects/PinTSimE/paper_PSCC2024/README.md new file mode 100644 index 0000000000..326d577136 --- /dev/null +++ b/pySDC/projects/PinTSimE/paper_PSCC2024/README.md @@ -0,0 +1,40 @@ +# Spectral Deferred Corrections with Discontinuity Handling for Dynamic Power System Simulation + +The Python file `paper_plots.py` creates all the plots contained in the publication: + +**Title:** Spectral Deferred Corrections with Discontinuity Handling for Dynamic Power System Simulation + +**Authors:** Junjie Zhang, Lisa Wimmer, Robert Speck, Matthias Bolten, Kyrill Ho, and Andrea Benigni + +**Conference:** [![PSCC 2024](https://pscc2024.fr/)](https://pscc2024.fr/) + +Current status of the submission: ***submitted*** + +## Plots for the discontinuous test DAE +In order to reproduce the plots for the discontinuous test DAE, the following setting is used: The test DAE `DiscontinuousTestDAE` +is simulated over the time domain with `t0=3.0` and `Tend=5.4` for different step sizes `dt_list = [1 / (2 ** m) for m in range(2, 9)]`. The fully implicit SDC-DAE sweeper `fully_implicit_DAE` solves the problem for different number of collocation +nodes `nnodes=[2, 3, 4, 5]` at Radau IIa nodes `quad_type='RADAU-RIGHT` with LU preconditioner `QI='LU'` using `tol_hybr=1e-6` to solve the nonlinear system. +SDC terminates if the maximum number of iterations `maxiter=45` or the residual tolerance `restol=1e-13` is satisfied. +For event detection, if an event is found the step sizes will be adapted using the factor `alpha=0.95`. A found event should satisfy the tolerance `epsilon_SE=1e-10`. + +Then, executing `make_plots_for_test_DAE()` creates the plots, where functions in the script + +- Fig. 1: `plot_functions_over_time()` for `dt_fix=1 / (2 ** 7)`, +- Fig. 2: `plot_state_function_detection()`, +- Fig. 3: `plot_event_time_error_before_restarts()` for `dt_fix=1 / (2 ** 7)` + +are used. Here, the routine contains additional functions `plot_error_norm()` and `plot_event_time_error()` to create further plots +not used in the publication. The interested applicant is referred to also consider these ones. + +## Plots for the WSCC 9-bus test case +To reproduce the plots for the WSCC 9-bus system, enable the function `make_plots_for_WSCC9_test_case()`. It is recomended to execute the script generating the plots for the WSCC 9-bus test case on a cluster, since the execution takes several hours. Use the following setup: The DAE `WSCC9BusSystem` is simulated over the time domain with `t0=0.0` and `Tend=0.7` for different step sizes `dt_list = [1 / (2 ** m) for m in range(5, 11)]`. The fully implicit SDC-DAE sweeper `fully_implicit_DAE` solves the problem for different number of collocation nodes `nnodes=[2, 3, 4, 5]` at Radau IIa nodes `quad_type='RADAU-RIGHT` with LU preconditioner `QI='LU'` using `tol_hybr=1e-10` to solve the nonlinear system. +SDC terminates if the maximum number of iterations `maxiter=50` or the residual tolerance `restol=5e-13` is satisfied. For event detection, if an event is found the step sizes will be adapted using the factor `alpha=0.95`. A found event should satisfy the tolerance `epsilon_SE=1e-10`. + +Then, executing `make_plots_for_WSCC9_test_case()` creates the plots, where functions in the script + +- Fig. 5: `plot_state_function_detection()`, +- Fig. 6: `plot_functions_over_time()` for `dt_fix=1 / (2 ** 8)` + +are used. + +It is important to note that the `SwitchEstimator` is under ongoing development. Using the newest version of `pySDC` could therefore lead to different results, and hence to different plots. For reproducing the plots in the paper, the commit `74794af5d338affd322a9e1879ae0a6d931f262b` should be used. diff --git a/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py b/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py new file mode 100644 index 0000000000..6db702465a --- /dev/null +++ b/pySDC/projects/PinTSimE/paper_PSCC2024/log_event.py @@ -0,0 +1,50 @@ +from pySDC.core.Hooks import hooks + + +class LogEventDiscontinuousTestDAE(hooks): + """ + Logs the data for the discontinuous test DAE problem containing one discrete event. + Note that this logging data is dependent from the problem itself. + """ + + def post_step(self, step, level_number): + super(LogEventDiscontinuousTestDAE, self).post_step(step, level_number) + + L = step.levels[level_number] + + L.sweep.compute_end_point() + + self.add_to_stats( + process=step.status.slot, + time=L.time + L.dt, + level=L.level_index, + iter=0, + sweep=L.status.sweep, + type='state_function', + value=2 * L.uend[0] - 100, + ) + + +class LogEventWSCC9(hooks): + """ + Logs the data for the discontinuous test DAE problem containing one discrete event. + Note that this logging data is dependent from the problem itself. + """ + + def post_step(self, step, level_number): + super(LogEventWSCC9, self).post_step(step, level_number) + + L = step.levels[level_number] + P = L.prob + + L.sweep.compute_end_point() + + self.add_to_stats( + process=step.status.slot, + time=L.time + L.dt, + level=L.level_index, + iter=0, + sweep=L.status.sweep, + type='state_function', + value=L.uend[10 * P.m] - P.psv_max, + ) diff --git a/pySDC/projects/PinTSimE/paper_PSCC2024/paper_plots.py b/pySDC/projects/PinTSimE/paper_PSCC2024/paper_plots.py new file mode 100644 index 0000000000..626e518d8b --- /dev/null +++ b/pySDC/projects/PinTSimE/paper_PSCC2024/paper_plots.py @@ -0,0 +1,623 @@ +from pathlib import Path +import numpy as np +import dill + +from pySDC.core.Errors import ParameterError + +from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE +from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE +from pySDC.projects.DAE.problems.WSCC9BusSystem import WSCC9BusSystem + +from pySDC.projects.PinTSimE.battery_model import generateDescription +from pySDC.projects.PinTSimE.battery_model import controllerRun +from pySDC.helpers.stats_helper import get_sorted +import pySDC.helpers.plot_helper as plt_helper + +from pySDC.projects.PinTSimE.paper_PSCC2024.log_event import LogEventDiscontinuousTestDAE, LogEventWSCC9 +from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostStep +from pySDC.implementations.hooks.log_restarts import LogRestarts + + +def make_plots_for_test_DAE(): # pragma: no cover + """ + Makes the plot for the discontinuous test DAE, i.e., + + - error over time for fixed time step size for different number of collocation nodes in + comparison with use of switch detection and not, + - error norm for different step sizes and different number of collocation nodes in + comparison with use of switch detection and not, additionally with number of restarts + for each case, + - absolute value of state function at end time for different number of collocation nodes + and different step sizes in comparison with use of switch detection and not, + - event error to exact event time for differen step sizes and different number of + collocation nodes, + - plots event time error of all founded events (not necessarily satisfying the tolerances) + and the maximum value of the state function in this time step for different number of + collocation nodes in comparison with use of switch detection and not. + + Thus, this function contains all the parameters used in the paper for this numerical example. + """ + + Path("data").mkdir(parents=True, exist_ok=True) + + problem_class = DiscontinuousTestDAE + prob_class_name = DiscontinuousTestDAE.__name__ + + sweeper = fully_implicit_DAE + nnodes = [2, 3, 4, 5] + quad_type = 'RADAU-RIGHT' + QI = 'LU' + maxiter = 45 + tol_hybr = 1e-6 + restol = 1e-13 + + hook_class = [LogGlobalErrorPostStep, LogEventDiscontinuousTestDAE, LogRestarts] + + problem_params = dict() + problem_params['newton_tol'] = tol_hybr + + use_detection = [False, True] + max_restarts = 200 + epsilon_SE = 1e-10 + alpha = 0.95 + + t0 = 3.0 + Tend = 5.4 + + dt_list = [1 / (2**m) for m in range(2, 9)] + dt_fix = 1 / (2**7) + + recomputed = False + + results_error_over_time = {} + results_error_norm = {} + results_state_function = {} + results_event_error = {} + results_event_error_restarts = {} + + for M in nnodes: + results_error_over_time[M], results_error_norm[M] = {}, {} + results_state_function[M], results_event_error[M] = {}, {} + results_event_error_restarts[M] = {} + + for dt in dt_list: + results_error_over_time[M][dt], results_error_norm[M][dt] = {}, {} + results_state_function[M][dt], results_event_error[M][dt] = {}, {} + results_event_error_restarts[M][dt] = {} + + for use_SE in use_detection: + results_error_over_time[M][dt][use_SE], results_error_norm[M][dt][use_SE] = {}, {} + results_state_function[M][dt][use_SE], results_event_error[M][dt][use_SE] = {}, {} + results_event_error_restarts[M][dt][use_SE] = {} + + description, controller_params = generateDescription( + dt, + problem_class, + sweeper, + M, + quad_type, + QI, + hook_class, + False, + use_SE, + problem_params, + restol, + maxiter, + max_restarts, + epsilon_SE, + alpha, + ) + + stats, t_switch_exact = controllerRun( + description, controller_params, t0, Tend, exact_event_time_avail=True + ) + + err_val = get_sorted(stats, type='e_global_post_step', sortby='time', recomputed=recomputed) + results_error_over_time[M][dt][use_SE] = err_val + + err_norm = max([item[1] for item in err_val]) + results_error_norm[M][dt][use_SE] = err_norm + + h_val = get_sorted(stats, type='state_function', sortby='time', recomputed=recomputed) + h_abs = abs([item[1] for item in h_val][-1]) + results_state_function[M][dt][use_SE]['h_abs'] = h_abs + + if use_SE: + switches = get_sorted(stats, type='switch', sortby='time', recomputed=recomputed) + + t_switch = [item[1] for item in switches][-1] + results_event_error[M][dt][use_SE] = abs(t_switch_exact - t_switch) + + restarts = get_sorted(stats, type='restart', sortby='time', recomputed=None) + sum_restarts = sum([item[1] for item in restarts]) + results_state_function[M][dt][use_SE]['restarts'] = sum_restarts + + switches_all = get_sorted(stats, type='switch_all', sortby='time', recomputed=None) + t_switches_all = [item[1] for item in switches_all] + event_error_all = [abs(t_switch_exact - t_switch) for t_switch in t_switches_all] + results_event_error_restarts[M][dt][use_SE]['event_error_all'] = event_error_all + h_val_all = get_sorted(stats, type='h_all', sortby='time', recomputed=None) + results_event_error_restarts[M][dt][use_SE]['h_max_event'] = [item[1] for item in h_val_all] + + plot_functions_over_time( + results_error_over_time, prob_class_name, r'Global error $|y(t) - y_{ex}(t)|$', 'upper left', dt_fix + ) + plot_error_norm(results_error_norm, prob_class_name) + plot_state_function_detection( + results_state_function, prob_class_name, r'Absolute value of $h$ $|h(y(T))|$', 'upper left' + ) + plot_event_time_error(results_event_error, prob_class_name) + plot_event_time_error_before_restarts(results_event_error_restarts, prob_class_name, dt_fix) + + +def make_plots_for_WSCC9_test_case(cwd='./'): # pragma: no cover + """ + Generates the plots for the WSCC 9-bus test case, i.e., + + - the values of the state function over time for different number of collocation nodes in comparison + with event detection and not, + - the values of the state function at end time for different number of collocation nodes and + different step sizes. + + Thus, this function contains all the parameters used for this numerical example. + + Parameters + ---------- + cwd : str, optional + Current working directory. + """ + + Path("data").mkdir(parents=True, exist_ok=True) + + problem_class = WSCC9BusSystem + prob_class_name = WSCC9BusSystem.__name__ + + sweeper = fully_implicit_DAE + nnodes = [2, 3, 4, 5] + quad_type = 'RADAU-RIGHT' + QI = 'LU' + maxiter = 50 + tol_hybr = 1e-10 + restol = 5e-13 + + hook_class = [LogEventWSCC9, LogRestarts] + + problem_params = dict() + problem_params['newton_tol'] = tol_hybr + + use_detection = [False, True] + max_restarts = 400 + epsilon_SE = 1e-10 + alpha = 0.95 + + t0 = 0.0 + Tend = 0.7 + + dt_list = [1 / (2**m) for m in range(5, 11)] + dt_fix = 1 / (2**8) + + recomputed = False + + results_state_function_over_time = {} + results_state_function_detection = {} + for M in nnodes: + results_state_function_over_time[M], results_state_function_detection[M] = {}, {} + + for dt in dt_list: + results_state_function_over_time[M][dt], results_state_function_detection[M][dt] = {}, {} + + for use_SE in use_detection: + results_state_function_over_time[M][dt][use_SE], results_state_function_detection[M][dt][use_SE] = ( + {}, + {}, + ) + + description, controller_params = generateDescription( + dt, + problem_class, + sweeper, + M, + quad_type, + QI, + hook_class, + False, + use_SE, + problem_params, + restol, + maxiter, + max_restarts, + epsilon_SE, + alpha, + ) + + # ---- either solution is computed or it is loaded from .dat file already created ---- + path = Path('data/{}_M={}_dt={}_useSE={}.dat'.format(prob_class_name, M, dt, use_SE)) + if path.is_file(): + f = open(cwd + 'data/{}_M={}_dt={}_useSE={}.dat'.format(prob_class_name, M, dt, use_SE), 'rb') + stats = dill.load(f) + f.close() + else: + stats, _ = controllerRun(description, controller_params, t0, Tend) + + fname = 'data/{}_M={}_dt={}_useSE={}.dat'.format(prob_class_name, M, dt, use_SE) + f = open(fname, 'wb') + dill.dump(stats, f) + f.close() + + h_val = get_sorted(stats, type='state_function', sortby='time', recomputed=recomputed) + results_state_function_over_time[M][dt][use_SE] = h_val + + h_abs_end = abs([me[1] for me in h_val][-1]) + results_state_function_detection[M][dt][use_SE]['h_abs'] = h_abs_end + + if use_SE: + restarts = get_sorted(stats, type='restart', sortby='time', recomputed=None) + sum_restarts = sum([me[1] for me in restarts]) + results_state_function_detection[M][dt][use_SE]['restarts'] = sum_restarts + + plot_functions_over_time( + results_state_function_over_time, + prob_class_name, + r'Absolute value of $h$ $|h(P_{SV,0}(t))|$', + 'lower left', + dt_fix, + ) + plot_state_function_detection( + results_state_function_detection, prob_class_name, r'Absolute value of $h$ $|h(P_{SV,0}(T))|$', 'upper right' + ) + + +def plot_styling_stuff(prob_class): # pragma: no cover + """ + Implements all the stuff needed for making the plots more pretty. + """ + + colors = { + 2: 'limegreen', + 3: 'firebrick', + 4: 'deepskyblue', + 5: 'purple', + } + + markers = { + 2: 's', + 3: 'o', + 4: '*', + 5: 'd', + } + + if prob_class == 'DiscontinuousTestDAE': + xytext = { + 2: (-15.0, 16.5), + 3: (-2.0, 55), + 4: (-13.0, -27), + 5: (-1.0, -40), + } + elif prob_class == 'WSCC9BusSystem': + xytext = { + 2: (-13.0, 16), + 3: (-13.0, 30), + 4: (-13.0, -17), + 5: (-1.0, -38), + } + else: + raise ParameterError(f"For {prob_class} no dictionary for position of data points is set up!") + + return colors, markers, xytext + + +def plot_functions_over_time( + results_function_over_time, prob_class, y_label, loc_legend, dt_fix=None +): # pragma: no cover + """ + Plots the functions over time for different numbers of collocation nodes in comparison with detection + and not. + + Parameters + ---------- + results_function_over_time : dict + Results of some function over time for different number of coll.nodes. + prob_class : str + Indicates of which problem class results are plotted (used to define the file name). + y_label : str + y-label used for the plot. + loc_legend : str + Location of the legend in the plot. + dt_fix : bool, optional + If it is set to a considered step size, only one plot will generated. + """ + + colors, _, _ = plot_styling_stuff(prob_class) + x0 = 3.5 if prob_class == 'DiscontinuousTestDAE' else 0.5 + + M_key = list(results_function_over_time.keys())[0] + dt_list = [dt_fix] if dt_fix is not None else results_function_over_time[M_key].keys() + for dt in dt_list: + fig, ax = plt_helper.plt.subplots(1, 1, figsize=(7.5, 5)) + for M in results_function_over_time.keys(): + for use_SE in results_function_over_time[M][dt].keys(): + err_val = results_function_over_time[M][dt][use_SE] + t, err = [item[0] for item in err_val], [abs(item[1]) for item in err_val] + + linestyle_detection = 'solid' if not use_SE else 'dashdot' + (line,) = ax.plot(t, err, color=colors[M], linestyle=linestyle_detection) + + if not use_SE: + line.set_label(r'$M={}$'.format(M)) + + if M == 5: # dummy plot for more pretty legend + ax.plot(x0, 0, color='black', linestyle=linestyle_detection, label='Detection: {}'.format(use_SE)) + + ax.tick_params(axis='both', which='major', labelsize=16) + ax.set_ylim(1e-15, 1e1) + ax.set_yscale('log', base=10) + ax.set_xlabel(r'Time $t$', fontsize=16) + ax.set_ylabel(y_label, fontsize=16) + ax.grid(visible=True) + ax.legend(frameon=True, fontsize=12, loc=loc_legend) + ax.minorticks_off() + + if prob_class == 'DiscontinuousTestDAE': + file_name = 'data/test_DAE_error_over_time_dt{}.png'.format(dt) + elif prob_class == 'WSCC9BusSystem': + file_name = 'data/wscc9_state_function_over_time_dt{}.png'.format(dt) + else: + raise ParameterError(f"For {prob_class} no file name is implemented!") + + fig.savefig(file_name, dpi=300, bbox_inches='tight') + plt_helper.plt.close(fig) + + +def plot_error_norm(results_error_norm, prob_class): # pragma: no cover + """ + Plots the error norm for different step sizes and different number of collocation nodes in comparison + with detection and not. + + Parameters + ---------- + results_error_norm : dict + Statistics containing the error norms and sum of restarts for all considered coll. nodes. + prob_class : str + Indicates of which problem class results are plotted (used to define the file name). + """ + + colors, markers, xytext = plot_styling_stuff(prob_class) + + fig, ax = plt_helper.plt.subplots(1, 1, figsize=(7.5, 5)) + for M in results_error_norm.keys(): + dt = list(results_error_norm[M].keys()) + for use_SE in results_error_norm[M][dt[0]].keys(): + err_norm_dt = [results_error_norm[M][k][use_SE] for k in dt] + + linestyle_detection = 'solid' if not use_SE else 'dashdot' + (line,) = ax.loglog( + dt, + err_norm_dt, + color=colors[M], + linestyle=linestyle_detection, + marker=markers[M], + ) + + if not use_SE: + line.set_label(r'$M={}$'.format(M)) + + if M == 5: # dummy plot for more pretty legend + ax.plot(0, 0, color='black', linestyle=linestyle_detection, label='Detection: {}'.format(use_SE)) + + ax.tick_params(axis='both', which='major', labelsize=16) + ax.set_ylim(1e-15, 1e3) + ax.set_xscale('log', base=10) + ax.set_yscale('log', base=10) + ax.set_xlabel(r'Step size $\Delta t$', fontsize=16) + ax.set_ylabel(r'Error norm $||y(t) - \tilde{y}(t)||_\infty$', fontsize=16) + ax.grid(visible=True) + ax.minorticks_off() + ax.legend(frameon=True, fontsize=12, loc='lower right') + + fig.savefig('data/test_DAE_error_norms.png', dpi=300, bbox_inches='tight') + plt_helper.plt.close(fig) + + +def plot_state_function_detection(results_state_function, prob_class, y_label, loc_legend): # pragma: no cover + """ + Plots the absolute value of the state function after the event which denotes how close it is to the zero. + + Parameters + ---------- + results_state_function : dict + Contains the absolute value of the state function for each number of coll. nodes, each step size and + detection and not. + prob_class : str + Indicates of which problem class results are plotted (used to define the file name). + y_label : str + y-label used for the plot. + loc_legend : str + Location of the legend in the plot. + """ + + colors, markers, xytext = plot_styling_stuff(prob_class) + + fig, ax = plt_helper.plt.subplots(1, 1, figsize=(7.5, 5)) + for M in results_state_function.keys(): + dt = list(results_state_function[M].keys()) + for use_SE in results_state_function[M][dt[0]].keys(): + h_abs = [results_state_function[M][k][use_SE]['h_abs'] for k in dt] + + linestyle_detection = 'solid' if not use_SE else 'dashdot' + (line,) = ax.loglog( + dt, + h_abs, + color=colors[M], + linestyle=linestyle_detection, + marker=markers[M], + ) + + if not use_SE: + line.set_label(r'$M={}$'.format(M)) + + if use_SE: + sum_restarts = [results_state_function[M][k][use_SE]['restarts'] for k in dt] + for m in range(len(dt)): + ax.annotate( + sum_restarts[m], + (dt[m], h_abs[m]), + xytext=xytext[M], + textcoords="offset points", + color=colors[M], + fontsize=16, + ) + + if M == 5: # dummy plot for more pretty legend + ax.plot(0, 0, color='black', linestyle=linestyle_detection, label='Detection: {}'.format(use_SE)) + + ax.tick_params(axis='both', which='major', labelsize=16) + ax.set_ylim(1e-17, 1e3) + ax.set_xscale('log', base=10) + ax.set_yscale('log', base=10) + ax.set_xlabel(r'Step size $\Delta t$', fontsize=16) + ax.set_ylabel(y_label, fontsize=16) + ax.grid(visible=True) + ax.minorticks_off() + ax.legend(frameon=True, fontsize=12, loc=loc_legend) + + if prob_class == 'DiscontinuousTestDAE': + file_name = 'data/test_DAE_state_function.png' + elif prob_class == 'WSCC9BusSystem': + file_name = 'data/wscc9_state_function_detection.png' + else: + raise ParameterError(f"For {prob_class} no file name is set up!") + + fig.savefig(file_name, dpi=300, bbox_inches='tight') + plt_helper.plt.close(fig) + + +def plot_event_time_error(results_event_error, prob_class): # pragma: no cover + """ + Plots the error between event time founded by detection and exact event time. + + Parameters + ---------- + results_event_error : dict + Contains event time error for each considered number of coll. nodes, step size and + event detection and not. + prob_class : str + Indicates of which problem class results are plotted (used to define the file name). + """ + + colors, markers, _ = plot_styling_stuff(prob_class) + + fig, ax = plt_helper.plt.subplots(1, 1, figsize=(7.5, 5)) + for M in results_event_error.keys(): + dt = list(results_event_error[M].keys()) + for use_SE in [True]: + event_error = [results_event_error[M][k][use_SE] for k in dt] + + linestyle_detection = 'solid' if not use_SE else 'dashdot' + ax.loglog( + dt, + event_error, + color=colors[M], + linestyle=linestyle_detection, + marker=markers[M], + label=r'$M={}$'.format(M), + ) + + ax.tick_params(axis='both', which='major', labelsize=16) + ax.set_ylim(1e-15, 1e1) + ax.set_xscale('log', base=10) + ax.set_yscale('log', base=10) + ax.set_xlabel(r'Step size $\Delta t$', fontsize=16) + ax.set_ylabel(r'Event time error $|t^*_{ex} - t^*_{SE}|$', fontsize=16) + ax.grid(visible=True) + ax.minorticks_off() + ax.legend(frameon=True, fontsize=12, loc='lower right') + + fig.savefig('data/test_DAE_event_time_error.png', dpi=300, bbox_inches='tight') + plt_helper.plt.close(fig) + + +def plot_event_time_error_before_restarts(results_event_error_restarts, prob_class, dt_fix=None): # pragma: no cover + """ + Plots all events founded by switch estimation, not necessarily satisfying the conditions. + + Parameters + ---------- + results_event_error_restarts : dict + Contains all events for each considered number of coll. nodes, step size and + event detection and not. + prob_class : str + Indicates of which problem class results are plotted (used to define the file name). + dt_fix : float, optional + Step size considered. + """ + + colors, markers, _ = plot_styling_stuff(prob_class) + + M_key = list(results_event_error_restarts.keys())[0] + dt_list = [dt_fix] if dt_fix is not None else results_event_error_restarts[M_key].keys() + for dt in dt_list: + fig, ax = plt_helper.plt.subplots(1, 1, figsize=(7.5, 5)) + h_ax = ax.twinx() + for M in results_event_error_restarts.keys(): + for use_SE in results_event_error_restarts[M][dt].keys(): + if use_SE: + event_error_all = results_event_error_restarts[M][dt][use_SE]['event_error_all'] + + (line,) = ax.semilogy( + np.arange(1, len(event_error_all) + 1), + event_error_all, + color=colors[M], + linestyle='solid', + # marker=markers[M], + ) + + line.set_label(r'$M={}$'.format(M)) + + h_max_event = results_event_error_restarts[M][dt][use_SE]['h_max_event'] + h_ax.semilogy( + np.arange(1, len(h_max_event) + 1), + h_max_event, + color=colors[M], + linestyle='dashdot', + marker=markers[M], + markersize=5, + alpha=0.4, + ) + + if M == 5: # dummy plot for more pretty legend + ax.plot( + 1, event_error_all[0], color='black', linestyle='solid', label=r'$|t^*_{ex} - t^*_{SE}|$' + ) + ax.plot( + 1, + 1e2, + color='black', + linestyle='dashdot', + marker=markers[M], + markersize=5, + alpha=0.4, + label=r'$||h(t)||_\infty$', + ) + + h_ax.tick_params(labelsize=16) + h_ax.set_ylim(1e-11, 1e0) + h_ax.set_yscale('log', base=10) + h_ax.set_ylabel(r'Maximum value of h $||h(t)||_\infty$', fontsize=16) + h_ax.minorticks_off() + + ax.tick_params(axis='both', which='major', labelsize=16) + ax.set_ylim(1e-11, 1e-1) + ax.set_yscale('log', base=10) + ax.set_xlabel(r'Restarted steps $n_{restart}$', fontsize=16) + ax.set_ylabel(r'Event time error $|t^*_{ex} - t^*_{SE}|$', fontsize=16) + ax.grid(visible=True) + ax.minorticks_off() + ax.legend(frameon=True, fontsize=12, loc='upper right') + + fig.savefig('data/test_DAE_event_time_error_restarts_dt{}.png'.format(dt), dpi=300, bbox_inches='tight') + plt_helper.plt.close(fig) + + +if __name__ == "__main__": + make_plots_for_test_DAE() + # make_plots_for_WSCC9_test_case() diff --git a/pySDC/tests/test_projects/test_DAE/test_problems.py b/pySDC/tests/test_projects/test_DAE/test_problems.py index 0b5e59ccd3..f291a70fdc 100644 --- a/pySDC/tests/test_projects/test_DAE/test_problems.py +++ b/pySDC/tests/test_projects/test_DAE/test_problems.py @@ -325,3 +325,556 @@ def test_synchgen_infinite_bus_main(): # check error err = np.linalg.norm(uend - uend_ref, np.inf) assert np.isclose(err, 0.0, atol=1e-4), "Error too large." + + +@pytest.mark.base +def test_DiscontinuousTestDAE_singularity(): + """ + Test if the event occurs at the correct time and proves if the right-hand side has with the correct values at the event. + """ + import numpy as np + from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE + + t_event = np.arccosh(50.0) + disc_test_DAE = DiscontinuousTestDAE() + + # test for t < t^* by setting t^* = t^* - eps + eps = 1e-3 + t_before_event = t_event - eps + u_before_event = disc_test_DAE.u_exact(t_before_event) + du_before_event = (np.sinh(t_before_event), np.cosh(t_before_event)) + f_before_event = disc_test_DAE.eval_f(u_before_event, du_before_event, t_before_event) + + assert np.isclose(f_before_event[0], 0.0) and np.isclose( + f_before_event[1], 0.0 + ), f"ERROR: Right-hand side after event does not match! Expected {(0.0, 0.0)}, got {f_before_event}" + + # test for t <= t^* + u_event = disc_test_DAE.u_exact(t_event) + du_event = (np.sinh(t_event), np.cosh(t_event)) + f_event = disc_test_DAE.eval_f(u_event, du_event, t_event) + + assert np.isclose(f_event[0], 7 * np.sqrt(51.0)) and np.isclose( + f_event[1], 0.0 + ), f"ERROR: Right-hand side at event does not match! Expected {(7 * np.sqrt(51), 0.0)}, got {f_event}" + + # test for t > t^* by setting t^* = t^* + eps + t_after_event = t_event + eps + u_after_event = disc_test_DAE.u_exact(t_after_event) + du_after_event = (np.sinh(t_event), np.cosh(t_event)) + f_after_event = disc_test_DAE.eval_f(u_after_event, du_after_event, t_after_event) + + assert np.isclose(f_after_event[0], 7 * np.sqrt(51.0)) and np.isclose( + f_after_event[1], 0.0 + ), f"ERROR: Right-hand side after event does not match! Expected {(7 * np.sqrt(51), 0.0)}, got {f_after_event}" + + +@pytest.mark.base +@pytest.mark.parametrize('M', [2, 3, 4, 5]) +def test_DiscontinuousTestDAE_SDC(M): + """ + Simulates one SDC run for different number of coll.nodes and compares if the error satisfies an approppriate value. + """ + + from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE + from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + + # large errors are expected since the simulation domain contains the event + err_tol = { + 2: 0.2025, + 3: 0.2308, + 4: 0.2407, + 5: 0.245, + } + + level_params = { + 'restol': 1e-13, + 'dt': 1e-1, + } + + problem_params = { + 'newton_tol': 1e-6, + } + + sweeper_params = { + 'quad_type': 'RADAU-RIGHT', + 'num_nodes': M, + 'QI': 'IE', + } + + step_params = { + 'maxiter': 45, + } + + controller_params = { + 'logger_level': 30, + } + + description = { + 'problem_class': DiscontinuousTestDAE, + 'problem_params': problem_params, + 'sweeper_class': fully_implicit_DAE, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + } + + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + t0 = 4.6 + Tend = 4.7 + + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + uex = P.u_exact(Tend) + + uend, _ = controller.run(u0=uinit, t0=t0, Tend=Tend) + + err = abs(uex[0] - uend[0]) + assert err < err_tol[M], f"ERROR: Error is too large! Expected {err_tol[M]}, got {err}" + + +@pytest.mark.base +@pytest.mark.parametrize('M', [2, 3, 4, 5]) +def test_DiscontinuousTestDAE_SDC_detection(M): + """ + Test for one SDC run with event detection if the found event is close to the exact value and if the global error + can be reduced. + """ + + from pySDC.helpers.stats_helper import get_sorted + from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE + from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + + err_tol = { + 2: 5.3952e-9, + 3: 2.6741e-9, + 4: 1.9163e-8, + 5: 2.4791e-8, + } + + event_err_tol = { + 2: 3.6968e-5, + 3: 1.3496e-8, + 4: 0.02, + 5: 0.0101, + } + + level_params = { + 'restol': 1e-13, + 'dt': 1e-2, + } + + problem_params = { + 'newton_tol': 1e-6, + } + + sweeper_params = { + 'quad_type': 'RADAU-RIGHT', + 'num_nodes': M, + 'QI': 'IE', + } + + step_params = { + 'maxiter': 45, + } + + controller_params = { + 'logger_level': 30, + } + + switch_estimator_params = { + 'tol': 1e-10, + 'alpha': 0.95, + } + + restarting_params = { + 'max_restarts': 200, + 'crash_after_max_restarts': False, + } + + convergence_controllers = { + SwitchEstimator: switch_estimator_params, + BasicRestartingNonMPI: restarting_params, + } + + description = { + 'problem_class': DiscontinuousTestDAE, + 'problem_params': problem_params, + 'sweeper_class': fully_implicit_DAE, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + 'convergence_controllers': convergence_controllers, + } + + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + t0 = 4.6 + Tend = 4.7 + + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + uex = P.u_exact(Tend) + + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + + err = abs(uex[0] - uend[0]) + assert err < err_tol[M], f"ERROR for M={M}: Error is too large! Expected {err_tol[M]}, got {err}" + + switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) + assert len(switches) >= 1, 'ERROR for M={M}: No events found!' + t_switches = [item[1] for item in switches] + t_switch = t_switches[-1] + + t_switch_exact = P.t_switch_exact + event_err = abs(t_switch_exact - t_switch) + assert ( + event_err < event_err_tol[M] + ), f"ERROR for M={M}: Event error is too large! Expected {event_err_tol[M]}, got {event_err}" + + +@pytest.mark.base +def test_WSCC9_evaluation(): + r""" + Test for WSCC9 bus test case. The class is written for components :math:`m = 3`, :math:`n = 9`. + """ + from pySDC.core.Errors import ParameterError + from pySDC.projects.DAE.problems.WSCC9BusSystem import WSCC9BusSystem + + problem_params = { + 'newton_tol': 1e-10, + } + + WSCC9 = WSCC9BusSystem(**problem_params) + m, n = WSCC9.m, WSCC9.n + nvars = 13 * m + 2 * n + + # test if right-hand side of does have the correct length + t0 = 0.0 + u0 = WSCC9.u_exact(t0) + du0 = np.zeros(len(u0)) + + f = WSCC9.eval_f(u0, du0, t0) + + assert len(f) == nvars, 'Shape of f does not match with shape it is supposed to be!' + + # test if ParameterError is raised if m != 3 or n != 9 is set + problem_params.update( + { + 'm': 4, + 'n': 8, + } + ) + with pytest.raises(ParameterError): + WSCC9_test = WSCC9BusSystem(**problem_params) + + +@pytest.mark.base +def test_WSCC9_update_YBus(): + """ + Test if YBus is updated at time 0.05. For this SDC performs one time step. + """ + + from pySDC.projects.DAE.problems.WSCC9BusSystem import WSCC9BusSystem, get_initial_Ybus, get_event_Ybus + from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + + dt = 0.05 + level_params = { + 'restol': 5e-13, + 'dt': dt, + } + + problem_params = { + 'newton_tol': 1e-10, + } + + sweeper_params = { + 'quad_type': 'RADAU-RIGHT', + 'num_nodes': 2, + 'QI': 'LU', + } + + step_params = { + 'maxiter': 1, + } + + controller_params = { + 'logger_level': 30, + } + + description = { + 'problem_class': WSCC9BusSystem, + 'problem_params': problem_params, + 'sweeper_class': fully_implicit_DAE, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + } + + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + t0 = 0.0 + Tend = dt + + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + + YBus_initial = P.YBus + YBus_initial_ref = get_initial_Ybus() + + assert np.allclose(YBus_initial, YBus_initial_ref), 'YBus does not match with the YBus at initialization!' + + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + + YBus_line_outage = P.YBus + YBus_line6_8_outage = get_event_Ybus() + assert np.allclose( + YBus_line_outage, YBus_line6_8_outage + ), 'YBus after line outage does not match with the one it should supposed to!' + + +@pytest.mark.base +def test_WSCC9_SDC_detection(): + """ + Test if state function states a root. + """ + + from pySDC.helpers.stats_helper import get_sorted + from pySDC.projects.DAE.problems.WSCC9BusSystem import WSCC9BusSystem, get_initial_Ybus, get_event_Ybus + from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + + dt = 0.75 + level_params = { + 'restol': 5e-13, + 'dt': dt, + } + + problem_params = { + 'newton_tol': 1e-10, + } + + sweeper_params = { + 'quad_type': 'RADAU-RIGHT', + 'num_nodes': 2, + 'QI': 'LU', + } + + step_params = { + 'maxiter': 8, + } + + controller_params = { + 'logger_level': 30, + } + + switch_estimator_params = { + 'tol': 1e-10, + 'alpha': 0.95, + } + + restarting_params = { + 'max_restarts': 200, + 'crash_after_max_restarts': False, + } + + convergence_controllers = { + SwitchEstimator: switch_estimator_params, + BasicRestartingNonMPI: restarting_params, + } + + description = { + 'problem_class': WSCC9BusSystem, + 'problem_params': problem_params, + 'sweeper_class': fully_implicit_DAE, + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + 'convergence_controllers': convergence_controllers, + } + + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + + t0 = 0.0 + Tend = dt + + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + + switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) + assert len(switches) >= 1, 'ERROR: No events found!' + t_switch = [me[1] for me in switches][0] + assert np.isclose(t_switch, 0.6103290792685618, atol=1e-3), 'Found event does not match a threshold!' + + +# @pytest.mark.base +# def test_WSCC9_SDC_detection(): +# """ +# Test for one SDC run with event detection if the found event is close to the exact value and if the global error +# can be reduced. +# """ + +# from pySDC.helpers.stats_helper import get_sorted +# from pySDC.projects.DAE.problems.WSCC9BusSystem import WSCC9BusSystem +# from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE +# from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +# from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator +# from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI +# from pySDC.implementations.hooks.log_solution import LogSolution + +# ref = { +# 'Eqp': [1.02565963, 0.87077674, 0.75244422], +# 'Si1d': [1.08501824, 0.8659459, 0.59335005], +# 'Edp': [-9.89321652e-26, 3.55754231e-01, 5.66358724e-01], +# 'Si2q': [0.03994074, -0.380399, -0.67227838], +# 'Delta': [-2.13428304, 10.32368025, -1.48474241], +# 'w': [370.54298062, 398.85092866, 368.59989826], +# 'Efd': [1.33144618, 2.11434102, 2.38996818], +# 'RF': [0.22357495, 0.35186554, 0.36373663], +# 'VR': [1.3316767, 2.48163506, 2.97000777], +# 'TM': [0.98658474, 0.63068939, 1.12527586], +# 'PSV': [1.0, 0.52018862, 1.24497292], +# 'Id': [1.03392984e00, -7.55033973e-36, 1.39602103e00], +# 'Iq': [1.80892723e00, 1.15469164e-30, 6.38447393e-01], +# 'V': [ +# 0.97014097, +# 0.94376174, +# 0.86739643, +# 0.9361775, +# 0.88317809, +# 0.92201319, +# 0.83761267, +# 0.85049254, +# 0.85661891, +# ], +# 'TH': [ +# -2.30672821, +# 9.90481234, +# -2.45484121, +# -2.42758466, +# -2.57057159, +# -2.4746599, +# -2.67639373, +# -2.62752952, +# -2.5584944, +# ], +# 't_switch': [0.5937503078440701], +# } + +# level_params = { +# 'restol': 5e-13, +# 'dt': 1 / (2**5), +# } + +# problem_params = { +# 'newton_tol': 1e-10, +# } + +# sweeper_params = { +# 'quad_type': 'RADAU-RIGHT', +# 'num_nodes': 2, +# 'QI': 'LU', +# } + +# step_params = { +# 'maxiter': 50, +# } + +# controller_params = { +# 'logger_level': 30, +# 'hook_class': LogSolution, +# } + +# switch_estimator_params = { +# 'tol': 1e-10, +# 'alpha': 0.95, +# } + +# restarting_params = { +# 'max_restarts': 400, +# 'crash_after_max_restarts': False, +# } + +# convergence_controllers = { +# SwitchEstimator: switch_estimator_params, +# BasicRestartingNonMPI: restarting_params, +# } + +# description = { +# 'problem_class': WSCC9BusSystem, +# 'problem_params': problem_params, +# 'sweeper_class': fully_implicit_DAE, +# 'sweeper_params': sweeper_params, +# 'level_params': level_params, +# 'step_params': step_params, +# 'convergence_controllers': convergence_controllers, +# } + +# controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + +# t0 = 0.0 +# Tend = 0.7 + +# P = controller.MS[0].levels[0].prob +# uinit = P.u_exact(t0) +# m, n = P.m, P.n + +# uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + +# Eqp = np.array([me[1][0:m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# Si1d = np.array([me[1][m : 2 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# Edp = np.array([me[1][2 * m : 3 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# Si2q = np.array([me[1][3 * m : 4 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# Delta = np.array([me[1][4 * m : 5 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# w = np.array([me[1][5 * m : 6 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# Efd = np.array([me[1][6 * m : 7 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# RF = np.array([me[1][7 * m : 8 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# VR = np.array([me[1][8 * m : 9 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# TM = np.array([me[1][9 * m : 10 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# PSV = np.array([me[1][10 * m : 11 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# Id = np.array([me[1][11 * m : 11 * m + m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# Iq = np.array([me[1][11 * m + m : 11 * m + 2 * m] for me in get_sorted(stats, type='u', sortby='time')])[-1, :] +# V = np.array([me[1][11 * m + 2 * m : 11 * m + 2 * m + n] for me in get_sorted(stats, type='u', sortby='time')])[ +# -1, : +# ] +# TH = np.array( +# [me[1][11 * m + 2 * m + n : 11 * m + 2 * m + 2 * n] for me in get_sorted(stats, type='u', sortby='time')] +# )[-1, :] + +# switches = get_sorted(stats, type='switch', sortby='time', recomputed=False) +# assert len(switches) >= 1, "ERROR: No events found!" +# t_switch = np.array([item[1] for item in switches])[-1] + +# num = { +# 'Eqp': Eqp, +# 'Si1d': Si1d, +# 'Edp': Edp, +# 'Si2q': Si2q, +# 'Delta': Delta, +# 'w': w, +# 'Efd': Efd, +# 'RF': RF, +# 'VR': VR, +# 'TM': TM, +# 'PSV': PSV, +# 'Id': Id, +# 'Iq': Iq, +# 'V': V, +# 'TH': TH, +# 't_switch': t_switch, +# } + +# for key in ref.keys(): +# assert ( +# all(np.isclose(ref[key], num[key], atol=1e-4)) == True +# ), "For {}: Values not equal! Expected {}, got {}".format(key, ref[key], num[key])