1 Introduction

Non-linear dynamics has received a lot of attention since the discovery of the chaotic Lorenz attractor [1]. It opened Pandora’s box that led to a series of further discoveries of other phenomena such as, additional kinds of bifurcations [2, 3], multi-stability and its control [4,5,6,7,8], various routes to chaos and its control [9,10,11,12,13], transient chaotic behaviour [14, 15] or the characterisation of non-linear resonance phenomena [16,17,18,19,20,21], to name a few. Investigating a large number of classical low-dimensional equations, the above mentioned phenomena turned out to be universal features of non-linear systems. The corresponding emerging theories still play an important role in the qualitative understanding of many real-life phenomena in a large variety of scientific fields, for instance, in climate dynamics [22], social sciences [23], neurobiology [24], fluid dynamics [25], mechanical engineering [26] or in laser physics [27].

Although the aforementioned studies are important, they are carried out usually on low-dimensional systems by performing investigations only in low-dimensional parameter spaces or in the local flow of the state space. That is, they require relatively low computational resources compared to an up to date personal computer. However, in order to explore the complex bifurcation structure in parameter space with high resolution [28,29,30,31,32], the necessary computational power can increase by orders of magnitude. For instance, even in a two dimensional parameter plane—employing an initial value problem solver (IVP) with a resolution of \(1000 \times 1000\)—the computational requirements are increased by three orders of magnitude compared to conventional 1D bifurcation plots with the same resolution of 1000. Not to mention if other important, “secondary” control parameters are involved or the application of several initial conditions is mandatory (e.g. to investigate multi-stability). The total number of the parameter combinations can easily blow-up to tens or even hundreds of millions; for instance, see our recent paper about control of multi-stability [4].

At first, it might seem impractical to try to solve a two-dimensional problem with high-resolution IVP computations, since many clever techniques exist (e.g. the pseudo-arclength continuation using a boundary value problem solver (BVP) [33]) that can explore the evolution of bifurcation points even in two dimensions fast and easily. Indeed, in this way, valuable information can be obtained about the bifurcation structures [4, 20, 34,35,36,37]. Nevertheless, these techniques need an already found orbit to initiate the computation. Moreover, they are usually not capable to find a new set of co-existing solutions. Thus, the BVP computations are always combined with IVP simulations, see the aformentioned references. In the present paper, we demonstrate that the application of parameter scans with quite high resolution using IVP solvers can be the source for new ideas and discoveries. For this purpose, computations are carried out on two quite different test models, for details see Sect. 2.

2 The history of the choice of the test models

The first test model is the periodically excited Keller–Miksis equation that is a second order ordinary differential equation describing the radial pulsation of a single spherical gas bubble placed in an infinite domain of liquid [38]. During the radial oscillation of the bubble, due to the external forcing, its contraction phase can be so rapid (collapse) that the temperature inside can reach thousands of degrees of Kelvin inducing chemical reactions [39,40,41]. Therefore, this model is extensively used in the field of sonochemistry [42,43,44,45,46,47,48,49,50,51] to estimate the collapse strength and the chemical yield of a single bubble. In one of our previous papers [4], we extended the investigation to dual-frequency driving using two harmonic components in the external excitation. Therefore, the number of control parameters was increased to four: two driving amplitudes and two driving frequencies (for simplicity, the phase shift between the components was assumed to be zero). Our purpose was to investigate the effect of dual-frequency driving on the dynamics and the collapse strength of a bubble. The main strategy was to create high-resolution bi-parametric maps in the parameter plane of the amplitudes at several fixed frequency pairs. However, during the evaluation of the results, due to the high resolution of the parameter space, special features of the bifurcation structure could be observed. They helped to reveal that with a special choice of the frequencies, specific periodic orbits can be smoothly transformed into each other; for instance, a period-2 and a period-3 attractor. This observation inspired us to develop a non-feedback technique to control multi-stability, in which direct selection of the desired attractors is possible. To the best knowledge of the authors, such a technique was not proposed in the literature before. The present study presents the procedure of the discovery of the technique via an extension of our original work [4].

