Memory Consumption Check

[cindytsai]

Memory Consumption/Leakage Check

• Check each libyt API memory consumption and see if there is any memory leakage.

• Memory Consumption Check:

  • Using large for loops to run each API.
  • Use valgrind tool massif.
  • Check with MPI.
  • Execute bash script:

    # run with "sh run.sh <log_folder_name>"
    # naming
    LOG_FOLDER=$1
    EXE=API-test.out
    # execute
    export LD_PRELOAD=/home/calab912/software/valgrind-3.19.0/lib/valgrind/libmpiwrap-amd64-linux.so
    export PYTHONMALLOC=malloc
    mpirun -np 4 --output-filename $LOG_FOLDER /home/calab912/software/valgrind-3.19.0/bin/valgrind -v --tool=massif --time-unit=B --detailed-freq=1 ./$EXE 

  • Makefile:

    
    # path
    MPI_PATH   := /home/calab912/software/openmpi/3.1.5-gnu
    LIBYT_PATH := /home/calab912/Documents/GitHub/libyt

  output name

  BIN := API-test.out FILE := main.cpp COMPILER := $(MPI_PATH)/bin/mpic++

  command

  $(BIN): $(FILE) $(COMPILER) -o $(BIN) $(FILE) -I$(LIBYT_PATH)/include -L$(LIBYT_PATH)/lib -lyt

  clean: rm -f $(BIN) rm -rf log rm -rf pycache rm -f .png .gif RecordTime_ massif.

• Python Memory Leakage Check:
  • Count every libyt data member's reference count using sys.getrefcount.
  • sys.getrefcount will increase original count by 1, since it refers to that value to get the count.(link)
  • Every value should have reference count 2.

Memory Consumption Caused By Each API

• yt_init and yt_finalize
• yt_set_parameter
• yt_add_user_parameter_*
• yt_set_parameter and yt_get_gridsPtr
• yt_set_parameter, yt_get_gridsPtr, yt_get_fieldsPtr, and yt_commit_grids
• yt_set_parameter, yt_get_fieldsPtr and yt_commit_grids
• yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and RMA Process
• yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and load yt
• yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and load yt with GAMER hierarchy

Other yt Operations

Covering Grid

• yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and then draw covering_grids
• yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, load gamer hierarchy, and then draw covering_grids
• [yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, load gamer hierarchy, and then slice target covering_grid to smaller unit]()

Python API check

• [x] yt_init()
  • [x] PyInit_libyt (init_libyt_module.cpp)
  • libyt, grid_data, hierarchy, param_yt, param_user
  • [x] libyt_field_derived_func
  • Return 3d NumPy array owned by Python.
  • [x] libyt_particle_get_attr
  • Return 1d NumPy array owned by Python.
  • [x] libyt_field_get_field_remote
  • Return data[gid]["field"][3d-NumPyArray]
  • [x] libyt_particle_get_attr_remote
  • Return data[gid]["ptype"]["attr"][1d-NumPyArray]
• [x] yt_set_parameter()
  • [x] add_dict_string (add_dict.cpp)
  • param_yt
    • frontend, fig_basename
  • [x] add_dict_scalar (add_dict.cpp)
  • param_yt
    • current_time, current_redshift, omega_lambda, omega_matter, hubble_constant, length_unit, mass_unit, time_unit, magnetic_unit, cosmological_simulation, dimensionality, refine_by, num_grids, domain_left_edge, domain_right_edge, periodicity, domain_dimensions
  • [x] allocate_hierarchy() (allocate_hierarchy.cpp)
  • hierarchy
    • grid_left_edge, grid_right_edge, grid_dimensions, grid_particle_count, grid_parent_id, grid_levels, proc_num
• [x] yt_add_user_parameter_*
  • [x] add_dict_string
  • [x] add_dict_scalar
  • [x] add_dict_vector3
• [x] yt_commit_grids()
  • [x] add_dict_field_list (add_dict.cpp)
  • Dictionary param_yt["field_list"]
    • Key-Value: field name-Dictionary
    • "attribute": [ , [,], ]
    • "field_define_type":
    • "swap_axes": true/false
    • "ghost_cell": [ ]
  • [x] add_dict_particle_list (add_dict.cpp)
  • Dictionary param_yt["particle_list"]
    • Key-Value: ptype-Dictionary
    • "attribute": Dictionary
      • : [ , [], ]
    • "particle_coor_label": [ , , ]
  • [x] append_grid (append_grid.cpp)
  • grid_data
    • grid_data[gid]["field"]
• [x] yt_free()
• [x] yt_finalize()

Python Reference Count Check

Under these structure, sys.getrefcount are all 2.

• libyt: 9
  • Dictionary
  • grid_data: 3
  • hierarchy: 3
  • param_yt: 3
  • param_user: 3
  • Method
  • derived_func: 3
  • get_attr: 3
  • get_field_remote: 3
  • get_attr_remote: 3

[cindytsai]

yt_init and yt_finalize only

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

import yt
yt.enable_parallelism()
def yt_inline():
    pass

• Result:

  • Total Memory Peak Caused by yt_init:

  

[cindytsai]

yt_set_parameter

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* testing API. */
    long iter = 1;
    int  num_grids = 2000000;

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "FigName";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 1;   
    for (int d=0; d<3; d++)
    {
        param_yt.domain_dimensions[d] = 40;
        param_yt.domain_left_edge [d] = 0.0;
        param_yt.domain_right_edge[d] = 1.0;
        param_yt.periodicity      [d] = 0;
    }
    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

