Merge pull request #1493 from ExtremeFLOW/martin/map12d

Add map_2d

contrib/average_field_in_space/average_field_in_space.f90

@@ -9,18 +9,22 @@ program average_field_in_space
   real(kind=rp) :: start_time, el_h, el_dim(3,3), domain_height
   real(kind=rp), allocatable :: temp_el(:,:,:)
   type(fld_file_data_t) :: field_data
+  type(fld_file_data_t) :: output_data
   type(coef_t) :: coef
-  type(dofmap_t) :: dof
+  type(dofmap_t), target :: dof
   type(space_t) :: Xh
   type(mesh_t) :: msh
   type(gs_t) :: gs_h
   type(map_1d_t) :: map_1d
-  type(field_t), pointer :: u, avg_u, old_u, el_heights
-  type(vector_ptr_t), allocatable :: fields(:)
+  type(map_2d_t) :: map_2d
+  type(vector_ptr_t), allocatable :: fields(:), fields2d(:)
   type(matrix_t) :: avg_matrix
   type(vector_t) :: volume_per_gll_lvl
-  integer :: argc, i, n, lx, j, e, n_levels, dir, ierr, n_1d, tstep
+  integer :: argc, i, n, lx, j, e, n_levels, dir, ierr, n_1d, tstep, k
   logical :: avg_to_1d = .false.
+  integer :: nelv_2d, glb_nelv_2d, offset_el_2d, lxy, n_2d
+  integer, allocatable :: idx_2d(:)
+  real(kind=rp), pointer, dimension(:,:,:,:) :: x_ptr, y_ptr
   real(kind=rp) :: coord
 
   argc = command_argument_count()
@@ -97,14 +101,6 @@ program average_field_in_space
   call gs_h%init(dof)
   call coef%init(gs_h)
 
-  call neko_field_registry%add_field(dof, 'u')
-  u => neko_field_registry%get_field('u')
-  call neko_field_registry%add_field(dof, 'avg_u')
-  avg_u => neko_field_registry%get_field('avg_u')
-  call neko_field_registry%add_field(dof, 'old_u')
-  old_u => neko_field_registry%get_field('old_u')
-  call neko_field_registry%add_field(dof, 'el_heights')
-  el_heights => neko_field_registry%get_field('el_heights')
   ! 1 corresponds to x, 2 to y, 3 to z
   if (trim(hom_dir) .eq. 'x') then
      dir = 1
@@ -127,193 +123,30 @@ program average_field_in_space
   else
      call neko_error('homogenous direction not supported')
   end if
-  call map_1d%init(dof, gs_h,  dir, 1e-7_rp)
-  n_levels = map_1d%n_el_lvls
-  n = u%dof%size()
-
+  if (avg_to_1d) then
+     call map_1d%init(coef, dir, 1e-7_rp)
+  else
+     call map_2d%init(coef, dir, 1e-7_rp)
+  end if
+
   !allocate array with pointers to all vectors in the file
-  allocate(fields(field_data%size()))
-  ! Compute average in two direction directly and store in a csv file
+
   output_file = file_t(trim(output_fname))
   do tstep = 0, field_data%meta_nsamples-1
+    if (pe_rank .eq. 0) write(*,*) 'Averaging field:', tstep
      if (tstep .gt. 0) call field_file%read(field_data)
-     call field_data%get_list(fields,field_data%size())
-     if (pe_rank .eq. 0) write(*,*) 'Averaging field:', tstep
      if (avg_to_1d) then