The second test model (adapted from [52]) describes the dynamics of a pressure relief valve that can exhibit impact dynamics. It is a system of three first-order ordinary differential equations. Our main purpose was to test the special features of the numerical GPU code for non-smooth dynamical systems and reproduce some of the results presented in the original paper [52]. For some additional information about the code, the reader is referred to Sect. 3. There is a special type of impact called grazing impact related to the oscillation of the valve body. It means that the valve body approaches the valve seat, makes contact with the valve seat with zero velocity and then moves away from the seat. At a specific parameter set, the sets of initial conditions from which the pressure relief valve exhibit grazing impact are called grazing lines. They have a focal point in the initial condition space, at which an impacting Shil’nikov-like orbit exist. The grazing lines are computed by means of a BVP solver in the paper of Hős and Champneys [52], which was a cumbersome task that needed special care due to the discontinuous trajectories caused by the impact dynamics. According to the personal communication with the authors, the assembly of their MATLAB code took weeks. Comparing their grazing lines with our GPU accelerated IVP solver, the simulation time is reduced from a couple of hours to seconds. Moreover, the high-resolution scan of the initial conditions revealed a second focal point of the grazing lines that had been overlooked before.

It must be stressed, that in both cases, the original objective was to investigate the collapse strength of a single bubble or to reproduce some results corresponding to a pressure relief valve. The aforementioned discoveries are the “side effects” of the computations of high-resolution multi-dimensional parameter/initial condition scans.

3 The GPU accelerated solver: MPGOS

The usage of an IVP solver performing high-resolution parametric scans is sometimes called the “brute force” technique. It is easy to tune up the number of the parameters, their resolution and the number of the initial conditions; however, to write efficient computer code to do the task within a reasonable time is far from obvious. It is especially true in our case, as we intend to employ the high processing power of professional graphics cards (GPUs). It is not trivial how to use their massively parallel hardware architecture and distribute the workload evenly to tens of thousands of parallel threads.

The developed program package (also used here) is called Massively-Parallel-GPU-ODE-Solver (MPGOS) written in C++ and CUDA C software environments and capable to distribute the tasks to multiple GPUs. It supports explicit solvers: the classic Runge–Kutta solver with fixed time-stepping, and the adaptive Runge–Kutta–Cash–Karp method with embedded error estimation of orders 4 and 5. During the simulations of the present study, the adaptive solver is used. Event handling is also incorporated into the program package. It is mandatory to be able to detect the impact in case of the pressure-relief-valve test model. In addition, with specialized user-defined functions, it is possible to manipulate the trajectories by the user after every successful time step or event detection during the GPU computations. In this way, the impact law can be immediately applied upon the detection of an impact and the integration can be continued. Thus, it is not necessary to stop the integration or perform expensive memory transactions to apply the impact law via the CPU. The code is quite efficient, a simulation is approximately about 50 times faster on an Nvidia GTX GeForce Titan Black card (1707 GFLOPS peak performance) than on a four-core Intel Core i7-4790 CPU (115 GFLOPS peak performance) using double precision floating point arithmetic. The parallelisation strategy in the GPU code follows the “per-thread” approach; that is, to each GPU thread, a different instance of the investigated system is associated having different initial conditions or parameter sets. In the case of the CPU code, the different instances of the system were distributed amongst the CPU cores via the OpenMP application programming interface (API). A single CPU core solved a single instance of a system at a time. It must be emphasised that the proposed speed-up is an estimation using the Keller–Miksis equation introduced in Sect. 4; the achievable factor in the reduction of the runtime can highly depend on the investigated ODE system (handling of special events like impact, or the number of the evaluation of transcendental functions or divisions). The detailed description of the code is beyond the scope of the present paper; however, it has to be stressed that such a fast and efficient solver was the key to achieve the aforementioned discoveries. For more details, the interested reader is referred to the official website of the program package: www.gpuode.com or to its GitHub repository [53]. It is free to use under an MIT license and it has a detailed manual [54] with tutorial examples.

4 Mathematical model of the dual-frequency driven single bubble

The first test model is the Keller–Miksis equation [38]) describing the radial pulsation of a single spherical bubble placed in an infinite domain of liquid. The equation reads as

$$\begin{aligned}&\left( 1-\frac{{\dot{R}}}{c_L} \right) R{\ddot{R}} + \left( 1-\frac{{\dot{R}}}{3c_L} \right) \frac{3}{2} {\dot{R}}^2\nonumber \\&\quad = \left( 1+\frac{{\dot{R}}}{c_L} + \frac{R}{c_L}\frac{d}{dt} \right) \frac{\left( p_L - p_{\infty }(t) \right) }{\rho _L}, \end{aligned}$$
(1)

