Documentation du code de simulation numérique SUNFLUIDH

You find below a python package you can use or include in your own post-processing program to load data computed by Sunlfuidh. The data files are :

• The instantaneous fields (res*.d)
• The statistical files (rst*.d)
• The slice files (slice*.d)
• The space-averaged data file (spav*.d)

Information about these files can be found here (Sunfluidh output files). Information about the file reader is included in the following file.

fct_readsunfluidhdata.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
 
#===========================================================================
#  Author   : Yann Fraigneau, CNRS-LIMSI (yann.fraigneau@limsi.fr)
#  date     : July 2020 - last release : May 2022
#  comments : Function to read the binary data provided by the code SUNFLUIDH
#
#             - ReadSunfluidhDataset : main function to read the binary data of sunfluidh
#             - Check_Dict           : To read a python dictionnary 
#                                      useful to read file information provided by 'ReadSunfluidhDataset'
#             (see an example at the end)
#
#             Keep in mind :
#
#                   - Fields have been computed on a staggered grid (MAC)
#                     (to see some scketches : https://sunfluidh.lisn.upsaclay.fr/doku.php?id=sunfluidh:sunfluidh_tools)
#                            - scalar quantities are defined at the cell-centers
#                            - for snapshot files only the velocity components are located at the cell-faces. For other types of file, they are located at the cell-centers
#                            - domain ends are placed at the cell-faces 
#                            - Cells surrounding the computational domain are "ghost-cells" only used to define the boundary conditions
#                              They can be considered outside of the computational domain
#                   - Coordinates provided in files are associated to the cell-center points
#                   - Velocity components can be easily placed at the cell-centers (as scalar quantities).
#                     For that, set flag_center= True (see the example at the end)
#
#===========================================================================
 
#------------------------------
#--- modules
#------------------------------
 
import numpy as np
import os
import sys
import glob
 
#------------------------------
#--- Global variables
#------------------------------
 
_lst_filetype=['snapshots','slices','statistics','spavslices']
 
#------------------------------
#--- Functions
#------------------------------
 
#---------------------------------------------------------------------------------
#---------------------------------------------------------------------------------
 
def _ExistingFileList (dir_name,dset_type) :
    """
       Check te existence of the access path to dataa, the existence of files related to the type of data
       List the present files
 
       Args IN :
             dset_type (str) : Type of data ("snapshots","slices,"statistics")
             dir_name  (str) : Acces path (where files are stored)
 
       return  : Sorted list of present files
 
    """
 
    #--- Directory exists ?
 
    if os.path.exists(dir_name) :
        print("The path acces '{}' is present".format(dir_name))
    else :
        print("The path acces '{}' is not present".format(dir_name))
        print("--- Program stops ---")
        exit ()
 
    if dset_type == _lst_filetype[0] :
        root="res_*"
    elif dset_type == _lst_filetype[1] :
        root="slices_*"
    elif dset_type == _lst_filetype[2] :
        root="rst_*"
    else :
        print("The selected file type does not exit")
        print("Selected a file type among this list : {}".format(_lst_filetype))
        print("--- Program stops ---")
        exit ()
 
    os.chdir(dir_name)
 
    return sorted(glob.glob(root))
 
#---------------------------------------------------------------------------------
#---------------------------------------------------------------------------------
 
