Hi!

I think that there might be something strange going on with your field definition, but it's also entirely possible that I'm misinterpreting.

The way you have it set up, I *think* that the results of the derived field (gas, elec_den) will always have the same shape as your all_data_level_max.  I believe what you may want to do is change your field definition to use data['gas','density'] instead.

Additionally, can you try:

cutout_obj['index','ones']

and comparing the shape of that to the others?  If that shape is correct then we can narrow it down to the way elec_den is defined.

Thanks,

Matt

On Sun, Jan 31, 2021 at 8:41 AM Salome Mtchedlidze <salomchedlidze@gmail.com> wrote:
Hello!

I want to create 3D field within a spherical object. For example, I create 3D electron density in the following manner:

        ds = yt.load("DATA")
        dd = ds.all_data()

        max_level  = ds.index.max_level
        ref        = 2**max_level
        low        = ds.domain_left_edge
        dims       = ds.domain_dimensions*ref
        nx, ny, nz = dims
        L  = (ds.domain_right_edge - ds.domain_left_edge).d
        all_data_level_max = ds.covering_grid(level=max_level,left_edge=low,dims=dims)

        def _elec_den(field,data):
            N_A = 1/1.6726219e-24       
            X  = 0.70                   
            return ((all_data_level_max['gas','density'])*N_A*(1+X)/2/ds.quan(1,'g') )     
        ds.add_field(('gas','elec_den'), function=_elec_den, units='auto', dimensions=dimensions.number_density)

Then if I do:
        dd['elec_den']
it will give me a 3D field which will have the same dimensions as the root grid of the dataset. However, if I create spheres and then a spherical object between these two spheres:

        rvir = 4.2454307068129566e+24                                        #virial radius of a halo [cm].
        xyz = [0.4208984375, 0.2216796875, 0.5185546875]     #centre of a halo.
        sph1 = ds.sphere(xyz, (2*rvir, 'cm'))
        sph2 = ds.sphere(xyz, (4*rvir, 'cm'))
        cutout_obj = sph2 - sph1

cutout_obj['elec_den'] or sph1['elec_den'] and sph2['elec_den']  will have the same dimensions as dd['elec_den'] ; but this cannot be true since sph1 and sph2 include less grid points.
Is there a way I can create 3D field and use it within a different yt objects?

Thanks,
Salome
_______________________________________________
yt-users mailing list -- yt-users@python.org
To unsubscribe send an email to yt-users-leave@python.org
https://mail.python.org/mailman3/lists/yt-users.python.org/
Member address: matthewturk@gmail.com