where R(t) is the time dependent bubble radius. The values of the material properties of the employed liquid (water) are \(c_L=1497.3\,{\mathrm {m/s}}\) (sound speed) and \(\rho _L=997.1\,{\mathrm {kg/m^3}}\) (density). According to the general, dual-frequency treatment, the pressure far away from the bubble,

$$\begin{aligned} p_{\infty }(t) = P_{\infty } + P_{A1} \sin (\omega _1 t) + P_{A2} \sin (\omega _2 t + \theta ), \end{aligned}$$
(2)

is the sum of a static ambient pressure, \(P_{\infty }\), and periodic components with pressure amplitudes \(P_{A1}\) and \(P_{A2}\), angular frequencies \(\omega _1\) and \(\omega _2\), and with a phase shift \(\theta \). The connection between the pressures at the bubble interface can be written as

$$\begin{aligned} p_G + p_V = p_L + \frac{2 \sigma }{R} + 4 \mu _L \frac{{\dot{R}}}{R}, \end{aligned}$$
(3)

where the total pressure inside the bubble is the sum of the partial pressures of the non-condensable gas, \(p_G\), and the vapour, \(p_V=3166.8\,{\mathrm {Pa}}\) at ambient temperature of 25 °C. The surface tension is \(\sigma =0.072\,{\mathrm {N/m}}\) and the liquid kinematic viscosity is \(\mu _L=8.902^{-4}\,{\mathrm {Pa\,s}}\). The gas inside the bubble obeys a simple polytropic relationship

$$\begin{aligned} p_G = \left( P_{\infty } - p_V + \frac{2 \sigma }{R_E} \right) \left( \frac{R_E}{R}\right) ^{3 \gamma }, \end{aligned}$$
(4)

where the polytropic exponent \(\gamma =1.4\) (adiabatic behaviour), the equilibrium bubble radius is RE μm and the static pressure is \(P_{\infty }=1\,{\mathrm {bar}}\).

System (1)–(4) is written into a dimensionless form by introducing the dimensionless variables

$$\begin{aligned} \tau&= \frac{\omega _1}{2 \pi } t, \end{aligned}$$
(5)
$$\begin{aligned} y_1&= \frac{R}{R_E}, \end{aligned}$$
(6)
$$\begin{aligned} y_2&= {\dot{R}} \frac{2\pi }{R_E \omega _1}. \end{aligned}$$
(7)

The equations are rearranged in order to minimize the number of its coefficients. The final form is

$$\begin{aligned} {\dot{y}}_1&= y_2, \end{aligned}$$
(8)
$$\begin{aligned} {\dot{y}}_2&= \frac{N_{KM}}{D_{KM}}, \end{aligned}$$
(9)

where

$$\begin{aligned} N_{KM} =&\left( C_0 + C_1 y_2 \right) \left( \frac{1}{y_1} \right) ^{C_{10}} - C_2 \left( 1 + C_9 y_2 \right) \nonumber \\&- \,C_3 \frac{1}{y_1} - C_4 \frac{y_2}{y_1} - \left( 1 - C_9 \frac{y_2}{3} \right) \frac{3}{2} y_2^2 \nonumber \\&-\,\left( C_5 \sin (2 \pi \tau ) + C_6 \sin (2 \pi C_{11} \tau + C_{12}) \right) \left( 1 + C_9 y_2 \right) \nonumber \\&-\,y_1 \left( C_7 \cos (2 \pi \tau ) + C_8 \cos (2 \pi C_{11} \tau + C_{12}) \right) , \end{aligned}$$
(10)

and

$$\begin{aligned} D_{KM} = y_1 - C_9 y_1 y_2 + C_4 C_9. \end{aligned}$$
(11)

For completeness and reproducibility, the coefficients are summarised below

