Fortran-Contiguous points yield different results (with fix)

[MaceKuailv]

Description

Hi, I run into this benign little bug while doing Charles' homework. Basically, when I define the points in two different ways I get different results. Here is what I did:

from giverny.turbulence_dataset import *
from giverny.turbulence_toolkit import *
import numpy as np

# auth_token = 'edu.jhu.pha.turbulence.testing-201406'
auth_token = 'edu.jhu.meneveauEN530.625-37da5699'
dataset_title = 'isotropic1024coarse'
output_path = ''

dataset = turb_dataset(dataset_title = dataset_title, output_path = output_path, auth_token = auth_token)

variable = 'velocity'
temporal_method = 'none'
spatial_method = 'none'
spatial_operator = 'field'

# Now create the points in my way
n = 1024
time = 1.0

y_ = 0.9
z_ = 0.0
x = np.linspace(0,2*np.pi,n)
y = np.ones_like(x)*y_
z = np.ones_like(x)*z_
pointsa = np.vstack([x,y,z]).T.astype('float64')

# This is creating the points following the example
nx = 1024
nz = 1
n_points = nx * nz

x_points = np.linspace(0.0, 2 * np.pi, nx, dtype = np.float64)
y_points = 0.9
z_points = np.linspace(0.0, 2 * np.pi, nz, dtype = np.float64)
        
points = np.zeros((n_points, 3), dtype = np.float64)
for i in range(nx):
    for j in range(nz):
        points[i * nz + j, 0] = x_points[i]  
        points[i * nz + j, 1] = y_points
        points[i * nz + j, 2] = z_points[j]
        
# For python users, those two arrays are identical
print(np.allclose(points,pointsa))
print((points == pointsa).all())
print(points.dtype == pointsa.dtype)
print(points.shape == pointsa.shape)

# But the Result are not the same
resulta = getData(dataset, variable, time, temporal_method, spatial_method, spatial_operator, pointsa)
result  = getData(dataset, variable, time, temporal_method, spatial_method, spatial_operator, points)

For both cases, you get 4 True. However, when you plot the results, you get:

plt.plot(result[0]['ux'],label = 'Points')
plt.plot(resulta[0]['ux'], label = 'Points_my')
plt.legend()

Fix

If I just add a single line:

points_my = np.ascontiguousarray(pointsa)

this problem goes away. When you are using python, 99% of the time you don't have to worry about whether the array is F-contiguous or C-contiguous. Therefore, I think we should include this line inside the getData function.

This could also simplify the tutorials a bit, instead of writing for loops, we can do something like

points = np.array([i.ravel() for i in np.meshgrid(x_points, y_points, z_points)])