0

I am having data visualization issue. Short summup, I'm working on a project involving a polar spherical coordinate mesh, and trying to solved coupled system of ODE (chemical reactions) for each cell. For a specific reason I need my state vector to be of the form (r^2*sin(theta)*n_i), i={1,2,3...}.

I rewrote a quick example of my issue below, you can run it as is. Why is it that cst2, which I would have assumed equal to a np.ones(a[0].shape), doesn't show a uniform pcolormesh. And even more surprising to me, why does adding the colorbar make this issue vanish ?

My best guess is that dividing by r^2*sin(theta) causes numerical issues, but how can I workaround this ? (I need to see my data without the curvature term in order to interpret it -> the jacobian division seems mandatory to me at first thought).

import numpy as np
import matplotlib.pylab as plt
fig, ax = plt.subplots()

### Edges
r = np.logspace(np.log10(1), np.log10(4.6), num=14) #cell edges
theta = np.linspace(0+0.001,np.pi-0.001,num=10) 
b = np.meshgrid(r,theta)

### Center
r_c = r[0:-1] + np.ediff1d(r)/2 #get the cell center
theta_c = theta[0:-1] + np.ediff1d(theta)/2
a = np.meshgrid(r_c,theta_c)

### The jacobian division
cst = pow(a[0],2)*np.sin(a[1])
cst2 = np.copy(cst)/pow(a[0],2)/np.sin(a[1])
pcm = ax.pcolormesh(b[0]*np.cos(b[1]),\
                    b[0]*np.sin(b[1]), \
                    cst2,cmap='seismic',edgecolor='black')
# clb = fig.colorbar(pcm, ax=ax, orientation='horizontal')
Slyphlamen
  • 26
  • 6

1 Answers1

1

Tweaking the formula slightly seems to work. Also used a lighter color for cmap.

cst2 = np.copy(cst) / (pow(a[0],2)*np.sin(a[1]))
pcm = ax.pcolormesh(b[0]*np.cos(b[1]), \
                    b[0]*np.sin(b[1]), \
                    cst2,cmap='coolwarm',edgecolor='black')

uniform color mesh

viggnah
  • 1,709
  • 1
  • 3
  • 12
  • That is not what I am aiming at, I want to divide cst by r**2sin(theta) here, and what you wrote is equivalent to cst2 = cst*np.sin(a[1])/pow(a[0], 2) From my searches I guess the issue comes from the regular machine precision issue (https://stackoverflow.com/questions/588004/is-floating-point-math-broken) which matplotlib does not treshold. But I would still be interested in knowing if there are any workarounds not involving tresholding. (Especially since colorbar fixes max and min values I've been told) – Slyphlamen Jul 09 '22 at 09:09
  • Sorry, didn't realise. I think you are right with the precision. However, replacing the two divisions with one seems to help. – viggnah Jul 12 '22 at 14:26
  • That's actually true, interesting. I guess it's definitely related to numerical precision then. It would be nice if someone could elaborate on the intrinsic mechanism causing c/a/b =! a/(a*b) before ticking the question as answered. – Slyphlamen Jul 13 '22 at 09:23
  • 1
    It comes down to limitations of trying to represent floating numbers in base 2 ([representation error](https://docs.python.org/3/tutorial/floatingpoint.html#:~:text=Representation%20error%20refers%20to%20the,exact%20decimal%20number%20you%20expect.)). Even though doing 1 division seems more _accurate_ than doing 2 divisions that may only be due to the representation errors matching up because different intermediate numbers are generated in both cases! – viggnah Jul 13 '22 at 11:13