$$\begin{aligned} C_0&= \frac{1}{\rho _L} \left( P_{\infty } - p_V + \frac{2 \sigma }{R_E} \right) \left( \frac{2 \pi }{R_E \omega _1} \right) ^2, \end{aligned}$$
(12)
$$\begin{aligned} C_1&= \frac{1-3\gamma }{\rho _L c_L} \left( P_{\infty } - p_V + \frac{2 \sigma }{R_E} \right) \frac{2 \pi }{R_E \omega _1}, \end{aligned}$$
(13)
$$\begin{aligned} C_2&= \frac{P_{\infty } - p_V}{\rho _L} \left( \frac{2 \pi }{R_E \omega _1} \right) ^2, \end{aligned}$$
(14)
$$\begin{aligned} C_3&= \frac{2 \sigma }{\rho _L R_E} \left( \frac{2 \pi }{R_E \omega _1} \right) ^2, \end{aligned}$$
(15)
$$\begin{aligned} C_4&= \frac{4 \mu _L}{\rho _L R_E^2} \frac{2 \pi }{\omega _1}, \end{aligned}$$
(16)
$$\begin{aligned} C_5&= \frac{P_{A1}}{\rho _L} \left( \frac{2 \pi }{R_E \omega _1} \right) ^2, \end{aligned}$$
(17)
$$\begin{aligned} C_6&= \frac{P_{A2}}{\rho _L} \left( \frac{2 \pi }{R_E \omega _1} \right) ^2, \end{aligned}$$
(18)
$$\begin{aligned} C_7&= R_E \frac{\omega _1 P_{A1}}{\rho _L c_L} \left( \frac{2 \pi }{R_E \omega _1} \right) ^2, \end{aligned}$$
(19)
$$\begin{aligned} C_8&= R_E \frac{\omega _1 P_{A2}}{\rho _L c_L} \left( \frac{2 \pi }{R_E \omega _1} \right) ^2, \end{aligned}$$
(20)
$$\begin{aligned} C_9&= \frac{R_E \omega _1}{2 \pi c_L}, \end{aligned}$$
(21)
$$\begin{aligned} C_{10}&= 3\gamma , \end{aligned}$$
(22)
$$\begin{aligned} C_{11}&= \frac{\omega _2}{\omega _1}, \end{aligned}$$
(23)
$$\begin{aligned} C_{12}&= \theta . \end{aligned}$$
(24)

The angular frequencies \(\omega _1\) and \(\omega _2\) are normalized by the linear, undamped eigenfrequency [55]

$$\begin{aligned} \omega _0 = \sqrt{ \frac{3\gamma (P_{\infty }-p_V)}{\rho _L R_E^2} - \frac{2(3\gamma -1)\sigma }{\rho _L R_E^3} } = 340\,{\mathrm {kHz}} \end{aligned}$$
(25)

of the unexcited system that defines the relative frequencies as

$$\begin{aligned} \omega _{R1}&= \frac{\omega _1}{\omega _0}, \end{aligned}$$
(26)
$$\begin{aligned} \omega _{R2}&= \frac{\omega _2}{\omega _0}. \end{aligned}$$
(27)

4.1 The global Poincaré section

Due to the dual-frequency driving, the external forcing is not purely harmonic. In Eq. (10), the two dimensionless angular frequencies are \(2 \pi \) and \(2 \pi C_{11}\), here \(C_{11} = \omega _2 / \omega _1 = \omega _{R2} / \omega _{R1}\) is the frequency ratio. The corresponding periods are \(T_1=1\) and \(T_2=1/C_{11}=\omega _{R1} / \omega _{R2}\). For simplicity, the relative phase shift between the harmonic components is set to \(\theta = C_{12} = 0\). During the computations, the main control parameters are the pressure amplitudes while the frequency combinations are kept fixed. The ratio of the employed frequency pairs is always rational; thus, the dual-frequency driving is still periodic (quasiperiodic forcing is excluded). This period T, which is the smallest common multiple of \(T_1\) and \(T_2\) can be used as the global Poincaré section of the system. That is, the trajectories are sampled at time instances \(\tau _n=n \cdot T\) (\(n=0,1,2,\ldots \)).

5 The discovery of a non-feedback technique to directly control multi-stability

In order to represent the dynamical properties of a bubble in a four-dimensional parameter space, our strategy is to compute high-resolution bi-parametric plots with the pressure amplitudes \(P_{A1}\) and \(P_{A2}\) as control parameters applying fixed relative frequency pairs (\(\omega _{R1}\), \(\omega _{R2}\)). The pressure amplitudes are varied between 0 and \(5\,{\mathrm {bar}}\) with 501, uniformly distributed values. In order to explore the co-existing attractors, 10 randomly chosen initial conditions are used. In our experience, it was enough to find the most relevant attractors to draw meaningful conclusions. Thus, a single bi-parametric computation consists of approximately 2.5 million initial value problems. In the first part of the investigation, the relative frequencies are selected from the following set of values:

