viscid-hub / Viscid

Visualizing Plasma Simulations in Python
http://viscid-hub.github.io/Viscid-docs/docs/master/index.html
Other
6 stars 6 forks source link

slicing amr_field along single-cell dimension #15

Closed liangwang0734 closed 6 years ago

liangwang0734 commented 7 years ago

For example, for a field fld of size nx ny nz = 100 1 50 and y_cc=0.

Slicing fld["y=0.f"] would return nothing from the logic in amr_field.py:_prepare_amr_slice:

            if (not np.any(np.logical_or(fld.crds.xl_nc - atol > extent[1],
                                         fld.crds.xh_nc <= extent[0]))):
                if np.any(np.isclose(fld.crds.xh_nc, extent[0], atol=atol)):
                    maybe.append(i)
                else:
                    inds.append(i)
            np.seterr(invalid=invalid_err_level)

Maybe we should replace the <= in fld.crds.xh_nc <= extent[0] to <?

liangwang0734 commented 7 years ago

I just thought a little bit more, may this was not a problem with amr_field itself but that xh is not properly passed along in vpic reader (I got this problem using the vpic reader). For example, from

ex.patches[0].crds.xl_nc, ex.patches[0].crds.xh_nc

I got

(array([  0. ,  -2.5, -32. ]), array([ 32. ,  -2.5,  32. ]))

Here, the yh and yl are both -2.5. But the global.vpc file provided delta_y=5, thus yh should be 2.5.

For example, in the `vpic.py:VPIC_File:_parse':


        block_crds = []
        for k in range(_gfile.topology[2]):
            for j in range(_gfile.topology[1]):
                for i in range(_gfile.topology[0]):
                    idx = [i, j, k]
                    xl = [_gfile.xl[d] + _gfile.ldims[d] * idx[d] * _gfile.dx[d]
                          for d in range(3)]
                    xh = [_gfile.xl[d] + _gfile.ldims[d] * (idx[d] + 1) * _gfile.dx[d]
                          for d in range(3)]
                    arrs = [np.linspace(_xl, _xh, _n)
                            for _xl, _xh, _n in zip(xl, xh, _gfile.ldims)]
                    block_crds += [coordinate.arrays2crds(arrs)]

Was line arrs = [np.linspace(_xl, _xh, _n) for _xl, _xh, _n in zip(xl, xh, _gfile.ldims)] intended to get the cc or nc coordinates? If cc, then the linspace should be passed a endpoint=false option; if nc, then the number of grid should be _gfile.ldims+1 (assuming _gfile.ldims is cell number not node number).

liangwang0734 commented 7 years ago

Kris suggested a temporality fix to mark all fields as 'cell' center and fix the crds. It would be ideal, though to specify centering depending on the variable name.