• Result:

  • num_grids: 2e6
  • Hierarhcy array with length num_grids allocated by libyt: ~ 2e6 * 100 byte ~ 186 MB
  • grid_left_edge: ( num_grids, 3 ) NPY_DOUBLE
  • grid_right_edge: ( num_grids, 3 ) NPY_DOUBLE
  • grid_dimensions: ( num_grids, 3 ) NPY_LONG
  • grid_particle_count: ( num_grids, 1 ) NPY_LONG
  • grid_parent_id: ( num_grids, 1 ) NPY_LONG
  • grid_levels: ( num_grids, 1 ) NPY_LONG
  • proc_num: ( num_grids, 1) NPY_INT
  • Run with 1 iteration:

  

[cindytsai]

yt_add_user_parameter_*

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* Testing api */
    long iter = 1;

    for(long t=0; t<iter; t++){
        const  int   user_int        = 1;
        const long   user_long       = 2;
        const uint   user_uint       = 3;
        const ulong  user_ulong      = 4;
        const float  user_float      = 5.0;
        const double user_double     = 6.0;
        const char  *user_string     = "my_string";
        const int    user_int3   [3] = { 7, 8, 9 };
        const double user_double3[3] = { 10.0, 11.0, 12.0 };

        yt_add_user_parameter_int   ( "user_int",     1, &user_int     );
        yt_add_user_parameter_long  ( "user_long",    1, &user_long    );
        yt_add_user_parameter_uint  ( "user_uint",    1, &user_uint    );
        yt_add_user_parameter_ulong ( "user_ulong",   1, &user_ulong   );
        yt_add_user_parameter_float ( "user_float",   1, &user_float   );
        yt_add_user_parameter_double( "user_double",  1, &user_double  );
        yt_add_user_parameter_string( "user_string",      user_string  );

        yt_add_user_parameter_int   ( "user_int3",    3,  user_int3    );
        yt_add_user_parameter_double( "user_double3", 3,  user_double3 );

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

• Result:

  • Memory consumption does not change much.
  • Run with 1 iteration.

  

[cindytsai]

yt_set_parameter and yt_get_gridsPtr

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* testing API. */
    long iter = 1;
    int  num_grids = 2000000;
    int  num_grids_local = num_grids / nrank;

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "FigName";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 3;   
    for (int d=0; d<3; d++)
    {
        param_yt.domain_dimensions[d] = 40;
        param_yt.domain_left_edge [d] = 0.0;
        param_yt.domain_right_edge[d] = 1.0;
        param_yt.periodicity      [d] = 0;
    }
    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local;

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );

        yt_grid *grids_local;
        yt_get_gridsPtr( &grids_local );

        for(long lid=0; lid<num_grids_local; lid++){
            for(int d=0; d<3; d++){
                grids_local[lid].left_edge[d]  = 0.0;
                grids_local[lid].right_edge[d] = 1.0;
                grids_local[lid].grid_dimensions[d] = 16;
            }
            grids_local[lid].id = lid;
            grids_local[lid].parent_id = -1;
            grids_local[lid].level = 0;
        }

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

• Result:

  • num_grids: 2e6
  • num_grids_local: 2e6 / 4
  • This libyt API is responsible for allocating yt_grid array with length num_grids_local for user to fill in: ~ (2e6 / 4) * 116 byte ~ 54 MB
  • left_edge[3]: double
  • right_edge[3]: double
  • particle_count_list: long*
  • grid_particle_count: long
  • id: long
  • parent_id: long
  • grid_dimensions[3]: int
  • level: int
  • proc_num: int
  • field_data: yt_data*
  • Run with 1 iteration.
  • Only 54 MB of memory used by this API, the other are used by yt_set_parameter.

  

[cindytsai]

yt_set_parameter, yt_get_gridsPtr, yt_get_fieldsPtr, and yt_commit_grids

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

void DerivedFunc(int len, long *list_gid, yt_array *data_array){
    for(int l=0; l<len; l++){
        for(int idx=0; idx<(16*16*16); idx++){ ((float*)data_array[l].data_ptr)[idx] = 1.0; }
    }
}

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* testing API. */
    long iter = 1;
    int  num_grids = 2000000;
    int  num_grids_local = num_grids / nrank;

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "FigName";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 1;   
    for (int d=0; d<3; d++)
    {
        param_yt.domain_dimensions[d] = 40;
        param_yt.domain_left_edge [d] = 0.0;
        param_yt.domain_right_edge[d] = 1.0;
        param_yt.periodicity      [d] = 0;
    }
    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local;

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );

        yt_field *field_list;
        yt_get_fieldsPtr( &field_list );
        field_list[0].field_name = "Dens";
        field_list[0].field_define_type = "derived_func";
        field_list[0].swap_axes = true; // default is true, which means data is represented in [z][y][x] coordinate.
        field_list[0].field_dtype = YT_DOUBLE;
        char *field_name_alias[] = {"Name Alias 1", "Name Alias 2", "Name Alias 3"};
        field_list[0].field_name_alias = field_name_alias;
        field_list[0].num_field_name_alias = 3;
        field_list[0].derived_func = DerivedFunc;

        yt_grid *grids_local;
        yt_get_gridsPtr( &grids_local );

        for(long lid=0; lid<num_grids_local; lid++){
            for(int d=0; d<3; d++){
                grids_local[lid].left_edge[d]  = 0.0;
                grids_local[lid].right_edge[d] = 1.0;
                grids_local[lid].grid_dimensions[d] = 16;
            }
            grids_local[lid].id = lid + myrank * num_grids_local;
            grids_local[lid].parent_id = -1;
            grids_local[lid].level = 0;
        }

        yt_commit_grids();

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