$$\begin{aligned} \frac{1}{10}, \, \frac{1}{5}, \, \frac{1}{3}, \, \frac{1}{2}, \, \frac{1}{1}, \, \frac{2}{1}, \, \frac{3}{1}, \frac{5}{1} \,\, {\text {and}} \,\, \frac{10}{1}. \end{aligned}$$
(28)

Bi-parametric computations are performed at every possible relative frequency combination, meaning a total number of 36 frequency pairs (taking into account the symmetry property of the driving). In order to explore the subharmonic resonance region in more detail, an additional series of simulations were performed with every possible combination of the frequency values

$$\begin{aligned} \frac{2}{1}, \, \frac{3}{1}, \, \frac{4}{1} \cdots \frac{9}{1}. \end{aligned}$$
(29)

This means 22 additional high-resolution bi-parametric plots (taking into account again the symmetry property and the already computed pairs of frequencies during the first computation period). Thus, the overall number of the solved initial value problems is approximately 145 millions.

At each parameter combination, the first 2048 iterations are regarded as transients and discarded. Then the system is integrated further by additional 8192 iterations to achieve convergence of averaged quantities like the Lyapunov exponent or the winding number. One iteration means the integration of the system from 0 to the period of the excitation T, see Sect. 4.1. To avoid code complexity, the numbers of the iterations mentioned above are the same for all instances of the initial value problems being solved, and they turned out to be enough according to our preliminary calculations. Thus, the convergence of the transients and the average quantities are not monitored. Besides the aforementioned averaged quantities, the period, the maximum bubble radius expansion and the subsequent minimum bubble radius (important to calculate the collapse strength of the bubble oscillation) are also stored. Furthermore, 32 points of the Poincaré section of the last 32 iterations are also recorded. From the various quantities, only the period and the points of the Poincaré section are used in the present study.

Fig. 1
figure 1

Periodicity diagram of bi-parametric plots with pressure amplitudes as control parameters at different relative frequency pairs. The colour code represents the highest period (up to period-6) found at a given parameter set. Inside the black domains, there are chaotic solutions or obits having period higher than six. In the case of co-existing attractors, only the highest period is plotted

Full size image

Figure 1 shows four typical bi-parametric periodicity diagrams at different relative frequency combinations. The colour code represents the maximum period up to period-6 found at a given parameter set. Chaotic oscillations or orbits with periodicity higher than six occupy the black regions. In the case of co-existing attractors, only the highest period is plotted. Keep in mind that the axes in the figures represent single frequency driving since one of the pressure amplitudes is zero in these cases. The bifurcation structure in many of such diagrams shows extreme complexity, where it is hard to find a clear regularity in the bifurcation patterns, see e.g. the upper panels of Fig. 1. However, at specific frequency combinations, bridge shaped structures appear connecting periodic segments from the vertical axis to the horizontal axis, or vice-versa. Such bifurcation structure can be clearly seen in the bottom panels of Fig. 1. Consequently, periodic orbits of single frequency driving at different relative frequencies can be transformed into each other via a temporary dual-frequency driving. The bottom-left panel of Fig. 1 is investigated in more detail in the following to give an in-depth description of the phenomenon.

Figure 2 shows a 3D representation of the period-1 orbits (yellow and gray surfaces) corresponding to relative frequencies \(\omega _{R1}=4\) and \(\omega _{R2}=3\), where the second component of the points of the Poincaré section \(\varPi (y_2)\) is presented as a function of the pressure amplitudes \(P_{A1}\) and \(P_{A2}\). Keep in mind that the global Poincaré section is chosen according to the period of the dual-frequency driving T that is different from the period of the individual components \(T_1\) and \(T_2\). For the present frequency combination, \(T=4\), \(T_1=1\) and \(T_2=4/3\approx 1.333\) in terms of the dimensionless time \(\tau \). That is, the simulation defines every orbit as period-1 that repeats itself after every \(\Delta \tau = 4\). Therefore, for single frequency driving using \(\omega _{R1}=4\) (\(T=4 T_1\)), all period-1 and period-4 orbits are treated as period-1 solutions in the dual-frequency simulations. Similarly, in case of single frequency driven system with \(\omega _{R2}=3\) (\(T=3 T_2\)), all period-1 and period-3 orbits are regarded as period-1 solutions if the dual-frequency Poincaré map is applied. For an exhaustive discussion of the “period reduction” described above, the reader is referred to our previous paper [4].

Fig. 2
figure 2

The second component of the Poincaré section \(\varPi (y_2)\) of the period-1 orbits versus the pressure-amplitude parameter plane of the dual-frequency driving. (Color figure online)

