Compare results of `to_coprecessing_frame` before and after CoM correction with old and new versions of scri/sxs/...

make_plots.out

=======================================
Running in CoprecessingFrameTransients1
=======================================
/Users/boyle/.continuum/mambaforge/envs/CoprecessingFrameTransients1/lib/python3.9/site-packages/quaternion/calculus.py:310: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def fd_indefinite_integral(f, t):

python.__version__ = 3.9.16 | packaged by conda-forge | (main, Feb  1 2023, 21:42:20) 
[Clang 14.0.6 ]
numpy.__version__ = 1.21.0
scipy.__version__ = 1.7.0
h5py.__version__ = 3.8.0
numba.__version__ = 0.57.0
quaternion.__version__ = 2022.4.2
spherical_functions.__version__ = 2022.4.1
spinsfast.__version__ = 2022.4.1
scri.__version__ = 2022.8.3
quaternionic.__version__ = 1.0.6
spherical.__version__ = 1.0.13
sxs.__version__ = 2022.3.5

scri.__version__ = 2022.8.3
spherical_functions.__version__ = 2022.4.1
quaternion.__version__ = 2022.4.2
scipy.__version__ = 1.7.0
numba.__version__ = 0.57.0
numpy.__version__ = 1.21.0


=======================================
Running in CoprecessingFrameTransients2
=======================================
python.__version__ = 3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:07:22) [Clang 14.0.6 ]
numpy.__version__ = 1.24.3
scipy.__version__ = 1.10.1
h5py.__version__ = 3.8.0
numba.__version__ = 0.57.0
quaternion.__version__ = 2022.4.3
spherical_functions.__version__ = 2022.4.2
spinsfast.__version__ = 2022.4.2
scri.__version__ = 2022.8.4
quaternionic.__version__ = 1.0.6
spherical.__version__ = 1.0.13
sxs.__version__ = 2022.4.5

scri.__version__ = 2022.8.4
spherical_functions.__version__ = 2022.4.2
quaternion.__version__ = 2022.4.3
scipy.__version__ = 1.10.1
numba.__version__ = 0.57.0
numpy.__version__ = 1.24.3

make_plots.py

import sys
import numpy as np
import h5py
import sxs
import spherical_functions as sf
import scri
import quaternion
import matplotlib.pyplot as P

ident = sys.argv[1]

print("\n".join([f"{k}.__version__ = {v}" for k,v in sxs.version_info().items()]))
print()
print(scri.version_info())

def coprecessingQuatAndWaveform(t, h, lmax):
    """ Transfroms from inertial frame to coprecessing frame using scri.
    """

    w = scri.WaveformModes(
          dataType = scri.h,
          t = np.copy(t),
          data = np.copy(h.T),      # Because data gets modified internally
          ell_min = 2,
          ell_max = lmax,
          frameType = scri.Inertial,
          r_is_scaled_out = True,
          m_is_scaled_out = True,
    )

    w.to_coprecessing_frame()
    qc = w.frame
    return quaternion.as_float_array(qc).T, w.data.T


#-------------------------------------------------------------------------
def add_memory_correction(t, h, lmax, t_relax):

    # Convert the np.array to an sxs.WaveformModes object
    h_sxs = sxs.WaveformModes(
            input_array=np.copy(h.T),
            time=np.copy(t),
            time_axis=0,
            modes_axis=-1,
            ell_min=2,
            ell_max=lmax,
            spin_weight=-2)

    # Add memory corrections, but start the integral at t=t_relax
    # so it is not impacted by junk.
    h_sxs_mem = sxs.waveforms.memory.add_memory(h_sxs,
            integration_start_time=t_relax)

    #NOTE: This is currently required becase the ell=0,1 modes
    # get included by silly sxs when adding memory corrections.
    # So, we drop all modes before the first nonzero mode (2,-2).
    # This should eventually not be required if fixed in sxs, but
    # that should not break this code anyway.
    h_sxs_mem = h_sxs_mem[:, sf.LM_index(2, -2, h_sxs_mem.ell_min):]

    return h_sxs_mem.data.T


# Load raw extrapolated+com data
# SimAnnex/Public/HybTest/028/Lev3
f = h5py.File('rhOverM_Asymptotic_GeometricUnits_CoM.h5', 'r')
#f = h5py.File('HybTest/030/rhOverM_Asymptotic_GeometricUnits_CoM.h5', 'r')
h = []
lmax = 2
for ell in range(2, lmax+1):
    for emm in range(-ell, ell+1):
        t, hr, hi = f[f'Extrapolated_N2.dir/Y_l{ell}_m{emm}.dat'][()].T
        h.append(hr + 1j * hi)
h = np.array(h)

# Get copr before mem correction
qc, hcopr = coprecessingQuatAndWaveform(t, h, lmax)

# Do the mem correction
t_relax = 500
h_mem = add_memory_correction(t, h, lmax, t_relax)

# Get copr again
qc_mem, hcopr_mem = coprecessingQuatAndWaveform(t, h_mem, lmax)

q_idx = 0
P.plot(t, qc[q_idx], label='Copr quat without memory')
P.plot(t, qc_mem[q_idx], ls='dashed', label='Copr quat with memory')
P.ylabel(f'quat[{q_idx}]')
P.xlabel(f't')
P.legend()
#P.show()
P.title(f"CoprecessingFrameTransients{ident}")
P.savefig(f"plot_{ident}.png")

make_plots.sh

rsync -L meistri:/sxs-annex5/SimAnnexPreJune2023.git/Public/HybTest/030/Lev3/rhOverM_Asymptotic_GeometricUnits_CoM.h5 .

mamba create -y -n CoprecessingFrameTransients1 \
      numpy==1.21.0 scipy==1.7.0 sxs==2022.3.5 scri==2022.8.3 \
      spherical_functions==2022.4.1 quaternion==2022.4.2 matplotlib
mamba create -y -n CoprecessingFrameTransients2 \
      numpy==1.24.3 scipy==1.10.1 sxs==2022.4.5 scri==2022.8.4 \
      spherical_functions==2022.4.2 quaternion==2022.4.3 matplotlib

echo "======================================="
echo "Running in CoprecessingFrameTransients1"
echo "======================================="
mamba run -n CoprecessingFrameTransients1 python make_plots.py 1

echo
echo "======================================="
echo "Running in CoprecessingFrameTransients2"
echo "======================================="
mamba run -n CoprecessingFrameTransients2 python make_plots.py 2