• Result:

  • num_grids: 2e6
  • MPI Size: 1
  • Compare to MPI Size 4, see memory consumption difference between MPI size 4 and 1.

  

• Result:

  • num_grids: 2e6
  • num_grids_local: 2e6 / 4
  • MPI Size: 4
  • yt_commit_grids is responsible for collecting and broadcasting hierarchy stored in each rank.
  • Memory for collecting hierarchy in each rank, they will be freed after the job is done:
  • An additional full hierarchy is allocated. ~ 2e6 x 92 byte ~ 170 MB
  • A local hierarchy is also allocated for MPI data transmission. ~ (2e6 / 4) x 92 byte ~ 42.5 MB
  • Data structure of yt_hierarchy: 92 byte
    • left_edge[3]: double
    • right_edge[3]: double
    • grid_particle_count: long
    • id: long
    • parent_id: long
    • dimensions[3]: int
    • level: int
    • proc_num: int
  • Memory for creating key-value pairs in dictionary grid_data. They will stay alive, until yt_free_grids is called.
  • The test here does not input any data pointer, so Python None object is added. Even if we really input a data pointer, the memory won't grow. Because data pointer will first be wrapped by NumPy API, then pointed by dictionary.
  • Looks like creating Python objects takes quite some memory.
  • Run with 1 iteration.

  

[cindytsai]

yt_set_parameter, yt_get_fieldsPtr, and yt_commit_grids

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* testing API. */
    long iter = 1;
    int  num_grids = 2000000;
    int  num_grids_local = num_grids / nrank;

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "FigName";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 1;   
    for (int d=0; d<3; d++)
    {
        param_yt.domain_dimensions[d] = 40;
        param_yt.domain_left_edge [d] = 0.0;
        param_yt.domain_right_edge[d] = 1.0;
        param_yt.periodicity      [d] = 0;
    }
    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local;

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );

        yt_field *field_list;
        yt_get_fieldsPtr( &field_list );
        field_list[0].field_name = "Dens";
        field_list[0].field_define_type = "derived_func";
        field_list[0].swap_axes = true; // default is true, which means data is represented in [z][y][x] coordinate.
        field_list[0].field_dtype = YT_DOUBLE;
        char *field_name_alias[] = {"Name Alias 1", "Name Alias 2", "Name Alias 3"};
        field_list[0].field_name_alias = field_name_alias;
        field_list[0].num_field_name_alias = 3;

        yt_commit_grids();

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

• Result:

  • num_grids: 2e8
  • The memory consumption by this API is really low. Memory is used by yt_set_parameter. As you can see, the graph is identical in yt_set_parameter.

  

• Result:

  • num_grids: 2e6
  • The memory consumption by this API is really low. Memory is used by yt_set_parameter. As you can see, the graph is identical in yt_set_parameter.

  

[cindytsai]

yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and RMA Process

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

void DerivedFunc(int len, long *list_gid, yt_array *data_array){
    for(int l=0; l<len; l++){
        for(int idx=0; idx<(16*16*16); idx++){ ((float*)data_array[l].data_ptr)[idx] = 1.0; }
    }
}

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* testing API. */
    long iter = 1;
    int  num_grids = 2000000; // cannot alter arbitrary.
    int  num_grids_local = num_grids / nrank;

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "FigName";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 1;   
    param_yt.domain_dimensions[0]    = 100*16;
    param_yt.domain_dimensions[1]    = 100*16;
    param_yt.domain_dimensions[2]    = 200*16;

    double space[3] = { 100.0 / param_yt.domain_dimensions[0] ,
                        100.0 / param_yt.domain_dimensions[1] , 
                        100.0 / param_yt.domain_dimensions[2] };

    for (int d=0; d<3; d++)
    {
        param_yt.domain_left_edge [d] = 0.0;
        param_yt.domain_right_edge[d] = 100.0;
        param_yt.periodicity      [d] = 0;
    }

    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local;

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );
        const int mhd = 0;   
        yt_add_user_parameter_int("mhd", 1, &mhd);
        const int srhd = 0;
        yt_add_user_parameter_int("srhd", 1, &srhd);

        yt_field *field_list;
        yt_get_fieldsPtr( &field_list );
        field_list[0].field_name = "Dens";
        field_list[0].field_define_type = "derived_func";
        field_list[0].swap_axes = true; // default is true, which means data is represented in [z][y][x] coordinate.
        field_list[0].field_dtype = YT_DOUBLE;
        char *field_name_alias[] = {"Name Alias 1", "Name Alias 2", "Name Alias 3"};
        field_list[0].field_name_alias = field_name_alias;
        field_list[0].num_field_name_alias = 3;
        field_list[0].derived_func = DerivedFunc;

        yt_grid *grids_local;
        yt_get_gridsPtr( &grids_local );

        for(long lid=0; lid<num_grids_local; lid++){
            long gid = lid + myrank * num_grids_local;
            int grid_idx[3];
            grid_idx[2] = gid / (param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1]); 
            grid_idx[1] = (gid - (grid_idx[2] * param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1])) / param_yt.domain_dimensions[0];
            grid_idx[0] = gid - (grid_idx[2] * param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1]) - (grid_idx[1] * param_yt.domain_dimensions[0]);
            for(int d=0; d<3; d++){
                grids_local[lid].left_edge[d]  = space[d] * (double)grid_idx[d];
                grids_local[lid].right_edge[d] = grids_local[lid].left_edge[d] + space[d];
                grids_local[lid].grid_dimensions[d] = 16;
            }
            grids_local[lid].id = gid;
            grids_local[lid].parent_id = -1;
            grids_local[lid].level = 0;
        }

        yt_commit_grids();

        yt_inline("yt_inline");

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

