diff --git a/pySDC/helpers/spectral_helper.py b/pySDC/helpers/spectral_helper.py
index d74aab8733f256272f0cb99a5ea6e1fd13fe3a34..d2da9fbcbeb45ebeae4af27cb0ca30ea95f5cf90 100644
--- a/pySDC/helpers/spectral_helper.py
+++ b/pySDC/helpers/spectral_helper.py
@@ -1044,7 +1044,7 @@ class SpectralHelper:
         """
         Remove a BC from the matrix. This is useful e.g. when you add a non-scalar BC and then need to selectively
         remove single BCs again, as in incompressible Navier-Stokes, for instance.
-        Forward arguments for the boundary conditions using `kwargs`. Refer to documentation of 1D bases for details.
+        Forwards arguments for the boundary conditions using `kwargs`. Refer to documentation of 1D bases for details.
 
         Args:
             component (str): Name of the component the BC should act on
@@ -1311,15 +1311,15 @@ class SpectralHelper:
         """
         Get grid in spectral space
         """
-        grids = [self.axes[i].get_wavenumbers()[self.local_slice[i]] for i in range(len(self.axes))][::-1]
-        return self.xp.meshgrid(*grids)
+        grids = [self.axes[i].get_wavenumbers()[self.local_slice[i]] for i in range(len(self.axes))]
+        return self.xp.meshgrid(*grids, indexing='ij')
 
     def get_grid(self):
         """
         Get grid in physical space
         """
-        grids = [self.axes[i].get_1dgrid()[self.local_slice[i]] for i in range(len(self.axes))][::-1]
-        return self.xp.meshgrid(*grids)
+        grids = [self.axes[i].get_1dgrid()[self.local_slice[i]] for i in range(len(self.axes))]
+        return self.xp.meshgrid(*grids, indexing='ij')
 
     def get_fft(self, axes=None, direction='object', padding=None, shape=None):
         """
diff --git a/pySDC/implementations/problem_classes/Burgers.py b/pySDC/implementations/problem_classes/Burgers.py
index ef957f54d0dc33851b0688a73c1399c95ee6d96c..7ff48e16ad550f3edb7e2c9a1a451cd782246950 100644
--- a/pySDC/implementations/problem_classes/Burgers.py
+++ b/pySDC/implementations/problem_classes/Burgers.py
@@ -188,7 +188,7 @@ class Burgers2D(GenericSpectralLinear):
         components = ['u', 'v', 'ux', 'uz', 'vx', 'vz']
         super().__init__(bases=bases, components=components, spectral_space=False, **kwargs)
 
-        self.Z, self.X = self.get_grid()
+        self.X, self.Z = self.get_grid()
 
         # prepare matrices
         Dx = self.get_differentiation_matrix(axes=(0,))
diff --git a/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py b/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py
index 7a3423ec9792dd242d01494711b911879eecf1b2..718465180479cd9687a58743362d4161b4b4f349 100644
--- a/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py
+++ b/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py
@@ -264,7 +264,7 @@ class Heat2DChebychev(GenericSpectralLinear):
 
         super().__init__(bases, components, Dirichlet_recombination=False, spectral_space=False, **kwargs)
 
-        self.Y, self.X = self.get_grid()
+        self.X, self.Y = self.get_grid()
 
         I = self.get_Id()
         self.Dx = self.get_differentiation_matrix(axes=(0,))
@@ -372,7 +372,7 @@ class Heat2DUltraspherical(GenericSpectralLinear):
 
         super().__init__(bases, components, Dirichlet_recombination=False, spectral_space=False, **kwargs)
 
-        self.Y, self.X = self.get_grid()
+        self.X, self.Y = self.get_grid()
 
         self.Dxx = self.get_differentiation_matrix(axes=(0,), p=2)
         self.Dyy = self.get_differentiation_matrix(axes=(1,), p=2)
diff --git a/pySDC/implementations/problem_classes/RayleighBenard.py b/pySDC/implementations/problem_classes/RayleighBenard.py
index 69a93a3d0fd046d184dc207563d332f00e8dd122..e69a86d4962097738aa67d35cf1417203a248681 100644
--- a/pySDC/implementations/problem_classes/RayleighBenard.py
+++ b/pySDC/implementations/problem_classes/RayleighBenard.py
@@ -109,8 +109,8 @@ class RayleighBenard(GenericSpectralLinear):
         components = ['u', 'v', 'T', 'p']
         super().__init__(bases, components, comm=comm, **kwargs)
 
-        self.Z, self.X = self.get_grid()
-        self.Kz, self.Kx = self.get_wavenumbers()
+        self.X, self.Z = self.get_grid()
+        self.Kx, self.Kz = self.get_wavenumbers()
 
         # construct 2D matrices
         Dzz = self.get_differentiation_matrix(axes=(1,), p=2)
diff --git a/pySDC/tests/test_helpers/test_spectral_helper.py b/pySDC/tests/test_helpers/test_spectral_helper.py
index 627632b14b98132b28e1318eaf7465b3dafdada8..3a2f9d2136a3ae86294beffe0f7e7e453a68d222 100644
--- a/pySDC/tests/test_helpers/test_spectral_helper.py
+++ b/pySDC/tests/test_helpers/test_spectral_helper.py
@@ -22,7 +22,7 @@ def test_integration_matrix2D(nx, nz, variant, axes, useMPI=False, **kwargs):
     helper.add_axis(base='cheby', N=nz)
     helper.setup_fft()
 
-    Z, X = helper.get_grid()
+    X, Z = helper.get_grid()
 
     conv = helper.get_basis_change_matrix()
     S = helper.get_integration_matrix(axes=axes)
@@ -68,7 +68,7 @@ def test_differentiation_matrix2D(nx, nz, variant, axes, bx, bz, useMPI=False, *
     helper.add_axis(base=bz, N=nz)
     helper.setup_fft()
 
-    Z, X = helper.get_grid()
+    X, Z = helper.get_grid()
     conv = helper.get_basis_change_matrix()
     D = helper.get_differentiation_matrix(axes)
 
@@ -123,7 +123,7 @@ def test_identity_matrix2D(nx, nz, variant, bx, useMPI=False, **kwargs):
     helper.add_axis(base='cheby', N=nz)
     helper.setup_fft()
 
-    Z, X = helper.get_grid()
+    X, Z = helper.get_grid()
     conv = helper.get_basis_change_matrix()
     I = helper.get_Id()
 
@@ -215,9 +215,9 @@ def _test_transform_dealias(
     u2_hat_expect = helper.u_init_forward
     u_expect = helper.u_init
     u_expect_pad = helper_pad.u_init
-    Kz, Kx = helper.get_wavenumbers()
-    Z, X = helper.get_grid()
-    Z_pad, X_pad = helper_pad.get_grid()
+    Kx, Kz = helper.get_wavenumbers()
+    X, Z = helper.get_grid()
+    X_pad, Z_pad = helper_pad.get_grid()
 
     if axis == -2:
         f = nx // 3
@@ -476,7 +476,7 @@ def test_tau_method2D(variant, nz, nx, bc_val, bc=-1, useMPI=False, plotting=Fal
     helper.add_component(['u'])
     helper.setup_fft()
 
-    Z, X = helper.get_grid()
+    X, Z = helper.get_grid()
     x = X[:, 0]
     z = Z[0, :]
     shape = helper.init[0][1:]