Full size image

Let us summarise the colour code in Fig. 2. The red curves represent period-3 orbits using a single frequency Poincaré map if only the second frequency component is active (\(\omega _{R2}=3\), \(P_{A1}=0\)). The green curves represent period-4 orbits of single frequency driving (again using a single frequency Poincaré map) with relative frequency \(\omega _{R1}=4\) (\(P_{A2}=0\)). Finally, the yellow and grey surfaces and both the red and green curves are the second components of the Poincaré section of period-1 orbits corresponding to the dual-frequency driving (as already discussed above). The surfaces are presented with different colours (yellow and grey) only for the better visibility. It can be clearly seen how these surfaces make connections between the period-3 and period-4 orbits related to different relative frequency values. That is, these two kinds of orbits can be transformed into each other via a temporary dual-frequency excitation.

Although the above-described orbits are related to different relative frequency values, a special kind of control of multi-stability can be achieved in this way if the period-3 (red curves) and the period-4 (green curves) attractors have overlapping domains in the frequency-amplitude parameter plane in case of single frequency driving. However, the transformation works well even if such overlapping domains do not exist. Thus, one can still drive the system from one attractor to another regardless of their co-existence. Observe that such a control technique is a non-feedback method, but the direct selection of the desired attractor is nevertheless possible. Up to now, this was possible only by feedback control techniques [6]. A thorough discussion of the advantages and the drawbacks can be found in our already mentioned previous work [4]; however, only for the transformation between period-2 and period-3 orbits. Therefore, the results presented here indicate that the control technique can be generalised for other pairs of periodic orbits.

It must be emphasized that the high-resolution, multi-dimensional parameter scans have played a vital role in the discovery of the new non-feedback control technique. Since not all the bi-parametric plots show even the sign of the transformation possibility (see e.g. the top panels of Fig. 1), it is very likely that investigating only a few frequency combinations or using coarse resolutions for the pressure amplitudes, we might have missed the special bifurcation structure that led to the discovery. Moreover, as the total number of parameter combinations is of the order of a hundred million, the high-performance GPU computing was a prerequisite of this success.

6 Mathematical model of the pressure relief valve exhibiting impact dynamics

The second test case describes the behaviour of a pressure relief valve that can exhibit impact dynamics. The dimensionless governing equations are adopted from [52] and are written as

$$\begin{aligned} {\dot{y}}_1&= y_2, \end{aligned}$$
(30)
$$\begin{aligned} {\dot{y}}_2&= -\kappa y_2 - (y_1+\delta ) + y_3, \end{aligned}$$
(31)
$$\begin{aligned} {\dot{y}}_3&= \beta ( q - y_1 \sqrt{y_3} ), \end{aligned}$$
(32)

where \(y_1\) and \(y_2\) are the displacement and the velocity of the valve body, respectively. The pressure relief valve is attached to a reservoir chamber in which the dimensionless pressure is \(y_3\). The fixed parameters in the system during the computations are as follows: \(\kappa =1.25\) is the damping coefficient, \(\delta =10\) is the precompression parameter, \(\beta =20\) is the compressibility parameter and \(q=0.3\) is the dimensionless flow rate.

In Eqs. (30)–(32), the zero value of the displacement (\(y_1=0\)) means that the valve body is in contact with the seat of the valve. If the velocity of the valve body \(y_2\) has a non-zero, negative value at this point, the following impact law is applied:

$$\begin{aligned} y_1^+&= y_1^- = 0, \end{aligned}$$
(33)
$$\begin{aligned} y_2^+&= -r y_2^-, \end{aligned}$$
(34)
$$\begin{aligned} y_3^+&= y_3^- \end{aligned}$$
(35)

That is, the velocity of the valve body is reversed by the Newtonian coefficient of restitution \(r=0.8\) that approximates the loss of energy of the impact.

7 The discovery of a new focal point of grazing lines