import yt
from mpi4py import MPI
import time

yt.enable_parallelism()
myrank = MPI.COMM_WORLD.Get_rank()
num_grids = 2000000
num_exchange = 4500

def yt_inline():
    import libyt

    # Mock prepare gid list and nonlocal id list
    if myrank == 0:
        to_prepare = list(range(num_exchange))
        nonlocal_id, nonlocal_rank = [], []
    else:
        to_prepare = []
        nonlocal_id, nonlocal_rank = list(range(num_exchange)), [0] * num_exchange

    # libyt RMA
    print("RMA process...")
    nonlocal_data = libyt.get_field_remote([b"Dens"], 1, to_prepare, len(to_prepare), nonlocal_id, nonlocal_rank, len(nonlocal_id))
    time.sleep(3)

• Result:

  • num_grids: 2e6

  • Number of grids to exchange at RMA process: 4500

  • Root rank: prepare 4500 grids.

  • Non-root rank: Fetch 4500 grids from root rank.

  • RMA Process:

  • 1. Exchange data info: ~ 4500 x 34 byte ~ 0.00014 GB
       • id: long
       • address: MPI_Aint (long)
       • rank: int
       • data_dtype: short
       • data_dim[3]: int
  • 2. Prepare local data and fetched remote data: ~ 4500 x 16^3 x 8 byte ~ 137 MB

  • Run with 1 iteration:

  • Each rank has roughly the same memory consumption graph.

  • The first peak is caused by yt_commit_grids. The second memory peak is caused by RMA process.

  • Root rank:

  • See yt_rma_field::prepare_data label. It takes roughly 131 MiB

  • Non-root rank:

  • See yt_rma_field::fetch_remote_data label. It takes roughly 128 MiB.

[cindytsai]

yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and load yt

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

void DerivedFunc(int len, long *list_gid, yt_array *data_array){
    for(int l=0; l<len; l++){
        for(int idx=0; idx<(16*16*16); idx++){ ((float*)data_array[l].data_ptr)[idx] = 1.0; }
    }
}

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "print_stats";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* testing API. */
    long iter = 1;
    int  num_grids = 2000000; // cannot alter arbitrary.
    int  num_grids_local = num_grids / nrank;

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "FigName";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 1;   
    param_yt.domain_dimensions[0]    = 100*16;
    param_yt.domain_dimensions[1]    = 100*16;
    param_yt.domain_dimensions[2]    = 200*16;

    double space[3] = { 100.0 / param_yt.domain_dimensions[0] ,
                        100.0 / param_yt.domain_dimensions[1] , 
                        100.0 / param_yt.domain_dimensions[2] };

    for (int d=0; d<3; d++)
    {
        param_yt.domain_left_edge [d] = 0.0;
        param_yt.domain_right_edge[d] = 100.0;
        param_yt.periodicity      [d] = 0;
    }

    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local;

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );
        const int mhd = 0;   
        yt_add_user_parameter_int("mhd", 1, &mhd);
        const int srhd = 0;
        yt_add_user_parameter_int("srhd", 1, &srhd);

        yt_field *field_list;
        yt_get_fieldsPtr( &field_list );
        field_list[0].field_name = "Dens";
        field_list[0].field_define_type = "derived_func";
        field_list[0].swap_axes = true; // default is true, which means data is represented in [z][y][x] coordinate.
        field_list[0].field_dtype = YT_DOUBLE;
        char *field_name_alias[] = {"Name Alias 1", "Name Alias 2", "Name Alias 3"};
        field_list[0].field_name_alias = field_name_alias;
        field_list[0].num_field_name_alias = 3;
        field_list[0].derived_func = DerivedFunc;

        yt_grid *grids_local;
        yt_get_gridsPtr( &grids_local );

        for(long lid=0; lid<num_grids_local; lid++){
            long gid = lid + myrank * num_grids_local;
            int grid_idx[3];
            grid_idx[2] = gid / (param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1]); 
            grid_idx[1] = (gid - (grid_idx[2] * param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1])) / param_yt.domain_dimensions[0];
            grid_idx[0] = gid - (grid_idx[2] * param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1]) - (grid_idx[1] * param_yt.domain_dimensions[0]);
            for(int d=0; d<3; d++){
                grids_local[lid].left_edge[d]  = space[d] * (double)grid_idx[d];
                grids_local[lid].right_edge[d] = grids_local[lid].left_edge[d] + space[d];
                grids_local[lid].grid_dimensions[d] = 16;
            }
            grids_local[lid].id = gid;
            grids_local[lid].parent_id = -1;
            grids_local[lid].level = 0;
        }

        yt_commit_grids();

        yt_inline("yt_inline");

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

Put gc.collect() At the Back

import yt, sys, gc
import numpy as np

yt.enable_parallelism()

def yt_inline(  ):
    ds = yt.frontends.libyt.libytDataset()
    ds.print_stats()
    del(ds)
    gc.collect()

def yt_inline_inputArg( field ):
    pass