def _BuildDataAccess (dset_type,dir_name,id_subdom,id_time,slice_id,slice_dir,flag_checkexist=False):
    """ Private function used by ' ReadSunfluidhDataset'
 
        Build the data file name and its directory access
 
        Args IN :
             dset_type (str) : Type of data ("snapshots","slices,"statistics")
             dir_name  (str) : Acces path (where files are stored)
             id_subdom (str) : subdomain rank  : "0" or "1" or ... "0" if no MPI domain decomposition
             id_time   (str) : Time ID ( "0" or "1" or ...)
             slice_id  (str) : slice ID   (for slice files only)
             slice_dir (str) : file orientation  (for slice files only)
 
        return : data acces (directory/file name) , data file name
 
    """
 
    if dset_type not in _lst_filetype :
        print(" Define the dataset type to read : \n  - slices \n  - snapshots \n  - statistics \n - \spavslices")
        print(" --- Stop ---")
        exit()
    elif dset_type == _lst_filetype[0] :
        file_name='res_'+id_subdom.zfill(5)+'_'+id_time.zfill(7)+'.d'
    elif dset_type == _lst_filetype[1] :
        if slice_id == '0' :
            print("The slice ID is not set")
            print(" -- Stop --")
            exit()
        if slice_dir == '0' :
            print("The slice orientation is not set")
            exit()
        if id_subdom == '0' :
            f_name=dir_name+'/'+'slice_'+slice_id.zfill(2)+'_'+slice_dir.zfill(1)+'_'+id_time.zfill(7)+'.d'
            if os.path.exists(f_name) :
                file_name='slice_'+slice_id.zfill(2)+'_'+slice_dir.zfill(1)+'_'+id_time.zfill(7)+'.d'
            else :
                file_name='slice_'+slice_id.zfill(2)+'_'+slice_dir.zfill(1)+'_'+id_subdom.zfill(5)+'_'+id_time.zfill(7)+'.d'
        else :
            file_name='slice_'+slice_id.zfill(2)+'_'+slice_dir.zfill(1)+'_'+id_subdom.zfill(5)+'_'+id_time.zfill(7)+'.d'
    elif dset_type == _lst_filetype[2] :
        file_name='rst_'+id_subdom.zfill(5)+'_'+id_time.zfill(7)+'.d'
    elif dset_type == _lst_filetype[3] :
        file_name='spav_slice_'+slice_dir.zfill(1)+'_'+id_time.zfill(7)+'.d'
 
    data_access=dir_name+'/'+file_name
 
    flag= True
    if not os.path.exists(data_access) and flag_checkexist:
        print("The data file : ",data_access," does not exist. It is ignored.")
        flag= False
        exit()
 
    return data_access,file_name
 
#---------------------------------------------------------------------------------
#---------------------------------------------------------------------------------
 
def _CenteredVelocity (fieldname,field) :
 
    """ Computation of velocity components at the cell-centre 
 
        IN :
              fieldname (list of str) : list of field names
              field   (numpy array)   : fields
 
        centered velocity fields are stored in the array "field" (overwriting previous velocity fields)
 
    """
 
    ib=0 ; ie=field.shape[0]
    jb=0 ; je=field.shape[1]
    kb=0 ; ke=field.shape[2]
 
    for i,n in enumerate(fieldname):
        if n == "U" :
            field[ib+1:ie,jb:je,kb:ke,i]=(field[ib:ie-1,jb:je,kb:ke,i] \
                                              +field[ib+1:ie,jb:je,kb:ke,i])*5.e-01
        elif n == "V" :
            field[ib:ie,jb+1:je,kb:ke,i]=(field[ib:ie,jb:je-1,kb:ke,i] \
                                              +field[ib:ie,jb+1:je,kb:ke,i])*5.e-01
        elif n == "W" and ke > 1 :
            field[ib:ie,jb:je,kb+1:ke,i]=(field[ib:ie,jb:je,kb:ke-1,i] \
                                              +field[ib:ie,jb:je,kb+1:ke,i])*5.e-01
 
    return
 
#---------------------------------------------------------------------------------
#---------------------------------------------------------------------------------
 