-        n_1d = n_levels*Xh%lx
-        call avg_matrix%init(n_1d,field_data%size()+1)
-        call volume_per_gll_lvl%init(n_1d)
-        do i = 1, n
-           volume_per_gll_lvl%x(map_1d%pt_lvl(i,1,1,1)) = &
-           volume_per_gll_lvl%x(map_1d%pt_lvl(i,1,1,1)) + coef%B(i,1,1,1)
-        end do
-        call MPI_Allreduce(MPI_IN_PLACE,volume_per_gll_lvl%x, n_1d, &
-           MPI_REAL_PRECISION, MPI_SUM, NEKO_COMM, ierr)
-        !ugly way of getting coordinates, computes average
-        do i = 1, n
-           if (dir .eq. 1) coord = dof%x(i,1,1,1)
-           if (dir .eq. 2) coord = dof%y(i,1,1,1)
-           if (dir .eq. 3) coord = dof%z(i,1,1,1)
-           avg_matrix%x(map_1d%pt_lvl(i,1,1,1),1) = &
-           avg_matrix%x(map_1d%pt_lvl(i,1,1,1),1) + coord*coef%B(i,1,1,1) &
-           /volume_per_gll_lvl%x(map_1d%pt_lvl(i,1,1,1))
-        end do
-        do j = 2, field_data%size()+1
-           do i = 1, n
-              avg_matrix%x(map_1d%pt_lvl(i,1,1,1),j) = &
-              avg_matrix%x(map_1d%pt_lvl(i,1,1,1),j) + fields(j-1)%ptr%x(i)*coef%B(i,1,1,1) &
-              /volume_per_gll_lvl%x(map_1d%pt_lvl(i,1,1,1))
-           end do
-        end do 
-        call MPI_Allreduce(MPI_IN_PLACE,avg_matrix%x, (field_data%size()+1)*n_1d, &
-           MPI_REAL_PRECISION, MPI_SUM, NEKO_COMM, ierr)
-
+        call map_1d%average_planes(avg_matrix, fields)
         call output_file%write(avg_matrix,field_data%time)
-     ! Compute averages in 1 direction and store in a 3d field (lots of redundant data, sorry)
-     ! Should output a 2d field in principle
+        ! Compute averages in 1 direction and store in a 3d field (lots of redundant data, sorry)
+        ! Should output a 2d field in principle
      else
-        do i = 1, msh%nelv
-           !find height in hom-dir
-           !direction in local coords (r,s,t) that is hom is stored in map_1d%dir_el
-           !set element to height
-           !we assume elements are stacked on eachother...
-           el_dim(1,:) = abs(msh%elements(i)%e%pts(1)%p%x-msh%elements(i)%e%pts(2)%p%x)
-           el_dim(2,:) = abs(msh%elements(i)%e%pts(1)%p%x-msh%elements(i)%e%pts(3)%p%x)
-           el_dim(3,:) = abs(msh%elements(i)%e%pts(1)%p%x-msh%elements(i)%e%pts(5)%p%x)
-           ! 1 corresponds to x, 2 to y, 3 to z
-           el_h = el_dim(map_1d%dir_el(i),dir)
-           el_heights%x(:,:,:,i) = el_h
-        end do
-
-        call copy(u%x,el_heights%x,n)
-        call copy(old_u%x,el_heights%x,n)
-        call copy(avg_u%x,el_heights%x,n)
-        call perform_global_summation(u, avg_u, old_u, n_levels, &
-             map_1d%dir_el,gs_h, coef%mult, msh%nelv, lx)
-        domain_height = u%x(1,1,1,1)
-
-
-        do i = 1, field_data%size()
-           call copy(old_u%x,fields(i)%ptr%x,n)
-           call perform_local_summation(u,old_u, el_heights, domain_height, &
-                map_1d%dir_el, coef, msh%nelv, lx)
-           call copy(old_u%x,u%x,n)
-           call copy(avg_u%x,u%x,n)
-           call perform_global_summation(u, avg_u, old_u, n_levels, &
-                map_1d%dir_el,gs_h, coef%mult, msh%nelv, lx)
-           call copy(fields(i)%ptr%x,u%x,n)
-        end do 
-
-        call output_file%write(field_data,field_data%time)
+        call map_2d%average(output_data,field_data)
+        call output_file%write(output_data,field_data%time)
      end if
   end do
   if (pe_rank .eq. 0) write(*,*) 'Done'
-
   call neko_finalize
 
 end program average_field_in_space
 