• Result:

  • num_grids: 2e6

  • MPI Size: 4

  • Every grid is in root level.

  • Only root rank executes ds.print_stats.

  • valgrind massif:

  • Root rank: (peak ~ 12.4 GB)

    

  • Non-root rank: (peak ~ 1.1 GB)

    

Put gc.collect() At the Front

import yt, sys, gc
import numpy as np

yt.enable_parallelism()

def yt_inline(  ):
    gc.collect()
    ds = yt.frontends.libyt.libytDataset()
    ds.print_stats()

def yt_inline_inputArg( field ):
    pass

• Result:

  • num_grids: 2e6

  • MPI Size: 4

  • valgrind massif:

  • Root rank: (peak ~ 12.4 GB)

    • There is still a huge amount of memory on _PyObject_GC_Malloc.

  • Non-root rank: (peak ~ 1.1 GB)

    

No gc

import yt
yt.enable_parallelism()

def yt_inline(  ):
   ds = yt.frontends.libyt.libytDataset()
   ds.print_stats()

def yt_inline_inputArg( field ):
   pass

• Result:

  • Root rank: (peak ~ 12.4 GB)

  

  • Non-root rank: (peak ~ 1.1 GB)

  

[cindytsai]

yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and load yt with GAMER hierarchy

#include "libyt.h"
#include <mpi.h>
#include <iostream>
#include <fstream>

void DerivedFunc(int len, long *list_gid, yt_array *data_array){
    for(int l=0; l<len; l++){
        for(int idx=0; idx<(16*16*16); idx++){ ((float*)data_array[l].data_ptr)[idx] = 1.0; }
    }
}

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "inline";
    param_libyt.check_data = true;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* Load Hierarchy */
    std::fstream myfile("./Hierarchy.txt", std::ios_base::in);
    double left_edge[3], right_edge[3];
    int grid_dim[3], level;
    long id, parent_id;

    /* testing API. */
    long iter = 1;
    int  num_grids = 1823232; // TODO
    int  *num_grids_local = new int  [nrank];
    long *acc_grids       = new long [nrank+1];
    if( num_grids % nrank != 0 ){
        for(int r=0; r<nrank; r++){
            if( r == nrank - 1 ) num_grids_local[r] = num_grids - num_grids_local[0] * (nrank - 1);
            else                 num_grids_local[r] = (num_grids / nrank) + 1;
        }
    }
    else{
        for(int r=0; r<nrank; r++){ num_grids_local[r] = num_grids / nrank; }
    }

    for(int r=0; r<nrank+1; r++){
        acc_grids[r] = 0;
        for(int r1=0; r1<r; r1++){
            acc_grids[r] += num_grids_local[r1];
        }
    }

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "FigName";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 4;   
    param_yt.domain_dimensions[0]    = 576; // TODO
    param_yt.domain_dimensions[1]    = 576; // TODO
    param_yt.domain_dimensions[2]    = 576; // TODO

    for (int d=0; d<3; d++)
    {
        param_yt.domain_left_edge [d] = 0.0;         // TODO
        param_yt.domain_right_edge[d] = 0.15823011;  // TODO
        param_yt.periodicity      [d] = 0;
    }

    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local[myrank];

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );
        const int mhd = 0;   
        yt_add_user_parameter_int("mhd", 1, &mhd);
        const int srhd = 0;
        yt_add_user_parameter_int("srhd", 1, &srhd);

        yt_field *field_list;
        yt_get_fieldsPtr( &field_list );
        char *field_name[] = {"REAL", "IMAG", "DENS", "POTE"};
        for(int v=0; v<param_yt.num_fields; v++){
            field_list[v].field_name = field_name[v];
            field_list[v].field_define_type = "derived_func";
            field_list[v].swap_axes = true; // default is true, which means data is represented in [z][y][x] coordinate.
            field_list[v].field_dtype = YT_DOUBLE;
            field_list[v].derived_func = DerivedFunc;
        }

        yt_grid *grids_local;
        yt_get_gridsPtr( &grids_local );

        int lid = 0;
        for(long i=0; i<num_grids; i++){
            // Read file
            myfile >> id;
            myfile >> parent_id;
            myfile >> level;
            myfile >> left_edge[0]  >> left_edge[1]  >> left_edge[2];
            myfile >> right_edge[0] >> right_edge[1] >> right_edge[2];
            myfile >> grid_dim[0]   >> grid_dim[1]   >> grid_dim[2];

            // Write to grids_local in different rank respectively.
            if(acc_grids[myrank] <= i && i < acc_grids[myrank+1]){
                for(int d=0; d<3; d++){
                    grids_local[lid].left_edge[d]  = left_edge[d];
                    grids_local[lid].right_edge[d] = right_edge[d];
                    grids_local[lid].grid_dimensions[d] = grid_dim[d];
                }
                grids_local[lid].id = id;
                grids_local[lid].parent_id = parent_id;
                grids_local[lid].level = level;
                lid += 1;
            }
        }

        yt_commit_grids();

        yt_inline("yt_inline");

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();
    delete [] num_grids_local;
    delete [] acc_grids;

    return 0;
}

import yt, sys, gc
import numpy as np

yt.enable_parallelism()

def yt_inline(  ):
   ds = yt.frontends.libyt.libytDataset()
   ds.print_stats()
   del(ds)
   gc.collect()

def yt_inline_inputArg( field ):
   pass