def ReadSunfluidhDataset (dset_type="snapshots",dir_name=".",id_subdom=0,id_time=-1,slice_id=1,slice_dir=1,selected_fields=[],flag_centered_vel=False) :
    """
        Read the full binary dataset from the suitable file
 
      *  Synopsis :
 
          file_info,coordinates,field= ReadSunfluidhDataset (dset_type= (str), dir_name= (str), id_subdom= (str), id_time(str), .... )
 
      * Keyword arguments (IN) :
 
        - dset_type (str, optional) : type of dataset considered ("snapshots", "slices", "statistics"). Default value= "snapshots"
        - dir_name  (str)           : directory where files are located
        - id_subdom (int, optional) : ID number of the split data file (MPI computation with domain decomposition)
                                      The default value is zero
        - id_time   (int,optional)  : ID number associated to timestep of file recording. If id_time is < 0 -- or omitted--, 
                                                                                          is selected the last file created 
                                                                                          (with the greatest value of id_time)
        - slice_id  (int,optional)  : ID number of slice dataset. Useful for 'slice' files only. Default value is zero
        - slice_dir (int,optional)  : normal direction of slice dataset  (1= I,2= J or 3=K). Useful for 'slice' file only.  Default value is zero
        - selected_fields (list of str, optional) : list of selected fields 
                                                    (empty list by default : all fields are selected)
                                                    see the sunfluidh wiki (sunfluidh.limsi.fr) for more details on field names
        - flag_centered_vel (boolean) : if True, the velocity field components are interpolated at at the cell-center (as scalar qunatities)
                                        The default value is False
 
      * Return (OUT) :
 
        - fil_info (dict)                     : relevant data on the data file (grid size, fields name, ...)
                                                 contents of this dictionnary can be checked with the function "check_dict"
        - coordinates (list of 1D numpy array) : coordinates
                                                 coordinates[0]  : cell-center coordinates along the I-direction
                                                 coordinates[1]  : cell-center coordinates along the J-direction
                                                 coordinates[2]  : cell-center coordinates along the K-direction
        - field  (4D numpy arrays)             : arrays of fields
                                                 field (i,j,k,n) : with i,j,k : grid indices
                                                                            n : field ID ranked as the field names in file_info['field names'] 
 
      """
 
    file_info={}
 
    #---------------------------------------------------------------------------------
    #---- Data access
    #---------------------------------------------------------------------------------
 
    if not isinstance(id_time,int) :
       print("Variable 'id_time' must be a integer")
       print("--- Program stops ---")
       exit ()
    if not isinstance(id_subdom,int) :
       print("Variable 'id_subdom' must be a integer")
       print("--- Program stops ---")
       exit ()
    if not isinstance(slice_id,int) :
       print("Variable 'slice_id' must be a integer")
       print("--- Program stops ---")
       exit ()
    if not isinstance(slice_dir,int) :
       print("Variable 'slice_dir' must be a integer")
       print("--- Program stops ---")
       exit ()
 
    if id_time >= 0 :
        data_access,filename=_BuildDataAccess (dset_type,dir_name,str(id_subdom),str(id_time),str(slice_id),str(slice_dir))
    else :
        lst_file=_ExistingFileList (dir_name=dir_name,dset_type= dset_type)
        filename=lst_file[-1]
        data_access=dir_name+'/'+filename
        print()
        print("Automatic selection of the last file created {}".format(data_access))
        print()
 
    file_info['file name']=filename
    file_info['dataset type']=dset_type
 
    #---------------------------------------------------------------------------------
    # Structure and size of each data type
    #---------------------------------------------------------------------------------
 
    dt_int32=np.dtype('>i4') 
    dt_int64=np.dtype('>i8')
    dt_flt32=np.dtype('>f4')
    dt_flt64=np.dtype('>f8')
    dt_str06=np.dtype('a6')
 
    with open(data_access,"rb") as f :
        #------------------------------------------------------------------------------------------------
        #--- Extract The three 1st integer 
        #    (Binary file version, number of fields, number of elements in the file header)
        #------------------------------------------------------------------------------------------------
 
        data=np.fromfile(f,dtype=dt_int32,count=3,offset=0)
        file_info['file version']=np.abs(data[0])
        file_info['number of fields']=data[1]
        #file_info['number of useless elmt']=data[2]
 
        #------------------------------------------------------------------------------------------------
        #--- Extract the data location (useless data at present)
        #------------------------------------------------------------------------------------------------
 
        data=np.fromfile(f,dtype=dt_int64,count=data[1]+data[2],offset=0)
        #print("useless = ",data)
 
        #------------------------------------------------------------------------------------------------
        #--- Extract the mesh size
        #------------------------------------------------------------------------------------------------
 
        data=np.fromfile(f,dtype=dt_int32,count=3,offset=0)
        file_info['mesh size']=data[0:3]
 
        #------------------------------------------------------------------------------------------------
        #--- Extract some simulation features
        #    single or double real precision, cylindrical geometry, centered velocity or not 
        #------------------------------------------------------------------------------------------------
 
        data=np.fromfile(f,dtype=dt_int32,count=3,offset=0)
 
        if dset_type == _lst_filetype[1] : data[2] = 1 #--- correction for slice (already centered 17/01/2023)
 
        file_info['real precision']=data[0]
        if file_info['real precision'] == 1 : #--- Single float precision
            dt_flt= dt_flt32
        else:    #--- Double float precision
            dt_flt= dt_flt64
 
        file_info['geometry type']=data[1]
        if data[2] == 1 :
            file_info['centered velocity']= True
        else:
            file_info['centered velocity']= False
 
        #------------------------------------------------------------------------------------------------
        #--- Extract list of field name
        #------------------------------------------------------------------------------------------------
 
        if file_info['geometry type'] == 0:
            file_info['coordinate names']= ["X","Y","Z"]
        elif file_info['geometry type'] == 1:
            file_info['coordinate names']= ["R","Theta","Z"]
        elif file_info['geometry type'] == 2:
            file_info['coordinate names']= ["Z","R","Theta"]
        elif file_info['geometry type'] == 3:
            file_info['coordinate names']= ["Theta","Z","R"]
        data=np.fromfile(f,dtype=dt_str06,count=file_info['number of fields'],offset=0)
        fieldname=[ data[i].decode('utf-8').replace(" ","") for i in range(len(data))]
 
        #------------------------------------------------------------------------------------------------
        #--- Extract Time
        #------------------------------------------------------------------------------------------------
 
        data=np.fromfile(f,dtype=dt_flt,count=1,offset=0)
        file_info['time']= data[0]
 
        #------------------------------------------------------------------------------------------------
        #--- Extract coordinates
        #------------------------------------------------------------------------------------------------
 
        if file_info['mesh size'][2] == 1 :
            file_info['geometry dimension']=2
        else:
            file_info['geometry dimension']=3
 
        coordinates=[0,0,0]
        for n in range(0,3):
            coordinates[n]=np.fromfile(f,dtype=dt_flt,count=file_info['mesh size'][n],offset=0)
            #print(file_info['coordinate names'][n] ,' start coord : ',coordinates[n][0],\
            #                'end coord : ',coordinates[n][-1],' size : ',np.size(coordinates[n]))
 
        #------------------------------------------------------------------------------------------------
        #--- Extract specific data related to the dataset type
        #------------------------------------------------------------------------------------------------
 
        if dset_type == _lst_filetype[1] :
            data=np.fromfile(f,dtype=dt_int32,count=2,offset=0)
            file_info['slice direction' ]= data[0]
            file_info['slice index' ]= data[1]
 
        elif dset_type == _lst_filetype[2] :
            data=np.fromfile(f,dtype=dt_int32,count=3,offset=0)
            file_info['average type' ]= data[0]
            file_info['sample' ]= data[1]
            file_info['range number' ]= data[2]
            data=np.fromfile(f,dtype=dt_flt,count=1,offset=0)
            file_info['statistic time' ]= data[0]
 
        #------------------------------------------------------------------------------------------------
        #---             FIELDS
        #------------------------------------------------------------------------------------------------
 
        nmax=file_info['mesh size'][0]*file_info['mesh size'][1]*file_info['mesh size'][2] \
                                      *file_info['number of fields']
        data=np.fromfile(f,dtype=dt_flt,count=nmax,offset=0)
        field=data.reshape(file_info['mesh size'][0],file_info['mesh size'][1],file_info['mesh size'][2],file_info['number of fields'],order='F')
 
 
        #--- Collect fields selected by the user
 
        if len(selected_fields) != 0 :
            fieldname,field=  SunfluidhData.FieldSelection (fulllist=fieldname,selectlist=selected_fields,fields=field)
            file_info['number of fields']= len(fieldname)
 
        #--- Velocity location at the cell-centre
 
        if flag_centered_vel and not file_info['centered velocity'] :
            _CenteredVelocity (fieldname=fieldname,field=field)
            file_info['centered velocity']= True
 
        file_info['field names']= fieldname
 
        #print("NEW READ")
        #print("SYS byte order     ", sys.byteorder)
        #print("CHECK name       : ",field.dtype.name)
        #print("CHECK byte order : ",field.dtype.byteorder)
 
        #for n in range(file_info['number of fields']) :
        #    print ("var {}, min {}, max {}".format(fieldname[n],np.min(field[:,:,:,n]),np.max(field[:,:,:,n])) )
 
 
    #exit ()
 
    return file_info,coordinates,field
 