-subroutine perform_global_summation(u, avg_u, old_u, n_levels, hom_dir_el, gs_h, mult, nelv,lx)
-  use neko
-  implicit none
-  type(field_t), intent(inout) :: u, avg_u, old_u
-  type(gs_t), intent(inout) :: gs_h
-  integer, intent(in) :: n_levels, nelv, lx
-  integer, intent(in) :: hom_dir_el(nelv)
-  real(kind=rp), intent(in) :: mult(nelv*lx**3)
-  real(kind=rp) :: temp_el(lx,lx,lx)
-  integer :: n, i, j, e, k
-
-  n = u%dof%size()
-
-  do i = 1, n_levels-1
-     !compute average
-     if (NEKO_BCKND_DEVICE .eq. 1) &
-          call device_memcpy(u%x, u%x_d, n, &
-                             HOST_TO_DEVICE, sync=.false.)
-     call gs_h%op(u,GS_OP_ADD)
-     if (NEKO_BCKND_DEVICE .eq. 1) &
-          call device_memcpy(u%x, u%x_d, n, DEVICE_TO_HOST, sync=.true.)
-     call col2(u%x,mult,n)
-     do e = 1, nelv
-        temp_el = 2.0*u%x(:,:,:,e)-old_u%x(:,:,:,e)
-        if (hom_dir_el(e) .eq. 1) then
-           u%x(1,:,:,e) = temp_el(lx,:,:)
-           avg_u%x(1,:,:,e) = avg_u%x(1,:,:,e)+temp_el(1,:,:)
-           u%x(lx,:,:,e) = temp_el(1,:,:)
-        else if (hom_dir_el(e) .eq. 2) then
-           u%x(:,1,:,e) = temp_el(:,lx,:)
-           avg_u%x(:,1,:,e) = avg_u%x(:,1,:,e)+temp_el(:,1,:)
-           u%x(:,lx,:,e) = temp_el(:,1,:)
-        else if (hom_dir_el(e) .eq. 3) then
-           u%x(:,:,1,e) = temp_el(:,:,lx)
-           avg_u%x(:,:,1,e) = avg_u%x(:,:,1,e)+temp_el(:,:,1)
-           u%x(:,:,lx,e) = temp_el(:,:,1)
-        end if
-        old_u%x(:,:,:,e) = u%x(:,:,:,e)
-     end do
-  end do
-  do e = 1, nelv
-     do i = 1,lx
-        do j = 1, lx
-           if (hom_dir_el(e) .eq. 1) then
-              u%x(:,i,j,e) = avg_u%x(1,i,j,e)
-           else if (hom_dir_el(e) .eq. 2) then
-              u%x(i,:,j,e) = avg_u%x(i,1,j,e)
-           else if (hom_dir_el(e) .eq. 3) then
-              u%x(i,j,:,e) = avg_u%x(i,j,1,e)
-           end if
-        end do
-     end do
-  end do
-end subroutine perform_global_summation
-
-subroutine perform_local_summation(u_out, u, el_heights,domain_height, hom_dir_el, coef, nelv,lx)
-  use neko
-  implicit none
-  type(field_t), intent(inout) :: u, u_out, el_heights
-  type(coef_t), intent(inout) :: coef
-  integer, intent(in) :: nelv, lx
-  integer, intent(in) :: hom_dir_el(nelv)
-  real(kind=rp), intent(in) :: domain_height
-  real(kind=rp) :: wt
-  integer :: n, i, j, e, k
-
-  n = nelv*lx**3
-
-  call col2(u%x,el_heights%x,n)
-  call cmult(u%x, 1.0_rp/(2.0*domain_height),n)
-  call rzero(u_out%x,n)
-
-
-  do e = 1, nelv
-     do i = 1,lx
-        do j = 1, lx
-           do k = 1, lx
-              wt = coef%Xh%wx(k)
-              if (hom_dir_el(e) .eq. 1) then
-                 u_out%x(1,i,j,e) = u_out%x(1,i,j,e)+wt*u%x(k,i,j,e)
-              else if (hom_dir_el(e) .eq. 2) then
-                 u_out%x(i,1,j,e) = u_out%x(i,1,j,e)+wt*u%x(i,k,j,e)
-              else if (hom_dir_el(e) .eq. 3) then
-                 u_out%x(i,j,1,e) = u_out%x(i,j,1,e)+wt*u%x(i,j,k,e)
-              end if
-           end do
-        end do
-     end do
-  end do
-
-  do e = 1, nelv
-     do i = 1,lx
-        do j = 1, lx
-           if (hom_dir_el(e) .eq. 1) then
-              u_out%x(:,i,j,e) = u_out%x(1,i,j,e)
-           else if (hom_dir_el(e) .eq. 2) then
-              u_out%x(i,:,j,e) = u_out%x(i,1,j,e)
-           else if (hom_dir_el(e) .eq. 3) then
-              u_out%x(i,j,:,e) = u_out%x(i,j,1,e)
-           end if
-        end do
-     end do
-  end do
-end subroutine perform_local_summation