• Result:

  • MPI Size: 4
  • This test directly load a gamer hierarchy.
  • File IO for reading gamer's hierarchy only takes about 71 KB.
  • Adding more fields meaning adding more key-value pair to grid_data does not alter memory consumption that much, though it is still a waste.
  • ds.print_stats() executes on root rank.
  • Root rank: peak~ 12.4 GB

  

  • Non-root rank: peak~ 1.1 GB

  

[cindytsai]

yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, and then draw covering_grids

#include "libyt.h"
#include <mpi.h>
#include <unistd.h>
#include <stdlib.h>

void DerivedFunc(int len, long *list_gid, yt_array *data_array){
    for(int l=0; l<len; l++){
        for(int idx=0; idx<(16*16*16); idx++){ ((float*)data_array[l].data_ptr)[idx] = 1.0; }
    }
}

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "test_covering_grid_no_slabs";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* testing API. */
    long iter = 1;
    int  num_grids = 250000; // cannot alter arbitrary.
    int  num_grids_local = num_grids / nrank;

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "Fig";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 4;   
    param_yt.domain_dimensions[0]    = 100*8;  // cannot alter arbitrary.
    param_yt.domain_dimensions[1]    = 100*8;  // cannot alter arbitrary.
    param_yt.domain_dimensions[2]    = 200*8;  // cannot alter arbitrary.

    double space[3] = { 100.0 / param_yt.domain_dimensions[0] ,
                        100.0 / param_yt.domain_dimensions[1] , 
                        100.0 / param_yt.domain_dimensions[2] };

    for (int d=0; d<3; d++)
    {
        param_yt.domain_left_edge [d] = 0.0;
        param_yt.domain_right_edge[d] = 100.0;
        param_yt.periodicity      [d] = 0;
    }

    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local;

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );
        const int mhd = 0;   
        yt_add_user_parameter_int("mhd", 1, &mhd);
        const int srhd = 0;
        yt_add_user_parameter_int("srhd", 1, &srhd);

        yt_field *field_list;
        yt_get_fieldsPtr( &field_list );

        char *FieldLabels[] = {"Dens", "MomX", "MomY", "MomZ"};
        for(int v=0; v<param_yt.num_fields; v++){
            field_list[v].field_name = FieldLabels[v];
            field_list[v].field_define_type = "derived_func";
            field_list[v].swap_axes = true;
            field_list[v].field_dtype = YT_DOUBLE;;
            field_list[v].derived_func = DerivedFunc;
        }

        yt_grid *grids_local;
        yt_get_gridsPtr( &grids_local );

        for(long lid=0; lid<num_grids_local; lid++){
            long gid = lid + myrank * num_grids_local;
            int grid_idx[3];
            grid_idx[2] = gid / (param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1]); 
            grid_idx[1] = (gid - (grid_idx[2] * param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1])) / param_yt.domain_dimensions[0];
            grid_idx[0] = gid - (grid_idx[2] * param_yt.domain_dimensions[0] * param_yt.domain_dimensions[1]) - (grid_idx[1] * param_yt.domain_dimensions[0]);
            for(int d=0; d<3; d++){
                grids_local[lid].left_edge[d]  = space[d] * (double)grid_idx[d];
                grids_local[lid].right_edge[d] = grids_local[lid].left_edge[d] + space[d];
                grids_local[lid].grid_dimensions[d] = 16;
            }
            grids_local[lid].id = gid;
            grids_local[lid].parent_id = -1;
            grids_local[lid].level = 0;
        }

        yt_commit_grids();

        yt_inline("yt_inline");

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();

    return 0;
}

yt load ds

import yt, sys, gc
import numpy as np

yt.enable_parallelism()

def yt_inline(  ):
   ds = yt.frontends.libyt.libytDataset()

def yt_inline_inputArg( field ):
   pass

• Result:

  • num_grids: 2.5e5
  • MPI Size: 4

  

yt ds.print_stats()

import yt, sys, gc
import numpy as np

yt.enable_parallelism()

def yt_inline(  ):
   ds = yt.frontends.libyt.libytDataset()
   ds.print_stats()

def yt_inline_inputArg( field ):
   pass

• Result:

  • num_grids: 2.5e5
  • MPI Size: 4
  • Only root rank executes ds.print_stats().
  • Root rank:

  

  • Non-root rank:

  

Covering Grid No Slabs

import yt, sys, tracemalloc, gc
import numpy as np

data_id_shift = 0 # shifting for naming the data dump
lv = 0
dim_x = 100    # in unit of N_cell for covering grid
dim_y = 100    # in unit of N_cell for covering grid 
dim_z = 100    # in unit of N_cell for covering grid 

box_range = 1e2

center = [ box_range/2.,\
           box_range/2.,\
           box_range/2. ]

yt.enable_parallelism()

def yt_inline(  ):
   ds = yt.frontends.libyt.libytDataset()
   idx = int(str( ds ).split("Fig")[-1]) + data_id_shift
   ad = ds.covering_grid( level=lv, left_edge=[ center[0], center[1], center[2] ], dims=[dim_x,dim_y,dim_z] )
   real_part = ad["Dens"].in_units("code_mass / code_length**3").astype(np.float32)
   if yt.is_root():
      with open("./covering-grid_test_Data_%06d_lv=%d.cbin"%(idx,lv),"wb") as f_RI:
          np.swapaxes(real_part,0,2).tofile(f_RI)

def yt_inline_inputArg( field ):
   pass

• Result:

  • num_grids: 2.5e5
  • MPI Size: 4
  • Extract covering grid dimensions: 100x100x100
  • Every rank gives the same memory consumption result.
  • Compare this result to ds.print_stats(), it looks like covering_grids has the same issue of memory too large in print_stats.

  

