<jupyter_start><jupyter_text>First just insert a single particle at (0,0,0) and test the force against the analytic expectation<jupyter_code>fakesnap = SnapShot(positions=torch.zeros((1,3),dtype=torch.float32),
                    velocities=torch.zeros((1,3),dtype=torch.float32),
                    masses=torch.full((1,),0.3,dtype=torch.float32))
potential=None
n=512
nz=512
potential=ForceGrid(n=(n,n,nz),
                    gridedges=torch.tensor([10.,10.,10.],dtype=torch.float32),
                    smoothing=0.2*20/n)
_=potential.griddata(fakesnap.positions,weights=fakesnap.masses,method='cic')
potential.grid_accelerations()<jupyter_output><empty_output><jupyter_text>Errors look small, apart from near the center, which seems to just the expected errors due to smoothing and gridding<jupyter_code>positions=torch.zeros((1000,3),dtype=torch.float32)
positions[:,0]=torch.linspace(-10,10,1000)
positions[:,1]=0 
positions[:,2]=0
acc=potential.get_accelerations(positions)

f,ax = plt.subplots(2,2,sharex='col')
ax[0,0].plot(positions[:,0].numpy(),-acc[:,0].numpy(),label='astro-dynamo')
ax[0,0].plot(potential.x.numpy(),-(potential.x*fakesnap.masses/potential.x.abs()**3).numpy(),label='Analytic',
           linestyle='--')
ax[0,0].set_ylim(-2,2)
ax[0,0].set_ylabel('Force')
ax[0,0].set_xlabel('Position')
ax[0,0].legend()
ax[0,1].set_xlim(-5,5)

ax[1,0].plot(positions[:,0].numpy(),(acc[:,0].numpy()-(positions[:,0]*fakesnap.masses/positions[:,0].abs()**3).numpy())/acc[:,0].numpy(),label='Analytic')
ax[1,0].set_ylim(-0.05,0.05)
ax[1,0].set_ylabel('Fractional difference in force')


ax[0,1].plot(positions[:,0].numpy(),-acc[:,0].numpy(),label='astro-dynamo')
ax[0,1].plot(potential.x.numpy(),-(potential.x*fakesnap.masses/potential.x.abs()**3).numpy(),label='Analytic',
           linestyle='--')
#ax[0].set_ylim(-2,2)
ax[1,1].set_xlabel('Position')
ax[0,1].set_xlim(-0.2,0.2)

ax[1,1].plot(positions[:,0].numpy(),(acc[:,0].numpy()-(positions[:,0]*fakesnap.masses/positions[:,0].abs()**3).numpy())/acc[:,0].numpy(),label='Analytic')
ax[1,1].set_ylim(-1,1)<jupyter_output><empty_output><jupyter_text>Now try with a reasonable snapshot<jupyter_code>snap=SnapShot('../inputmodels/M85_0.gz',omega=1.)
particletype = torch.full((snap.n,),ParticleType.Star,dtype=torch.uint8)
particletype[snap.particletype==0]=ParticleType.DarkMatter
snap.particletype = particletype
omega,omegaerr = astro_dynamo.analysesnap.patternspeed(snap.stars)
snap.omega = torch.Tensor([omega]).type(torch.float32)
print(f'Snapshot has pattern speed {snap.omega}')
potential=None
n=512
nz=512
potential=ForceGrid(n=(n,n,nz),
                    gridedges=torch.tensor([10.,10.,10.],dtype=torch.float32),
                    smoothing=0.2*20/n)
_=potential.griddata(snap.positions,weights=snap.masses,method='cic')
potential.grid_accelerations()
grav = mwtools.nemo.gravity_cartesian_grid(snap.as_numpy_array(), potential.x,potential.y,[0.])
positions=torch.zeros((len(potential.x),3),dtype=torch.float32)
positions[:,0]=grav['x']
positions[:,1]=grav['y'][255]
positions[:,2]=0
acc=potential.get_accelerations(positions)<jupyter_output><empty_output><jupyter_text>Things look reasonable - force differences seem small<jupyter_code>f,ax = plt.subplots(2,1,sharex='col')
ax[0].plot(grav['x'].numpy(),grav['F'][:,255,0],label='gyrfalcON')
ax[0].plot(grav['x'].numpy(),-acc[:,0].numpy(),label='astro-dynamo')
ax[0].legend()
ax[0].set_xlim(-2,2)
ax[0].set_ylabel('Force')

ax[1].plot(grav['x'].numpy(),grav['F'][:,255,0]+acc[:,0].numpy(),'.',label='Difference')
ax[1].set_ylabel('Fractional\ndifference\nin force')
ax[1].set_xlabel('Position')<jupyter_output><empty_output><jupyter_text>Finally compare an entire potential slice - differences look pretty small. There's the expected issues at the edge because the grid is finite and truncates things. And differences at the center that are probably due to the different smoothing and finite grid.<jupyter_code>from mpl_toolkits.axes_grid1 import make_axes_locatable
f,axs = plt.subplots(1,3,sharex=True,sharey=True,figsize=(16,8))

ax=axs[0]
im=ax.contourf(grav['x'],grav['y'],grav['pot'])
ax.set_aspect('equal', 'box')
divider = make_axes_locatable(ax)
cax = divider.append_axes('top', size='5%', pad=0.05)
cb = f.colorbar(im, cax=cax, orientation='horizontal')
cb.ax.xaxis.set_ticks_position('top')
cb.ax.xaxis.set_label_position('top')
cax.set_title('GyrFalcon')

ax=axs[1]
potslice = 0.5*(potential.pot[:,:,255]+potential.pot[:,:,255]).numpy()
im=ax.contourf(potential.x,potential.y,potslice)
ax.set_aspect('equal', 'box')
divider = make_axes_locatable(ax)
cax = divider.append_axes('top', size='5%', pad=0.05)
cb = f.colorbar(im, cax=cax, orientation='horizontal')
cb.ax.xaxis.set_ticks_position('top')
cb.ax.xaxis.set_label_position('top')
cax.set_title('astro-dynamo')

ax=axs[2]
diff = potslice-grav['pot'] - np.mean(potslice-grav['pot'])
im=ax.contourf(potential.x,potential.y,diff,
               vmin=-np.max(np.abs(diff)),vmax=np.max(np.abs(diff)),
              cmap=plt.cm.get_cmap('coolwarm'))
ax.set_aspect('equal', 'box')
divider = make_axes_locatable(ax)
cax = divider.append_axes('top', size='5%', pad=0.05)
cb = f.colorbar(im, cax=cax, orientation='horizontal')
cb.ax.xaxis.set_ticks_position('top')
cb.ax.xaxis.set_label_position('top')
cax.set_title('Difference')

f.tight_layout()
np.max(np.abs(diff))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20560921248142644, 'fraction_numerical': 0.03956166419019316, 'mean_word_length': 7.134441087613293, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '193536', 'license': 'MIT', 'n_tokens_mistral': 2243, 'n_tokens_neox': 2134, 'n_words': 230, 'path': 'notebooks/TestGravity.ipynb', 'repo_name': 'cwegg/astro-dynamo', 'size': '5384'}