During the oscillation of the valve body of a pressure relief valve, it can exhibit impact dynamics (the valve body is in contact with the valve seat) that can be categorised as follows. The transversal impact has a non-zero velocity during the impact (\(y_2<0\)); that is, it is a “normal” impact. Whereas, the so-called grazing impact occurs when the impact happens with a zero velocity (\(y_2=0\)). In this case, the impact law has no real effect as the valve body only touches the valve seat. Figure 3 shows the \(y_1\) component of two trajectories that exhibit impacts (\(y_1=0\)). The red dots denote the grazing impacts. The simulations are stopped at the next impact. In both cases, the initial conditions for the first two components are \(y_{10}=0\) and \(y_{20}=0.4\). The only difference is in the third initial condition: \(y_{30}=8.66\) and \(y_{30}=8.58\) depicted also in the figure. The employed parameter set is summarized in Sect. 6. The grazing impact can also be labelled (for a specific initial condition) according to how many transversal impacts there were before. Thus, in Fig. 3, the grazing impacts are denoted as \(G^{0}\) (zero transversal impact) and \(G^{2}\) (two transversal impacts).

Fig. 3
figure 3

Time series exhibiting transversal and grazing impacts applying different initial conditions. The grazing impacts are marked by the red dots. (Color figure online)

Full size image

From a theoretical point of view, the generalization of the grazing impacts to the \(y_{20}-y_{30}\) initial condition plane is an interesting problem. The first component of the initial condition is always set to \(y_{10}=0\). In this way, \(G^{(k)}\) denotes a set of points in the \(y_{20}-y_{30}\) initial condition plane, which leads to a grazing impact after k transversal impacts. Throughout this paper, we shall call such a set of points as grazing lines of order k. The first seven grazing lines computed by Hős and Champneys [52] are shown in the bottom-right panel of Fig. 4. Their strategy was to use a BVP solver and to employ the pseudo-arclength continuation technique to follow the path of the curves initiated from the results of an IVP solver. This formalism is quite complex, as for a single BVP, one needs to define sub-BVPs for each of the \(k+1\) segments divided by the impacts. These are coupled via the impact law for the internal connections. At one side of the full-BVP, the grazing condition, while at the other side of the full-BVP, the condition \(y_{1}=0\) has to be prescribed. Furthermore, the time instances of the intermediate transversal impacts need to be tracked properly as well. The main drawback of this approach is that for different values of k, a different set of BVPs has to be set up and solved. These are the main reasons why the total computational time of a single grazing line was as high as several hours (according to personal communications with the authors). In addition, the implementation of the solver took weeks. The main outcome of the results is that the grazing lines are organized as spirals with a single focal point; and at this focal point, a Shil’nikov-like orbit exists with impacts, see again [52].

Fig. 4
figure 4

Grazing lines computed by the GPU accelerated IVP solver (colour-coded panels) and the BVP solver with the pseudo-arclength continuation technique (bottom-right panel, reprinted with permission from Hős and Champneys [52]). The pure white colour represents the zero value of velocity of the valve body of an impact (grazing impact). The pure red colour means the velocity of \(1.5\,{\mathrm {m/s}}\) or higher. Between \(1.5\,{\mathrm {m/s}}\) and \(0\,{\mathrm {m/s}}\), the transition is uniform in the colour code. (Color figure online)

Full size image

Another way to compute the grazing lines is to take an IVP solver (like our GPU accelerated solver), solve the system forward in time, stop the integration after \(k+1\) impacts and register the velocity of the endpoint \(y_{2E}\). If this velocity is zero, the corresponding initial condition lies on a grazing line of order k denoted as \(G^{(k)}\). With a fine resolution of the set of initial conditions in the \(y_{20}-y_{30}\) plane, the grazing lines can be drawn easily by creating a contour plot of the \(y_{2E}\) value. Theoretically, the zero iso-lines shall represent the corresponding grazing line.