contrib/calc_lift_from_field/calc_lift_from_field.f90

@@ -116,7 +116,7 @@ program calc_lift_from_field
   else 
      call neko_error('The homogeneous direction should be "x", "y"or "z"')
   end if
-  call map_1d%init(dof, gs_h, dir, 1e-7_rp)
+  call map_1d%init(coef, dir, 1e-7_rp)
 
 
   call u%init(dof)

src/.depends

@@ -22,7 +22,8 @@ math/math.o : math/math.f90 comm/comm.o config/num_types.o
 math/mathops.o : math/mathops.f90 config/num_types.o 
 math/fast3d.o : math/fast3d.f90 math/math.o sem/speclib.o config/num_types.o 
 sem/space.o : sem/space.f90 math/mxm_wrapper.o math/tensor.o math/math.o math/fast3d.o common/utils.o device/device.o sem/speclib.o config/num_types.o config/neko_config.o 
-sem/map_1d.o : sem/map_1d.f90 comm/mpi_types.o math/math.o common/utils.o common/log.o comm/comm.o device/device.o mesh/mesh.o gs/gather_scatter.o sem/dofmap.o sem/space.o config/num_types.o 
+sem/map_1d.o : sem/map_1d.f90 comm/mpi_types.o math/math.o common/utils.o common/log.o math/vector.o math/matrix.o field/field_list.o sem/coef.o comm/comm.o device/device.o mesh/mesh.o gs/gather_scatter.o sem/dofmap.o sem/space.o config/num_types.o 
+sem/map_2d.o : sem/map_2d.f90 field/field.o io/fld_file_data.o comm/mpi_types.o math/math.o common/log.o math/vector.o math/matrix.o sem/coef.o field/field_list.o comm/comm.o common/utils.o device/device.o mesh/mesh.o gs/gather_scatter.o sem/map_1d.o sem/dofmap.o sem/space.o config/num_types.o 
 sem/dofmap.o : sem/dofmap.f90 mesh/hex.o mesh/quad.o mesh/element.o math/math.o device/device.o math/tensor.o math/fast3d.o common/utils.o config/num_types.o adt/tuple.o sem/space.o mesh/mesh.o config/neko_config.o 
 sem/coef.o : sem/coef.f90 device/device.o math/mxm_wrapper.o sem/bcknd/device/device_coef.o math/bcknd/device/device_math.o mesh/mesh.o math/math.o sem/space.o sem/dofmap.o config/num_types.o config/neko_config.o gs/gs_ops.o gs/gather_scatter.o 
 sem/cpr.o : sem/cpr.f90 sem/dofmap.o common/log.o math/mxm_wrapper.o math/tensor.o sem/coef.o mesh/mesh.o math/math.o sem/space.o field/field.o config/num_types.o config/neko_config.o gs/gather_scatter.o 