[cindytsai]

yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, load gamer hierarchy, and then draw covering_grids

#include "libyt.h"
#include <mpi.h>
#include <iostream>
#include <fstream>

void DerivedFunc(int len, long *list_gid, yt_array *data_array){
    for(int l=0; l<len; l++){
        for(int idx=0; idx<(16*16*16); idx++){ ((float*)data_array[l].data_ptr)[idx] = 1.0; }
    }
}

int main (int argc, char *argv[]){

    /* initialize testing environment. */
    int myrank;
    int nrank;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &nrank );

    yt_param_libyt param_libyt;
    param_libyt.verbose    = YT_VERBOSE_DEBUG;
    param_libyt.script     = "test_covering_grid_no_slabs";
    param_libyt.check_data = false;
    if ( yt_init( argc, argv, &param_libyt ) != YT_SUCCESS ) printf("yt_init failed!\n");

    /* Load Hierarchy */
    std::fstream myfile("./Hierarchy.txt", std::ios_base::in);
    double left_edge[3], right_edge[3];
    int grid_dim[3], level;
    long id, parent_id;

    /* testing API. */
    long iter = 1;
    int  num_grids = 1823232; // TODO
    int  *num_grids_local = new int  [nrank];
    long *acc_grids       = new long [nrank+1];
    if( num_grids % nrank != 0 ){
        for(int r=0; r<nrank; r++){
            if( r == nrank - 1 ) num_grids_local[r] = num_grids - num_grids_local[0] * (nrank - 1);
            else                 num_grids_local[r] = (num_grids / nrank) + 1;
        }
    }
    else{
        for(int r=0; r<nrank; r++){ num_grids_local[r] = num_grids / nrank; }
    }

    for(int r=0; r<nrank+1; r++){
        acc_grids[r] = 0;
        for(int r1=0; r1<r; r1++){
            acc_grids[r] += num_grids_local[r1];
        }
    }

    yt_param_yt param_yt;
    param_yt.frontend                = "gamer";
    param_yt.fig_basename            = "Fig";
    param_yt.length_unit             = 3.0857e21;
    param_yt.mass_unit               = 1.9885e33;
    param_yt.time_unit               = 3.1557e13;

    param_yt.current_time            = 0.0;
    param_yt.dimensionality          = 3;
    param_yt.refine_by               = 2;
    param_yt.num_grids               = num_grids;
    param_yt.num_fields              = 4;   
    param_yt.domain_dimensions[0]    = 576; // TODO
    param_yt.domain_dimensions[1]    = 576; // TODO
    param_yt.domain_dimensions[2]    = 576; // TODO

    for (int d=0; d<3; d++)
    {
        param_yt.domain_left_edge [d] = 0.0;         // TODO
        param_yt.domain_right_edge[d] = 0.15823011;  // TODO
        param_yt.periodicity      [d] = 0;
    }

    param_yt.cosmological_simulation = 0;
    param_yt.current_redshift        = 0.5;
    param_yt.omega_lambda            = 0.7;
    param_yt.omega_matter            = 0.3;
    param_yt.hubble_constant         = 0.7;
    param_yt.num_grids_local         = num_grids_local[myrank];

    for(long t=0; t<iter; t++){
        yt_set_parameter( &param_yt );
        const int mhd = 0;   
        yt_add_user_parameter_int("mhd", 1, &mhd);
        const int srhd = 0;
        yt_add_user_parameter_int("srhd", 1, &srhd);

        yt_field *field_list;
        yt_get_fieldsPtr( &field_list );
        char *field_name[] = {"Real", "Imag", "Dens", "Pote"};
        for(int v=0; v<param_yt.num_fields; v++){
            field_list[v].field_name = field_name[v];
            field_list[v].field_define_type = "derived_func";
            field_list[v].swap_axes = true; // default is true, which means data is represented in [z][y][x] coordinate.
            field_list[v].field_dtype = YT_DOUBLE;
            field_list[v].derived_func = DerivedFunc;
        }

        yt_grid *grids_local;
        yt_get_gridsPtr( &grids_local );

        int lid = 0;
        for(long i=0; i<num_grids; i++){
            // Read file
            myfile >> id;
            myfile >> parent_id;
            myfile >> level;
            myfile >> left_edge[0]  >> left_edge[1]  >> left_edge[2];
            myfile >> right_edge[0] >> right_edge[1] >> right_edge[2];
            myfile >> grid_dim[0]   >> grid_dim[1]   >> grid_dim[2];

            // Write to grids_local in different rank respectively.
            if(acc_grids[myrank] <= i && i < acc_grids[myrank+1]){
                for(int d=0; d<3; d++){
                    grids_local[lid].left_edge[d]  = left_edge[d];
                    grids_local[lid].right_edge[d] = right_edge[d];
                    grids_local[lid].grid_dimensions[d] = grid_dim[d];
                }
                grids_local[lid].id = id;
                grids_local[lid].parent_id = parent_id;
                grids_local[lid].level = level;
                lid += 1;
            }
        }

        yt_commit_grids();

        yt_inline("yt_inline");

        yt_free_gridsPtr();
    }

    /* finalize testing environment. */
    if ( yt_finalize() != YT_SUCCESS ) printf("yt_finalize failed!\n");
    MPI_Finalize();
    delete [] num_grids_local;
    delete [] acc_grids;

    return 0;
}

