Data Mapping Issue

Hi,

I have written the following Fortran code that uses the cgns library. The objective of this code is to add a new solution field to input.cgns and swap the data of two fields. For swapping, it correctly identifies the fields and performs the data swapping. However, the mapping is not done correctly, even though the maximum and minimum of the data are correct. Do you have any idea why I am facing this problem?

program manipulate_cgns_2d
  use cgns

  implicit none

  !> Single precision real numbers, 6 digits, range 10⁻³⁷ to 10³⁷-1; 32 bits
  integer, parameter :: sp = selected_real_kind(6, 37)
  !> Double precision real numbers, 15 digits, range 10⁻³⁰⁷ to 10³⁰⁷-1; 64 bits
  integer, parameter :: dp = selected_real_kind(15, 307)
  !> Quadruple precision real numbers, 33 digits, range 10⁻⁴⁹³¹ to 10⁴⁹³¹-1; 128 bits
  integer, parameter :: qp = selected_real_kind(33, 4931)

  !> Char length for integers, range -2⁷ to 2⁷-1; 8 bits
  integer, parameter :: i1 = selected_int_kind(2)
  !> Short length for integers, range -2¹⁵ to 2¹⁵-1; 16 bits
  integer, parameter :: i2 = selected_int_kind(4)
  !> Length of default integers, range -2³¹ to 2³¹-1; 32 bits
  integer, parameter :: i4 = selected_int_kind(9)
  !> Long length for integers, range -2⁶³ to 2⁶³-1; 64 bits
  integer, parameter :: i8 = selected_int_kind(18)
  
  integer(i8), parameter           :: Ndim = 2
  integer(cgsize_t)                :: isize(Ndim*3) ! (vertex size: isize(1:Ndim,1), cell size: isize(1:Ndim,2), boundary vertex size: isize(1:Ndim,3) is 0 for structured mesh) 

  ! Input and output file names
  character(len=128) :: input_file  = "input.cgns"
  character(len=128) :: output_file = "output.cgns"
  ! character(len=128) :: NewVariable = "DummyVariable"
  character(len=256) :: cp_command
  
  integer                       :: idx_1, idx_2
  character(len=128), parameter :: str_1 = 'X_A', str_2 = 'X_B' 
  character(len=128)            :: datatype_str_1, datatype_str_2 
  double precision, allocatable :: data_str_1(:), data_str_2(:)
  integer                       :: nrows, ncols
  
  integer :: ier, file_number, nbases, nzones, base, zone, zonetype
  integer :: discr
  integer :: IndexDim, CellDim, PhysDim

  character(len=32) :: basename, zonename
  integer(cgsize_t)  :: dim_vals(12)
  real, allocatable  :: data(:)
  integer            :: pos, cmp = 0, i, j, k

  character(len=32) :: name
  integer           :: iarray, narrays, datatype

  real(sp) :: data_single(100000)
  real(dp) :: data_double(100000)  

  ! Construct the cp command
  cp_command = 'cp ' // trim(input_file) // ' ' // trim(output_file)
  ! Execute the cp command
  call execute_command_line(cp_command)

  ! *** open file
  call cg_open_f(output_file, CG_MODE_MODIFY, file_number, ier)
  if (ier .ne. CG_OK) call cg_error_exit_f
  write(*,*) "Opening CGNS file."
  write(*,*) ""

  ! Get the number of bases
  call cg_nbases_f(file_number, nbases, ier)
  if (ier .eq. ERROR) call cg_error_exit_f

  do base = 1, nbases

     ! Get the name and dimensionality of the base
     call cg_base_read_f(file_number, base, basename, CellDim, PhysDim, ier)
     if (ier .eq. ERROR) call cg_error_exit_f
     write(*,*) "Getting the name and dimensionality of the base."
     write(*,*) 'BaseName = "', trim(basename), '"'
     write(*,'(A,I0)') ' Cell dimension = ', CellDim
     write(*,*) ""

     ! Get the number of zones
     call cg_nzones_f(file_number, base, nzones, ier)
     if (ier .eq. ERROR) call cg_error_exit_f
     write(*,*) "Getting the number of zones and their sizes."
     write(*,*) ""

     do zone = 1, nzones

        call cg_zone_read_f(file_number, base, zone, zonename, isize, ier)
        if (ier .eq. ERROR) call cg_error_exit_f
        write(*,*) 'Name: ', '"', trim(zonename), '"'
        write(*,*) 'isize = ', isize(:)

        call cg_zone_type_f(file_number, base, zone, zonetype, ier)
        if (ier .eq. ERROR) call cg_error_exit_f
        write(*,*) 'Zone type is ', ZoneTypeName(zonetype)

        if (zonetype.eq.Structured) then
           IndexDim=CellDim
        else
           IndexDim=1
        endif
        write(*,'(A,I0)') ' IndexDimension = ', IndexDim
        write(*,*) ""

        ! Navigate to the 'FlowSolution_t' node under the 'Zone_t' node
        call cg_goto_f(file_number, base, ier, 'Zone_t', zone, 'FlowSolution_t', 1, 'end')
        if (ier .ne. CG_OK) call cg_error_exit_f
        write(*,*) "Navigating 'FlowSolution_t' node under the", zone       
        write(*,*) ""

        ! Create and write discrete data associated with a specific CGNS zone
        call cg_discrete_write_f(file_number, base, zone, 'discrete#1', discr, ier)
        if (ier .eq. ERROR) call cg_error_exit_f
        write(*,*) "Creating and writing discrete data associated with the", zone       
        write(*,*) ""

        ! Create a new data array (assuming a scalar real variable named "NewVariable")    
        do j=1, isize(2)
           do i=1, isize(1)
              cmp = cmp + 1
           enddo
        enddo

        print*, 'cmp = ', cmp
        allocate(data(cmp))

        do j=1, isize(2)
           do i=1, isize(1)
              pos = i + (j-1)*isize(1) 
              data(pos) = pos	! * make up some dummy data
              ! print*, 'data(', pos, ')=', data(pos) 
           enddo
        enddo

        ! Write a data array (data) to a CGNS file for a given variable (NewVariable) with specified data type, dimension, and size
        ! call cg_array_write_f(NewVariable, RealSingle, IndexDim, isize, data, ier)
        ! if (ier .eq. ERROR) call cg_error_exit_f
        ! write(*,*) "Writing a data array to the CGNS file for a given variable named ", '"', trim(NewVariable), '"', ' with specified data type, dimension, and size'       
        ! write(*,*) ""


        call cg_narrays_f(narrays, ier)
        if (ier .eq. ERROR) call cg_error_exit_f
        write(*,*) '  ReferenceState_t contains ', narrays,' array(s)'

        do iarray=1, narrays

           call cg_array_info_f(iarray, name, datatype, IndexDim, dim_vals, ier)
           if (ier .eq. ERROR) call cg_error_exit_f

           write(*,*)  '   DataArray #',iarray,':'
           write(*,*)  '    Name =',name
           write(*,*)  '    Datatype=',DataTypeName(datatype)

           write(*,*)'    Data:'
           if (datatype .eq. RealSingle) then
              ! call cg_array_read_f(iarray, data_single, ier)
              ! if (ier .eq. ERROR) call cg_error_exit_f
              ! write(*,*) data_single(1)
              write(*,*) 'datatype is RealSingle. The cg_array_read_f does not support this type.'              
           elseif (datatype .eq. RealDouble) then
              call cg_array_read_f(iarray, data_double, ier)
              if (ier .eq. ERROR) call cg_error_exit_f
              write(*,*) data_double(1)
           endif
           
           nrows = isize(1)
           ncols = isize(2)
           if (name==str_1) then
              idx_1 = iarray
              datatype_str_1 = DataTypeName(datatype)
              cmp = size(data_double)
              allocate(data_str_1(cmp))
              data_str_1 = 0.0d0
              data_str_1 = data_double

           else if (name==str_2) then
              idx_2 = iarray
              datatype_str_2 = DataTypeName(datatype)
              cmp = size(data_double)
              allocate(data_str_2(cmp))
              data_str_2 = 0.0d0
              data_str_2 = data_double  
           else
           endif
              
        enddo

        
        ! Swap the arrays              
        if ( (idx_1.ne.0).and.(idx_2.ne.0) ) then
           write(*,*) ""
           write(*,*) 'Swapping the arrays: ' // trim(str_1) // ' and ' // trim(str_2)
           write(*,*) trim(str_1), ' index number: ', idx_1
           write(*,*) trim(str_2), ' index number: ', idx_2

           if (datatype_str_1 == datatype_str_2 ) then 
              ! Use predefined constants for datatypes
              if (datatype_str_1 == 'RealSingle') then
                 call cg_array_write_f(trim(str_1), RealSingle, IndexDim, isize, data_str_2, ier)
                 if (ier .eq. ERROR) call cg_error_exit_f
                 call cg_array_write_f(trim(str_2), RealSingle, IndexDim, isize, data_str_1, ier)
                 if (ier .eq. ERROR) call cg_error_exit_f
              elseif (datatype_str_1 == 'RealDouble') then
                 call cg_array_write_f(trim(str_1), RealDouble, IndexDim, isize, data_str_2, ier)
                 if (ier .eq. ERROR) call cg_error_exit_f
                 call cg_array_write_f(trim(str_2), RealDouble, IndexDim, isize, data_str_1, ier)
                 if (ier .eq. ERROR) call cg_error_exit_f
              else
                 write(*,*) 'Unsupported datatype:', datatype_str_1
                 stop
              endif
           else
              write(*,*) 'String 1 data type does not match the one of string 2'
              stop
           endif
        endif



     enddo
  enddo
  

  ! Close the CGNS file
  call cg_close_f(file_number, ier)
  if (ier .ne. CG_OK) call cg_error_exit_f
  write(*,*) ""
  write(*,*) "Closing CGNS file."
  write(*,*) ""

end program manipulate_cgns_2d

Thanks in advance for your help.

1 Like

Hi @Reiziger,

welcome to Discourse.

Can you tell us roughly at which statement do things go wrong and what have you tried so far? Is there any easy way to install cgns for others to test it?

I found a few other things a bit odd-looking, such as the double assignment and empty else block:

              data_str_2 = 0.0d0
              data_str_2 = data_double  
           else
           endif

Or the loop-based counting:

        do j=1, isize(2)
           do i=1, isize(1)
              cmp = cmp + 1
           enddo
        enddo

which could be replaced with:

cmp = cmp + product(isize(1:2))
! or
cmp = cmp + isize(1)*isize(2)