@@ -304,5 +305,5 @@ wall_models/wall_model.o : wall_models/wall_model.f90 common/log.o comm/comm.o m
 wall_models/rough_log_law.o : wall_models/rough_log_law.f90 common/utils.o common/json_utils.o field/field_registry.o wall_models/wall_model.o config/neko_config.o sem/coef.o sem/dofmap.o config/num_types.o field/field.o 
 wall_models/spalding.o : wall_models/spalding.f90 common/utils.o common/log.o common/json_utils.o field/field_registry.o wall_models/wall_model.o config/neko_config.o sem/coef.o sem/dofmap.o config/num_types.o field/field.o 
 wall_models/wall_model_fctry.o : wall_models/wall_model_fctry.f90 common/json_utils.o common/utils.o wall_models/rough_log_law.o wall_models/spalding.o les/vreman.o wall_models/wall_model.o 
-neko.o : neko.f90 common/json_utils.o common/runtime_statistics.o bc/field_dirichlet_vector.o bc/field_dirichlet.o mesh/point_zone_registry.o mesh/point_zones/sphere_point_zone.o mesh/point_zones/box_point_zone.o mesh/point_zone.o sem/point_interpolator.o common/time_interpolator.o io/data_streamer.o simulation_components/simcomp_executor.o field/scratch_registry.o field/field_registry.o qoi/drag_torque.o common/system.o sem/spectral_error_indicator.o simulation_components/probes.o simulation_components/simulation_component.o math/tensor.o math/matrix.o math/vector.o scalar/scalar_user_source_term.o fluid/fluid_user_source_term.o field/field_list.o fluid/fluid_stats.o sem/cpr.o sem/map_1d.o math/bcknd/device/device_math.o device/device.o common/jobctrl.o common/signal.o common/user_intf.o common/projection.o math/mathops.o math/operators.o simulation.o io/output.o common/sampler.o case.o config/neko_config.o comm/parmetis.o math/ax.o bc/dirichlet.o bc/wall.o bc/bc.o sem/coef.o gs/gather_scatter.o comm/mpi_types.o field/field.o io/file.o common/global_interpolation.o math/mxm_wrapper.o io/format/map.o field/mesh_field.o mesh/point.o mesh/mesh.o adt/tuple.o adt/stack.o adt/uset.o adt/htable.o sem/space.o sem/dofmap.o sem/speclib.o math/math.o common/log.o common/utils.o comm/comm.o config/num_types.o 
+neko.o : neko.f90 common/json_utils.o common/runtime_statistics.o bc/field_dirichlet_vector.o bc/field_dirichlet.o mesh/point_zone_registry.o mesh/point_zones/sphere_point_zone.o mesh/point_zones/box_point_zone.o mesh/point_zone.o sem/point_interpolator.o common/time_interpolator.o io/data_streamer.o simulation_components/simcomp_executor.o field/scratch_registry.o field/field_registry.o qoi/drag_torque.o common/system.o sem/spectral_error_indicator.o simulation_components/probes.o simulation_components/simulation_component.o math/tensor.o math/matrix.o math/vector.o scalar/scalar_user_source_term.o fluid/fluid_user_source_term.o field/field_list.o fluid/fluid_stats.o sem/cpr.o sem/map_2d.o sem/map_1d.o math/bcknd/device/device_math.o device/device.o common/jobctrl.o common/signal.o common/user_intf.o common/projection.o math/mathops.o math/operators.o simulation.o io/output.o common/sampler.o case.o config/neko_config.o comm/parmetis.o math/ax.o bc/dirichlet.o bc/wall.o bc/bc.o sem/coef.o gs/gather_scatter.o comm/mpi_types.o field/field.o io/file.o common/global_interpolation.o math/mxm_wrapper.o io/format/map.o field/mesh_field.o mesh/point.o mesh/mesh.o adt/tuple.o adt/stack.o adt/uset.o adt/htable.o sem/space.o sem/dofmap.o sem/speclib.o math/math.o common/log.o common/utils.o comm/comm.o config/num_types.o 
 driver.o : driver.f90 neko.o 

src/Makefile.am

@@ -26,6 +26,7 @@ neko_fortran_SOURCES = \
 	math/fast3d.f90\
 	sem/space.f90\
 	sem/map_1d.f90\