import yt, sys, tracemalloc, gc
import numpy as np

data_id_shift = 0 # shifting for naming the data dump
lv = 6
dim_x = 1024    # in unit of N_cell for covering grid
dim_y = 1024    # in unit of N_cell for covering grid 
dim_z = 1024    # in unit of N_cell for covering grid 

box_range = 1.58230113e-01                 # in unit of Mpc/h
lv_max_resolution = 4.29226653e-06 # in unit of Mpc/h

center = [ box_range/2.,\
           box_range/2.,\
           box_range/2. ]                # in unit of Mpc/h

covering_grid_left_edge_shift_x = -dim_x/2.        # in unit of N_cell for covering grid, with respective to peak density 
covering_grid_left_edge_shift_y = -dim_y/2.        # in unit of N_cell for covering grid, with respective to peak density
covering_grid_left_edge_shift_z = -dim_z/2.        # in unit of N_cell for covering grid, with respective to peak density

yt.enable_parallelism()

def yt_inline(  ):
   ds = yt.frontends.libyt.libytDataset()
   idx = int(str( ds ).split("Fig")[-1]) + data_id_shift
   ad = ds.covering_grid( level=lv, left_edge=[ center[0] + covering_grid_left_edge_shift_x*lv_max_resolution,\
                                                center[1] + covering_grid_left_edge_shift_y*lv_max_resolution,\
                                                center[2] + covering_grid_left_edge_shift_z*lv_max_resolution ],\
                                                dims=[dim_x,dim_y,dim_z] )
   real_part = ad["Real"].in_units("sqrt(code_mass)/code_length**(3/2)").astype(np.float32)  # convert to code_unit
   if yt.is_root():
      with open("./covering-grid_test_Data_%06d_lv=%d.cbin"%(idx,lv),"wb") as f_RI:
          np.swapaxes(real_part,0,2).tofile(f_RI)

def yt_inline_inputArg( field ):
   pass

Change dim_x/y/z in the above script.

1024x1024x1024

800x800x800

512x512x512

256x256x256

• This is really slow, even much slower than 512x512x512. I don't know if this is caused by machine or what.

128x128x128

64x64x64

[cindytsai]

yt_set_parameter, yt_get_fieldsPtr, yt_get_gridsPtr, yt_commit_grids, load gamer hierarchy, and then slice target covering_grid to smaller unit

(same as above)

import yt
import numpy as np

# Target Covering Grid INFO
lv = 6
dim_x = 1024    # in unit of N_cell for covering grid
dim_y = 1024    # in unit of N_cell for covering grid 
dim_z = 1024    # in unit of N_cell for covering grid 

box_range = 1.58230113e-01         # in unit of Mpc/h
lv_max_resolution = 4.29226653e-06 # in unit of Mpc/h

center = [ box_range/2.,\
           box_range/2.,\
           box_range/2. ]                # in unit of Mpc/h

covering_grid_left_edge_shift_x = -dim_x / 2.        # in unit of N_cell for covering grid, with respective to peak density 
covering_grid_left_edge_shift_y = -dim_y / 2.        # in unit of N_cell for covering grid, with respective to peak density
covering_grid_left_edge_shift_z = -dim_z / 2.        # in unit of N_cell for covering grid, with respective to peak density

target_left_edge = [ center[0] + covering_grid_left_edge_shift_x * lv_max_resolution,\
                     center[1] + covering_grid_left_edge_shift_y * lv_max_resolution,\
                     center[2] + covering_grid_left_edge_shift_z * lv_max_resolution  ]

# Output with 128x128x128 Cubes
cube_dim_x, cube_dim_y, cube_dim_z = 128, 128, 128
stride_x = cube_dim_x * lv_max_resolution
stride_y = cube_dim_y * lv_max_resolution
stride_z = cube_dim_z * lv_max_resolution
cube_filename = "Data_{}_{}_{}_{}.npz"            # Data_{step}_{cube_x}_{cube_y}_{cube_z}.cbin

yt.enable_parallelism()

def yt_inline(  ):
    ds = yt.frontends.libyt.libytDataset()
    data_idx = int(str( ds ).split("Fig")[-1])

    for cube_z in range( int(dim_z / cube_dim_z) ):
        for cube_y in range( int(dim_y / cube_dim_y) ):
            for cube_x in range( int(dim_x / cube_dim_x) ):
                shift = [ cube_x * stride_x, cube_y * stride_y, cube_z * stride_z ]
                ad = ds.covering_grid( level=lv, 
                                       left_edge=[ target_left_edge[0] + shift[0],\
                                                   target_left_edge[1] + shift[1],\
                                                   target_left_edge[2] + shift[2]],
                                       dims=[cube_dim_x, cube_dim_y, cube_dim_z] )
                real_part = ad["Real"].in_units("sqrt(code_mass)/code_length**(3/2)").swapaxes(0, 2) 
                imag_part = ad["Imag"].in_units("sqrt(code_mass)/code_length**(3/2)").swapaxes(0, 2)
                dens_part = ad["Dens"].in_units("code_mass/code_length**3").swapaxes(0, 2)

                if yt.is_root():
                    np.savez(cube_filename.format("%09d" % data_idx, "%d" % cube_x, "%d" % cube_y, "%d" % cube_z), Real=real_part, Imag=imag_part, Dens=dens_part)

def yt_inline_inputArg( field ):
    pass

Slice with different unit by changing cube_dim_x/y/z.

In unit 128x128x128

In unit 512x512x512