The \(G^{(1)}\) curve computed with our GPU-ODE solver is presented in the top-left panel of Fig. 4 via a white-red colour-coded plot. Here the integrations are stopped at the second impact (k + 1 = 2). The resolution of the initial conditions is \(1024 \times 1024\) and the total computation time is merely \(4\,{\mathrm {s}}\). The pure white colour represents the zero value of \(y_{2E}\). The pure red colour means \(y_{2E}>1.5\,{\mathrm {m/s}}\). Between \(1.5>y_{2E}>0\), the transition is uniform in the colour code. Interestingly, the zero values always lie at a discontinuity, see the jump in the colour code labelled by \(G^{(1)}\) in the top-left panel of Fig. 4. Accordingly, the grazing lines can be easily identified as a jump in the value of \(y_{2E}\). In this sense, the task can be reduced to an edge detection problem; this is beyond the scope of the present study. The computations corresponding to the \(G^{(2)}\) and \(G^{(5)}\) curves are shown in the top-right and bottom-left panels of Fig. 4, respectively. In the case of \(G^{(2)}\), a second focal point already appears in the initial condition plane which was not observed in the BVP computations of Hős and Champneys [52]. The two focal points are also connected with an additional \(G^{(2)}\) curve. Interestingly, the \(G^{(1)}\) curve also appears as a discontinuity in the \(y_{2E}\) values; however, in either sides \(y_{2E} \ne 0\). Therefore, this curve can be seen only as a “pale” dark red-light red transition. The reason for the non-zero velocity is that the integration is stopped at the third impact for \(G^{(2)}\) instead of at the second one required for the detection of \(G^{(1)}\). Nevertheless, an edge detection algorithm can find both the \(G^{(1)}\) and \(G^{(2)}\) curves from a single computation with \(k=2\). The grazing lines corresponding to \(k=5\) are presented in the bottom-left panel in Fig. 4. Similarly, as in the case of \(k=2\), all the previous grazing lines (\(k=1 \dots 4\)) are visible in the figure making it extremely complex. Thus, to detect the edges properly, a suitably fine resolution is necessary. This is not a problem in our case, as a single computation with one million initial conditions takes only a couple of seconds. Observe that in the bottom-left panel, no further focal points are discovered apart from the second one.

In summary, high-resolution scans of the initial conditions using our GPU accelarated IVP solver have revealed an additional feature (second focal point) of the grazing lines in the \(y_{20}-y_{30}\) initial condition plane. This shows that fast “brute force” scanning is nowadays able to discover features otherwise not visible or overlooked—here by the available BVP approach. At first sight, high-performance computation seems to be exaggerated. Even without using GPUs, the above “brute force” computations can be done within a few hours using MATLAB on a CPU. However, the main message here is that considering the usage of a “brute force” approach can lead to unexpected discoveries. Although in this specific example, high-performance computing is not really mandatory, in general, to obtain results within reasonable time for a detailed “brute force” computation, the applications of high-performance GPU (and/or CPU) clusters is usually a must.

8 Summary

In this paper, the efficacy of “brute force” technique combined with high-performance GPU computing is demonstrated through two test cases. The first model, the Keller–Miksis equation, is related to the scientific topic of sonochemistry and bubble dynamics. Apart from mapping the dynamics of bubbles to obtain approximate information about their chemical activity, the bifurcation structure of the high-resolution plots led to a discovery of a new technique to control multi-stability. The second model describes the behaviour of a pressure relief valve that can exhibit non-smooth impact dynamics. Results in the literature revealed that the grazing lines—computed via a boundary value problem solver—in the initial condition plane are organized around a spiral hub. The high-resolution scans of the initial conditions using our GPU accelerated initial value problem solver led to the discovery of a new focal point of the grazing lines. In summary, “brute force” technique can play an important role in many fields of sciences, including non-linear dynamics.

In general, the prerequisite to employ high-resolution parameter scans is a fast solver. If high computational capacities are required, a natural choice is the usage of CPU clusters that are available in many research institutes. The advantage of this approach is that highly optimised libraries are available for CPUs. However, GPUs have outstanding computational capacity/price ratio, which makes them a good alternative over CPUs. Although the parallelisation strategy for parameter scans seems to be straightforward (assign a GPU thread to each parameter combination) and libraries supporting solution of ODEs on GPUs are already available, still there can be many special issues resulting in a large performance drop.

For instance, the extremely slow CPU-GPU memory transactions need to be avoided by all costs. This can be a cumbersome task for example for systems with impact dynamics, where thousands of parallel threads (each having its own instance of the ODE with a specific parameter combination) can encounter an impact at any time. What should the programmer do if a single thread is impacting? He/she can stop the whole computation, apply the impact law on that specific thread and continue the integration process. This can be quite inefficient if the programmer has to involve CPU computations (depending on the interface and data structure of the package used), and there is always an overhead to restart the simulation as well. Thus, an efficient solver has to be able to detect impact (via event handling) and manipulate the trajectory immediately “on the fly” on the GPU for each thread selectively.

There can be several other issues that may have a negative effect on code performance if GPUs are involved. Thus, tuning up the number of the parameters is easy, but a fast and efficient GPU solver usually needs a clever implementation. Such a detailed discussion is beyond the scope of the present study. Nevertheless, our GPU code is designed to efficiently address the majority of these issues. For more details, the reader is again referred to the manual of the program package [54] and to its website www.gpuode.com or to its GitHub repository [53].