+	sem/map_2d.f90\
 	sem/dofmap.f90\
 	sem/coef.f90\
 	sem/cpr.f90\

src/io/fld_file_data.f90

@@ -43,6 +43,8 @@ module fld_file_data
      procedure, pass(this) :: add => fld_file_data_add
      procedure, pass(this) :: size => fld_file_data_size
      procedure, pass(this) :: get_list => fld_file_get_list
+     procedure, pass(this) :: init_same => fld_file_init_same
+     procedure, pass(this) :: init_n_fields => fld_file_init_n_fields
   end type fld_file_data_t
 
 contains
@@ -70,6 +72,71 @@ function fld_file_data_size(this) result(i)
 
   end function fld_file_data_size
 
+  !> Genereate same fields as in another fld_file
+  subroutine fld_file_init_same(this, fld_file, n)
+    class(fld_file_data_t), target, intent(inout) :: this
+    class(fld_file_data_t), target, intent(in) :: fld_file
+    integer, intent(in) :: n
+    integer :: i, j
+
+    if(fld_file%u%n .gt. 0) then
+       call this%u%init(n)
+    end if
+    if(fld_file%v%n .gt. 0) then
+       call this%v%init(n)
+    end if
+    if(fld_file%w%n .gt. 0) then
+       call this%w%init(n)
+    end if
+    if(fld_file%p%n .gt. 0) then
+       call this%p%init(n)
+    end if
+    if(fld_file%t%n .gt. 0) then
+       call this%t%init(n)
+    end if
+    this%n_scalars = fld_file%n_scalars
+    allocate(this%s(fld_file%n_scalars))
+    do j = 1, fld_file%n_scalars
+       call this%s(j)%init(n)
+    end do
+
+  end subroutine fld_file_init_same
+
+  !> Genereate same fields as in another fld_file
+  subroutine fld_file_init_n_fields(this, n_fields, n)
+    class(fld_file_data_t), target, intent(inout) :: this
+    integer, intent(in) :: n, n_fields
+    integer :: i, j
+
+
+    if(n_fields .gt. 0) then
+       call this%u%init(n)
+    end if
+    if(n_fields .gt. 1) then
+       call this%v%init(n)
+    end if
+    if(n_fields .gt. 2) then
+       call this%w%init(n)
+    end if
+    if(n_fields .gt. 3) then
+       call this%p%init(n)
+    end if
+    if(n_fields .gt. 4) then
+       call this%t%init(n)
+    end if 
+    if (n_fields .gt. 5) then
+       this%n_scalars = n_fields-5
+       allocate(this%s(this%n_scalars))
+       do j = 1, this%n_scalars
+          call this%s(j)%init(n)
+       end do
+    end if
+
+  end subroutine fld_file_init_n_fields
+
+
+
+
   !> Get a list with pointers to the fields in the fld file
   subroutine fld_file_get_list(this, ptr_list, n)
     class(fld_file_data_t), target, intent(in) :: this
@@ -157,6 +224,7 @@ subroutine fld_file_data_free(this)
        do i = 1, this%n_scalars
           call this%s(i)%free()
        end do
+       deallocate(this%s)
     end if
     this%n_scalars = 0
     this%time = 0.0
@@ -169,6 +237,7 @@ subroutine fld_file_data_free(this)
     this%t_counter = 0
     this%meta_nsamples = 0
     this%meta_start_counter = 0
+    if(allocated(this%idx)) deallocate(this%idx)
   end subroutine fld_file_data_free
 
 end module fld_file_data

src/neko.f90

@@ -89,6 +89,7 @@ module neko
        device_glsc3_many, device_add2s2_many, device_glsc2, device_glsum, &
        device_masked_copy, device_cfill_mask, device_add3, device_cadd2
   use map_1d, only : map_1d_t
+  use map_2d, only : map_2d_t
   use cpr, only : cpr_t, cpr_init, cpr_free
   use fluid_stats, only : fluid_stats_t
   use field_list, only : field_list_t