#---------------------------------------------------------------------------------
#---------------------------------------------------------------------------------
 
def Check_Dict (title= "", dic={}) :
    """ check dictionnary contents """
 
    print("\nCheck dictionnary : ", title, "\n")
    for key,value in dic.items():
        print('{0:30} : {1:}'.format(key,value))
 
#---------------------------------------------------------------------------------
#             TEST
#---------------------------------------------------------------------------------
 
if __name__ == "__main__" :
 
    #--------------------------------------------------------------------------------------------------
    #--- Parameters to build the file name to read (arguments for the function "ReadSunfluidhDataset")
    #    See here for information about files created by Sunfluidh 
    #          (https://sunfluidh.lisn.upsaclay.fr/doku.php?id=sunfluidh:sunfluidh_output)
    #--------------------------------------------------------------------------------------------------
 
    dirname="/Users/yann/Projets/WORK_SUNFLUIDH/DATABASE_LUCOR/CAS_SQUARECYL/SNAPSHOTS"         #--- Path to the directory where files are stored
                        #    for example : dirname="/Users/path_to_datset"
                        #    This parameter is optional and it could be omitted. 
                        #    In this case, default value="." and the python script must be run in the directory where files are stored
 
    dsetype="snapshots" #--- Type of file : "snapshots"  --> instantaneous 3D or 2D fields.
                        #                                    Name :  res_idsdom_idtime.d , 
                        #                   "slices",    --> this file corresponds to a plane extracted from a 3D field (3D simulation only)
                        #                                    Name : slice_sliceid_slicedir_idsdom_idtime.d
                        #                   "statistics" --> (3D or 2D) averaged fields  
                        #                                    Name : rst_idsdom_idtime.d 
 
    idsdom = 0          #--- subdomain ID  (if the computation has been performed in MPI domain decomposition, 
                        #                   fields are split in space. id_sdom is the ID value related to the MPI subdomain ID) 
                        #
                        #  This parameter is optional and it could be omitted. In this case, default value=0 
                        #                                                      (unic domain ID for simulation without MPI domain decomposition)
 
    idtime=  250          #--- time ID of the file (time ID of file). 
                        #    This parameter is optional and it could be omitted.
 
    sliceid=0           #--- For slice files only : define the rank ID of the extracted plan as it has been defined in the input data file of sunfluidh
                        #    This parameter is optional and it could be omitted (default value=1).
 
    slicedir=1          #--- For slice files only : orientation of the extracted plan defined by its normal vector (1 : I-direction, 2 : J-direction, 3: K-direction)
                        #    This parameter is optional and it could be omitted (default value=1).
 
    lst_fields=[]       #--- Allow you to select the fields to read in the file. When the list is empty, all fields are read by default
                        #    ex : lst_fields ['T'] --> select the temperature only (if it is present)
                        #     see the sunfluidh wiki (sunfluidh.limsi.fr) for more details on field names
                        #    This parameter is optional and it could be omitted (default value=[] empty list).
 
    flag_center= True   #--- The velocity field components are interpolated at the cell-center (as the scalar quantities)
                        #    This parameter is optional and it could be omitted (default value=False).
 
    #--- Read file
 
    #.....  Put here loops on the variables idtime (or/and idsdom if needed) according to files to read
 
    file_info,coordinates,field= ReadSunfluidhDataset (dset_type=dsetype,dir_name=dirname,id_subdom=idsdom,id_time=idtime, \
                                       slice_id=sliceid,slice_dir=slicedir,selected_fields=lst_fields, flag_centered_vel= flag_center)
 
    # OUTPUTS : 
    #--- file_info   : dictionnary       : It contains information describing the data : grid size , fields ...
    #--- coordinates : list of 1D arrays : each element of the list contains the coordinates at the center of cell in a given direction
    #                                    : coordinates[0] : in the x-direction
    #                                    : coordinates[1] : in the y-direction
    #                                    : coordinates[2] : in the y-direction
    #--- field       : 4D array field(i,j,k,n)
    #                  i,j,k refers the grid indices. For 2D fields, k=0
    #                  n refers the field rank (in the same order as field names expressed in file_info['field names']
    #                  BEWARE : the array size about spatial indices (i,j,k) includes ghost-cells 
    #                           Ghost-Cells surround the discretized computational(sub)domain, 
    #                           they are used to treat the boundary conditions during the computation
    #                           for 'snapshot' or 'statistics' files, the ghost-cells are the first and last cells in each direction
    #                           for 'slice' files, the ghost-cells are only the last cells in each direction
    #                           These ghost-cells are usually not important except in the case of periodic direction
    #                           you can remove these ghost-cells for data processing 
    #                           be careful in cases where these ghost-cell can be relevant : i.e. periodic domains 
 
    #--- check the dictionnary  file_info which contains information about the file
 
    Check_Dict (title= "Information collected in dictionnary 'file_info' ", dic=file_info)