adapted for different floating point precisions (use with fpmath)

examples/test_cheb_mo_1p_1d_twostep.py

@@ -5,13 +5,20 @@ from rotagaporp.chebyshev import Chebyshev, TwoStepChebyshev
 from rotagaporp.symplectic import Symplectic
 from rotagaporp.systems import Morse
 from liouville import SPLPQP
+from generic import spfloat
 
-params = {'tstart': 0., 'tend': 500., 'tstep': 0.05, 'nterms': 3, 'tqdm': True}
-syspar0 = {'q0': 3., 'p0': 0.}
+prec = 8
+
+params = {'tstart': spfloat(0., prec), 'tend': spfloat(500., prec),
+          'tstep': spfloat(0.05, prec), 'nterms': 5, 'tqdm': True}
+syspar0 = {'q0': spfloat(3., prec), 'p0': spfloat(0., prec),
+           'fpprec': 'multi', 'prec': prec}
 system = Morse(**syspar0)
 
 # find the initial values for the two step propagator
 ini_params = params.copy()
+ini_params['fpprec'] = 'multi'
+ini_params['prec'] = prec
 ini_params['tend'] = params['tstart'] + params['tstep']
 ini_params['tstep'] = 0.001
 ini_params['nterms'] = 5
@@ -20,9 +27,10 @@ ini_prop.propagate()
 ini_prop.analyse()
 print('energy drift of initializer', np.max(np.fabs(ini_prop.er)))
 
-assert np.isclose(ini_prop.time[-1], params['tstart'] + params['tstep'])
+assert np.isclose(float(ini_prop.time[-1]), float(ini_params['tend']))
 assert len(ini_prop.sq) == len(ini_prop.time)
-syspar = {'q0': ini_prop.sq[-1], 'p0': ini_prop.sp[-1]}
+syspar = {'q0': ini_prop.sq[-1], 'p0': ini_prop.sp[-1], 'fpprec': 'multi',
+          'prec': prec}
 print(syspar)
 system = Morse(**syspar)
 
@@ -32,26 +40,26 @@ del params_verlet['nterms']
 vv_prop = Symplectic(syst=system, **params_verlet)
 vv_prop.propagate()
 vv_prop.analyse()
-edrift_vv = vv_prop.er
 
 # one-step propagator
+params['fpprec'] = 'multi'
+params['prec'] = prec
 ch_prop = Chebyshev(syst=system, **params)
 ch_prop.propagate()
 time_ch, q_ch, p_ch = ch_prop.get_trajectory()
 ch_prop.analyse()
-edrift_ch = ch_prop.er
 
 # two-step propagator
-params['q_1'] = syspar0['q0']
-params['p_1'] = syspar0['p0']
+params['q_1'] = spfloat(syspar0['q0'], prec)
+params['p_1'] = spfloat(syspar0['p0'], prec)
 params['liouville'] = SPLPQP
 params['parity'] = 0
 params['signed'] = False
+
 ts_prop = TwoStepChebyshev(syst=system, **params)
 ts_prop.propagate()
 time_ts, q_ts, p_ts = ts_prop.get_trajectory()
 ts_prop.analyse()
-edrift_ts = ts_prop.er
 
 assert len(time_ts) == len(time_ch)
 assert np.allclose(time_ts, time_ch)
@@ -59,9 +67,9 @@ q_an, p_an = system.solution(times=time_ts)
 
 fig = plt.figure()
 plot1 = fig.add_subplot(311)
-plot1.plot(time_ch, np.maximum.accumulate(edrift_ch), label='Chebyshev')
-plot1.plot(time_ts, np.maximum.accumulate(edrift_ts), label='Two-step Chebyshev')
-plot1.plot(time_ts, np.maximum.accumulate(edrift_vv), label='Velocity Verlet')
+plot1.plot(time_ch, np.maximum.accumulate(ch_prop.er), label='Chebyshev')
+plot1.plot(time_ts, np.maximum.accumulate(ts_prop.er), label='Two-step Chebyshev')
+plot1.plot(time_ts, np.maximum.accumulate(vv_prop.er), label='Velocity Verlet')
 plot1.set_ylabel('maximum energy drift')
 plot1.set_yscale('log')
 plot1.set_xscale('log')