Issue with orientation when solving permeability

[liamxu1]

I am wondering about the coordinate system used in the code.
I wrote a simple test code as below to calculate the permeability of a simple speicified structure:

import pumapy as puma
import pyvista as pv
import numpy as np
import os

notebook = False  # when running locally, actually open pyvista window

n = 32
v = np.ones([n,n,n], dtype=np.uint16)
for i in range(n):
    for j in range(n):
        if abs(i - j) < 3 or n - i < 5:
            v[i,j,:] = 0

ws = puma.Workspace.from_array(v)
ws.set_voxel_length(1 / n)
keff, (u_x, u_y, u_z) = puma.compute_permeability(ws, (1, 1), tol=1e-7, maxiter=10000, solver_type='minres', matrix_free=True, precondition=False)

ws.voxel_length = 1
px = puma.render_orientation(u_x, scale_factor=1e3, plot_directly=False)
puma.render_volume(ws, cutoff=(1, ws.max()), solid_color=(0, 0, 255), add_to_plot=px)

and the output reads:

'import TexGen.Core' failed (TexGen is only made available when installing puma with conda on UNIX).
Approximate memory requirement for simulation: 1.63 MB
Initializing indexing matrices ... Done
Creating A matrix ... Done
Time to setup system: 0.39247389999218285
Running x direction
Creating b vector ... Done
Solving Ax=b using minres solver
Iteration: 192, driving either residual (0.0000000886, 0.0120472855) --> target = 0.0000001000 ... Done
Running y direction
Creating b vector ... Done
Solving Ax=b using minres solver
Iteration: 199, driving either residual (0.0000000864, 0.0190863559) --> target = 0.0000001000 ... Done
Running z direction
Creating b vector ... Done
Solving Ax=b using minres solver
Iteration: 15, driving either residual (0.0000000473, 0.0263754847) --> target = 0.0000001000 ... Done

Effective permeability tensor:
[[ 7.30006300e-05 -3.79325314e-05  6.41813365e-18]
 [-3.79320670e-05  2.09240023e-04  4.20316401e-20]
 [-1.06532747e-12 -1.87436839e-12  2.73029620e-04]]
Time to solve: 10.284912599949166

However as shown in the figure, when a unit force is posed on +x direction, the velocity on +y direction should be positive, which leads to a positive P[0, 1] component. This is contradictory to the code output. I am wondering whether I made somthing wrong. And I find "orientation" functions in the documentation, is that related to this phenomenon?