From e4fe8c9222f112eec81ac862bc0132f66c5d1b3b Mon Sep 17 00:00:00 2001
From: Martin Karp <makarp@n173-p15.eduroam.kth.se>
Date: Thu, 19 Sep 2024 17:18:55 +0200
Subject: [PATCH 1/9] add things in map_1d

---
 src/sem/map_1d.f90 | 315 ++++++++++++++++++++++++++++++++-------------
 1 file changed, 226 insertions(+), 89 deletions(-)

diff --git a/src/sem/map_1d.f90 b/src/sem/map_1d.f90
index 28481d6f23c..cfb6056818c 100644
--- a/src/sem/map_1d.f90
+++ b/src/sem/map_1d.f90
@@ -8,14 +8,18 @@ module map_1d
   use mesh, only: mesh_t
   use device
   use comm
+  use coefs, only: coef_t
+  use field_list, only: field_list_t
+  use matrix, only: matrix_t
+  use vector, only: vector_ptr_t
   use logger, only: neko_log, LOG_SIZE
   use utils, only: neko_error, neko_warning
-  use math, only: glmax, glmin, glimax, glimin, relcmp
+  use math, only: glmax, glmin, glimax, glimin, relcmp, cmult, add2s1, col2
   use neko_mpi_types
   use, intrinsic :: iso_c_binding
   implicit none
   private
-  !> Type that encapsulates a mapping from each gll point in the mesh 
+  !> Type that encapsulates a mapping from each gll point in the mesh
   !! to its corresponding (global) GLL point index in one direction.
   !! @remark Could also be rather easily extended to say polar coordinates as well.
   type, public :: map_1d_t
@@ -24,70 +28,90 @@ module map_1d
      !> Checks which level an element belongs to.
      integer, allocatable :: el_lvl(:)
      !> Checks which level or id in the 1D GLL mapping each point in the dofmap is.
-     integer, allocatable :: pt_lvl(:,:,:,:) 
+     integer, allocatable :: pt_lvl(:,:,:,:)
      !> Number of elements stacked on top of eachother in the specified direction
      integer :: n_el_lvls
+     !> Number of total gll levels
+     integer :: n_gll_lvls
+     !> Dofmap
      !> Dofmap
      type(dofmap_t), pointer :: dof => null()
+     !> SEM coefs
+     type(coef_t), pointer :: coef => null()
      !> Mesh
      type(mesh_t), pointer :: msh => null()
      !> The specified direction in which we create the 1D mapping
      integer :: dir
      !> Tolerance for the mesh
      real(kind=rp) :: tol = 1e-7
+     !> Volume per level in the 1d grid
+     real(kind=rp), allocatable :: volume_per_gll_lvl(:)
    contains
      !> Constructor
-     procedure, pass(this) :: init => map_1d_init
+     procedure, pass(this) :: init_int => map_1d_init
+
+     procedure, pass(this) :: init_char => map_1d_init_char
+     generic :: init => init_int, init_char
      !> Destructor
      procedure, pass(this) :: free => map_1d_free
+     !> Average field list along planes
+     procedure, pass(this) :: average_planes_fld_lst =>  map_1d_average_field_list
+
+     procedure, pass(this) :: average_planes_vec_ptr =>  map_1d_average_vector_ptr
+     generic :: average_planes => average_planes_fld_lst, average_planes_vec_ptr
   end type map_1d_t
 
 
 contains
-  
-  subroutine map_1d_init(this, dof, gs, dir, tol)
+
+  subroutine map_1d_init(this, coef,  dir, tol)
     class(map_1d_t) :: this
-    type(dofmap_t), target, intent(in) :: dof
-    type(gs_t), intent(inout) :: gs
+    type(coef_t), intent(inout), target :: coef
     integer, intent(in) :: dir
     real(kind=rp), intent(in) :: tol
-    integer :: nelv, lx, n, i, e, lvl
+    integer :: nelv, lx, n, i, e, lvl, ierr
     real(kind=rp), contiguous, pointer :: line(:,:,:,:)
     real(kind=rp), allocatable :: min_vals(:,:,:,:)
+    real(kind=rp), allocatable :: min_temp(:,:,:,:)
     type(c_ptr) :: min_vals_d = c_null_ptr
     real(kind=rp) :: el_dim(3,3), glb_min, glb_max, el_min
+
     call this%free()
 
     if (NEKO_BCKND_DEVICE .eq. 1) then
-        if (pe_rank .eq. 0) then
-           call neko_warning('map_1d does not copy indices to device, but ok if used on cpu')
-        end if
+       if (pe_rank .eq. 0) then
+          call neko_warning('map_1d does not copy indices to device, but ok if used on cpu and for io')
+       end if
     end if
-    
+
     this%dir = dir
-    this%dof => dof
-    this%msh => dof%msh
+    this%dof => coef%dof
+    this%coef => coef
+    this%msh => coef%msh
     nelv = this%msh%nelv
     lx = this%dof%Xh%lx
-    n = dof%size()
-  
+    n = this%dof%size()
+
     if (dir .eq. 1) then
-        line => dof%x
+       line => this%dof%x
     else if (dir .eq. 2) then
-        line => dof%y
+       line => this%dof%y
     else if(dir .eq. 3) then
-        line => dof%z
+       line => this%dof%z
     else
-        call neko_error('Invalid dir for geopmetric comm')
+       call neko_error('Invalid dir for geopmetric comm')
     end if
-    
     allocate(this%dir_el(nelv))
     allocate(this%el_lvl(nelv))
-    allocate(min_vals(lx, lx, lx, nelv))
     allocate(this%pt_lvl(lx, lx, lx, nelv))
+    allocate(min_vals(lx, lx, lx, nelv))
+    allocate(min_temp(lx, lx, lx, nelv))
+    call MPI_BARRIER(NEKO_COMM)
     if (NEKO_BCKND_DEVICE .eq. 1) then
+
        call device_map(min_vals,min_vals_d,n)
     end if
+    call MPI_BARRIER(NEKO_COMM)
 
     do i = 1, nelv
        !store which direction r,s,t corresponds to speciifed direction, x,y,z
@@ -113,7 +137,7 @@ subroutine map_1d_init(this, dof, gs, dir, tol)
     do e = 1, nelv
        el_min = minval(line(:,:,:,e))
        min_vals(:,:,:,e) = el_min
-       ! Check if this element is on the bottom, in this case assign el_lvl = i = 1 
+       ! Check if this element is on the bottom, in this case assign el_lvl = i = 1
        if (relcmp(el_min, glb_min, this%tol)) then
           if(this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
        end if
@@ -122,88 +146,129 @@ subroutine map_1d_init(this, dof, gs, dir, tol)
     ! When the minumum value has propagated to the highest level this stops.
     ! Only works when the bottom plate of the domain is flat.
     do while (.not. relcmp(glmax(min_vals,n), glb_min, this%tol))
-      i = i + 1
-      do e = 1, nelv
-         !Sets the value at the bottom of each element to glb_max
-         if (this%dir_el(e) .eq. 1) then
-            if (line(1,1,1,e) .gt. line(lx,1,1,e)) then
-               min_vals(lx,:,:,e) = glb_max
-            else
-               min_vals(1,:,:,e) = glb_max
-            end if
-         end if
-         if (this%dir_el(e) .eq. 2) then
-            if (line(1,1,1,e) .gt. line(1,lx,1,e)) then
-               min_vals(:,lx,:,e) = glb_max
-            else
-               min_vals(:,1,:,e) = glb_max
-            end if
-         end if
-         if (this%dir_el(e) .eq. 3) then
-            if (line(1,1,1,e) .gt. line(1,1,lx,e)) then
-               min_vals(:,:,lx,e) = glb_max
-            else
-               min_vals(:,:,1,e) = glb_max
-            end if
-         end if
-      end do
-      if (NEKO_BCKND_DEVICE .eq. 1) &
+       i = i + 1
+       do e = 1, nelv
+          !Sets the value at the bottom of each element to glb_max
+          if (this%dir_el(e) .eq. 1) then
+             if (line(1,1,1,e) .gt. line(lx,1,1,e)) then
+                min_vals(lx,:,:,e) = glb_max
+             else
+                min_vals(1,:,:,e) = glb_max
+             end if
+          end if
+          if (this%dir_el(e) .eq. 2) then
+             if (line(1,1,1,e) .gt. line(1,lx,1,e)) then
+                min_vals(:,lx,:,e) = glb_max
+             else
+                min_vals(:,1,:,e) = glb_max
+             end if
+          end if
+          if (this%dir_el(e) .eq. 3) then
+             if (line(1,1,1,e) .gt. line(1,1,lx,e)) then
+                min_vals(:,:,lx,e) = glb_max
+             else
+                min_vals(:,:,1,e) = glb_max
+             end if
+          end if
+       end do
+       !Make sketchy min as GS_OP_MIN is not supported with device mpi
+       min_temp = min_vals
+       if (NEKO_BCKND_DEVICE .eq. 1) &
          call device_memcpy(min_vals, min_vals_d, n,&
                             HOST_TO_DEVICE, sync=.false.)
-      !Propagates the minumum value along the element boundary.
-      call gs%op(min_vals,n,GS_OP_MIN)
-      if (NEKO_BCKND_DEVICE .eq. 1) &
+       !Propagates the minumum value along the element boundary.
+       call coef%gs_h%op(min_vals, n, GS_OP_ADD)
+       if (NEKO_BCKND_DEVICE .eq. 1) &
           call device_memcpy(min_vals, min_vals_d, n,&
                              DEVICE_TO_HOST, sync=.true.)
-      !Checks the new minimum value on each element
-      !Assign this value to all points in this element in min_val
-      !If the element has not already been assinged a level, 
-      !and it has obtained the minval, set el_lvl = i
-      do e = 1, nelv
-         el_min = minval(min_vals(:,:,:,e))
-         min_vals(:,:,:,e) = el_min
-         if (relcmp(el_min, glb_min, this%tol)) then
-            if (this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
-         end if
-      end do
+       !Obtain average along boundary
+
+       call col2(min_vals, coef%mult, n)
+       call cmult(min_temp, -1.0_rp, n)
+       call add2s1(min_vals, min_temp, 2.0_rp, n)
+
+
+       !Checks the new minimum value on each element
+       !Assign this value to all points in this element in min_val
+       !If the element has not already been assinged a level,
+       !and it has obtained the minval, set el_lvl = i
+       do e = 1, nelv
+          el_min = minval(min_vals(:,:,:,e))
+          min_vals(:,:,:,e) = el_min
+          if (relcmp(el_min, glb_min, this%tol)) then
+             if (this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
+          end if
+       end do
     end do
     this%n_el_lvls = glimax(this%el_lvl,nelv)
+    this%n_gll_lvls = this%n_el_lvls*lx
     if ( pe_rank .eq. 0) then
        write(*,*) 'Number of element levels', this%n_el_lvls
+       write(*,*) 'Total number of levels', this%n_gll_lvls
     end if
     !Numbers the points in each element based on the element level
     !and its orientation
     do e = 1, nelv
        do i = 1, lx
-         lvl = lx * (this%el_lvl(e) - 1) + i
-         if (this%dir_el(e) .eq. 1) then
-            if (line(1,1,1,e) .gt. line(lx,1,1,e)) then
-               this%pt_lvl(lx-i+1,:,:,e) = lvl
-            else
-               this%pt_lvl(i,:,:,e) = lvl
-            end if
-         end if
-         if (this%dir_el(e) .eq. 2) then
-            if (line(1,1,1,e) .gt. line(1,lx,1,e)) then
-               this%pt_lvl(:,lx-i+1,:,e) = lvl
-            else
-               this%pt_lvl(:,i,:,e) = lvl
-            end if
-         end if
-         if (this%dir_el(e) .eq. 3) then
-            if (line(1,1,1,e) .gt. line(1,1,lx,e)) then
-               this%pt_lvl(:,:,lx-i+1,e) = lvl
-            else
-               this%pt_lvl(:,:,i,e) = lvl
-            end if
-         end if
+          lvl = lx * (this%el_lvl(e) - 1) + i
+          if (this%dir_el(e) .eq. 1) then
+             if (line(1,1,1,e) .gt. line(lx,1,1,e)) then
+                this%pt_lvl(lx-i+1,:,:,e) = lvl
+             else
+                this%pt_lvl(i,:,:,e) = lvl
+             end if
+          end if
+          if (this%dir_el(e) .eq. 2) then
+             if (line(1,1,1,e) .gt. line(1,lx,1,e)) then
+                this%pt_lvl(:,lx-i+1,:,e) = lvl
+             else
+                this%pt_lvl(:,i,:,e) = lvl
+             end if
+          end if
+          if (this%dir_el(e) .eq. 3) then
+             if (line(1,1,1,e) .gt. line(1,1,lx,e)) then
+                this%pt_lvl(:,:,lx-i+1,e) = lvl
+             else
+                this%pt_lvl(:,:,i,e) = lvl
+             end if
+          end if
        end do
     end do
-    call device_deassociate(min_vals)
-    call device_free(min_vals_d)
-    deallocate(min_vals)
+    if(allocated(min_vals)) deallocate(min_vals)
+    if(c_associated(min_vals_d)) call device_free(min_vals_d)
+    if(allocated(min_temp)) deallocate(min_temp)
+    allocate(this%volume_per_gll_lvl(this%n_gll_lvls))
+    this%volume_per_gll_lvl =0.0_rp
+    do i = 1, n
+       this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1)) = &
+       this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1)) + coef%B(i,1,1,1)
+    end do
+    call MPI_Allreduce(MPI_IN_PLACE,this%volume_per_gll_lvl, this%n_gll_lvls, &
+         MPI_REAL_PRECISION, MPI_SUM, NEKO_COMM, ierr)
+
   end subroutine map_1d_init
 
+  subroutine map_1d_init_char(this, coef,  dir, tol)
+    class(map_1d_t) :: this
+    type(coef_t), intent(inout), target :: coef
+    character(len=*), intent(in) :: dir
+    real(kind=rp), intent(in) :: tol
+    integer :: idir
+
+    if (trim(dir) .eq. 'yz' .or. trim(dir) .eq. 'zy') then
+       idir = 1
+    else if (trim(dir) .eq. 'xz' .or. trim(dir) .eq. 'zx') then
+       idir = 2
+    else if (trim(dir) .eq. 'xy' .or. trim(dir) .eq. 'yx') then
+       idir = 3
+    else
+       call neko_error('homogenous direction not supported')
+    end if
+
+    call  this%init(coef,idir,tol)
+
+  end subroutine map_1d_init_char
+
   subroutine map_1d_free(this)
     class(map_1d_t) :: this
 
@@ -212,9 +277,81 @@ subroutine map_1d_free(this)
     if(allocated(this%pt_lvl)) deallocate(this%pt_lvl)
     if(associated(this%dof)) nullify(this%dof)
     if(associated(this%msh)) nullify(this%msh)
+    if(associated(this%coef)) nullify(this%coef)
+    if(allocated(this%volume_per_gll_lvl)) deallocate(this%volume_per_gll_lvl)
     this%dir = 0
     this%n_el_lvls = 0
+    this%n_gll_lvls = 0
 
   end subroutine map_1d_free
 
+  subroutine map_1d_average_field_list(this, avg_planes, field_list)
+    class(map_1d_t), intent(inout) :: this
+    type(field_list_t), intent(inout) :: field_list
+    type(matrix_t), intent(inout) :: avg_planes
+    integer :: n, ierr, j, i
+    real(kind=rp) :: coord
+    call avg_planes%free()
+    call avg_planes%init(this%n_gll_lvls, field_list%size()+1)
+    avg_planes = 0.0_rp
+    !ugly way of getting coordinates, computes average
+    n = this%dof%size()
+    do i = 1, n
+       if (this%dir .eq. 1) coord = this%dof%x(i,1,1,1)
+       if (this%dir .eq. 2) coord = this%dof%y(i,1,1,1)
+       if (this%dir .eq. 3) coord = this%dof%z(i,1,1,1)
+       avg_planes%x(this%pt_lvl(i,1,1,1),1) = &
+       avg_planes%x(this%pt_lvl(i,1,1,1),1) + coord*this%coef%B(i,1,1,1) &
+       /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
+    end do
+    do j = 2, field_list%size() + 1
+       do i = 1, n
+          avg_planes%x(this%pt_lvl(i,1,1,1),j) = &
+          avg_planes%x(this%pt_lvl(i,1,1,1),j) + field_list%items(j-1)%ptr%x(i,1,1,1)*this%coef%B(i,1,1,1) &
+          /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
+       end do
+    end do
+    if (pe_size .gt. 1) then
+       call MPI_Allreduce(MPI_IN_PLACE, avg_planes%x, (field_list%size()+1)*this%n_gll_lvls, &
+            MPI_REAL_PRECISION, MPI_SUM, NEKO_COMM, ierr)
+    end if
+
+
+  end subroutine map_1d_average_field_list
+
+  subroutine map_1d_average_vector_ptr(this, avg_planes, vector_ptr)
+    class(map_1d_t), intent(inout) :: this
+    !Observe is an array...
+    type(vector_ptr_t), intent(inout) :: vector_ptr(:)
+    type(matrix_t), intent(inout) :: avg_planes
+    integer :: n, ierr, j, i
+    real(kind=rp) :: coord
+
+    call avg_planes%free()
+    call avg_planes%init(this%n_gll_lvls,size(vector_ptr)+1)
+    !ugly way of getting coordinates, computes average
+    avg_planes = 0.0_rp
+
+    n = this%dof%size()
+    do i = 1, n
+       if (this%dir .eq. 1) coord = this%dof%x(i,1,1,1)
+       if (this%dir .eq. 2) coord = this%dof%y(i,1,1,1)
+       if (this%dir .eq. 3) coord = this%dof%z(i,1,1,1)
+       avg_planes%x(this%pt_lvl(i,1,1,1),1) = &
+       avg_planes%x(this%pt_lvl(i,1,1,1),1) + coord*this%coef%B(i,1,1,1) &
+       /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
+    end do
+    do j = 2, size(vector_ptr)+1
+       do i = 1, n
+          avg_planes%x(this%pt_lvl(i,1,1,1),j) = &
+          avg_planes%x(this%pt_lvl(i,1,1,1),j) + vector_ptr(j-1)%ptr%x(i)*this%coef%B(i,1,1,1) &
+          /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
+       end do
+    end do
+    call MPI_Allreduce(MPI_IN_PLACE,avg_planes%x, (size(vector_ptr)+1)*this%n_gll_lvls, &
+       MPI_REAL_PRECISION, MPI_SUM, NEKO_COMM, ierr)
+
+
+  end subroutine map_1d_average_vector_ptr
+
 end module map_1d

From 6381244713923ac72ba394593035c0f9582ecbd6 Mon Sep 17 00:00:00 2001
From: Martin Karp <makarp@n173-p15.eduroam.kth.se>
Date: Thu, 19 Sep 2024 17:55:53 +0200
Subject: [PATCH 2/9] add map_2d

---
 .../average_field_in_space.f90                | 212 +--------
 .../calc_lift_from_field.f90                  |   4 +-
 src/.depends                                  |   3 +-
 src/Makefile.am                               |   1 +
 src/neko.f90                                  |   1 +
 src/sem/map_1d.f90                            |  13 +-
 src/sem/map_2d.f90                            | 415 ++++++++++++++++++
 7 files changed, 457 insertions(+), 192 deletions(-)
 create mode 100644 src/sem/map_2d.f90

diff --git a/contrib/average_field_in_space/average_field_in_space.f90 b/contrib/average_field_in_space/average_field_in_space.f90
index 73174cdf22b..e783fadede5 100644
--- a/contrib/average_field_in_space/average_field_in_space.f90
+++ b/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,33 @@ 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
+     print *,'whatup in loop', field_data%meta_nsamples
+    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)
+        print *,'lol'
+        call output_file%write(output_data,field_data%time)
+        print *,'lol2'
      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
diff --git a/contrib/calc_lift_from_field/calc_lift_from_field.f90 b/contrib/calc_lift_from_field/calc_lift_from_field.f90
index 0718b9fe7f7..b62beefeb86 100644
--- a/contrib/calc_lift_from_field/calc_lift_from_field.f90
+++ b/contrib/calc_lift_from_field/calc_lift_from_field.f90
@@ -1,4 +1,4 @@
-!> Program to calculate the force and acting on a single boundary zone as well as
+	!> Program to calculate the force and acting on a single boundary zone as well as
 !! the torque around a point if one changes the value of center (defaults to 0,0,0).
 !! Outputs the x,y,z prjections of the pressure and viscous forces and torques and
 !! additionally saves the distribution of these quantities along a selected homogenous
@@ -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)
diff --git a/src/.depends b/src/.depends
index f84fb56e866..fbf5f876676 100644
--- a/src/.depends
+++ b/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 
diff --git a/src/Makefile.am b/src/Makefile.am
index bd592944573..9c639b61ff5 100644
--- a/src/Makefile.am
+++ b/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\
diff --git a/src/neko.f90 b/src/neko.f90
index a9c27e76303..68b04f9f760 100644
--- a/src/neko.f90
+++ b/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
diff --git a/src/sem/map_1d.f90 b/src/sem/map_1d.f90
index cfb6056818c..278cdd31fb6 100644
--- a/src/sem/map_1d.f90
+++ b/src/sem/map_1d.f90
@@ -285,6 +285,12 @@ subroutine map_1d_free(this)
 
   end subroutine map_1d_free
 
+
+  !> Computes average if field list in two directions and outputs matrix
+  !! with averaged values
+  !! avg_planes contains coordinates in first row, avg. of fields in the rest
+  !! @param avg_planes output averages
+  !! @param field_list list of fields to be averaged
   subroutine map_1d_average_field_list(this, avg_planes, field_list)
     class(map_1d_t), intent(inout) :: this
     type(field_list_t), intent(inout) :: field_list
@@ -318,7 +324,12 @@ subroutine map_1d_average_field_list(this, avg_planes, field_list)
 
 
   end subroutine map_1d_average_field_list
-
+  
+  !> Computes average if vector_pt in two directions and outputs matrix
+  !! with averaged values
+  !! avg_planes contains coordinates in first row, avg. of fields in the rest
+  !! @param avg_planes output averages
+  !! @param vector_pts to vectors to be averaged
   subroutine map_1d_average_vector_ptr(this, avg_planes, vector_ptr)
     class(map_1d_t), intent(inout) :: this
     !Observe is an array...
diff --git a/src/sem/map_2d.f90 b/src/sem/map_2d.f90
new file mode 100644
index 00000000000..bb18fc7fb4b
--- /dev/null
+++ b/src/sem/map_2d.f90
@@ -0,0 +1,415 @@
+!> Maps a 3D dofmap to a 2D spectral element grid.
+
+module map_2d
+  use num_types, only: rp
+  use space, only: space_t
+  use dofmap, only: dofmap_t
+  use map_1d
+  use gather_scatter
+  use mesh, only: mesh_t
+  use device
+  use utils
+  use comm
+  use field_list
+  use coefs, only: coef_t
+  use field_list, only: field_list_t
+  use matrix, only: matrix_t
+  use vector, only: vector_ptr_t
+  use logger, only: neko_log, LOG_SIZE
+  use utils, only: neko_error, neko_warning
+  use math, only: glmax, glmin, glimax, glimin, relcmp, cmult, add2s1, col2, copy, rzero
+  use neko_mpi_types
+  use fld_file_data
+  use field
+  use, intrinsic :: iso_c_binding
+  implicit none
+  private
+  
+  type, public :: map_2d_t
+     integer :: nelv_2d = 0 !< Number of elements in 2D mesh on this rank
+     integer :: glb_nelv_2d = 0 !< global number of elements in 2d
+     integer :: offset_el_2d = 0 !< element offset for this rank
+     integer :: lxy = 0 !< number of gll points per 2D element
+     integer :: n_2d = 0 !< total number of gll points (nelv_2d*lxy)
+     integer, allocatable :: idx_2d(:) !< Mapping of GLL point from 3D to 2D
+     integer, allocatable :: el_idx_2d(:) !< Mapping of element in 3D to 2D
+     type(map_1d_t) :: map_1d !< 1D map in normal direction to 2D plane
+     type(mesh_t), pointer :: msh !< 3D mesh
+     type(dofmap_t), pointer :: dof => null() !< 3D dofmap
+     type(coef_t), pointer :: coef => null() !< 3D SEM coefs
+     type(field_t) :: u !< Work array 1, naming based on when they are used
+     type(field_t) :: old_u !< Work array 2
+     type(field_t) :: avg_u !< Work array 3
+     type(field_t) :: el_heights !< Weight elements by their size in integral
+     integer :: dir !< direction normal to 2D plane
+     real(kind=rp) :: domain_height !< total height of 3D domain
+   contains
+     procedure, pass(this) :: init_int => map_2d_init
+     procedure, pass(this) :: init_char => map_2d_init_char
+     generic :: init => init_int, init_char
+     procedure, pass(this) :: average_file => map_2d_average
+     procedure, pass(this) :: average_list => map_2d_average_field_list
+     generic :: average => average_list, average_file
+  end type map_2d_t
+
+contains
+
+  subroutine map_2d_init(this,coef, dir, tol)
+    class(map_2d_t), intent(inout) :: this
+    type(coef_t), intent(inout), target :: coef
+    integer, intent(in) :: dir
+    real(kind=rp), intent(in) :: tol
+    real(kind=rp) :: el_dim(3,3), el_h
+    integer :: i, e, j, ierr, k, lx, lxy, n
+    call this%map_1d%init(coef,dir,tol)
+    this%msh => coef%msh
+    this%coef => coef
+    this%dof => coef%dof
+    
+    call this%u%init(this%dof)
+    call this%old_u%init(this%dof)
+    call this%avg_u%init(this%dof)
+    call this%el_heights%init(this%dof)
+    this%dir = dir
+
+    n = this%dof%size()
+    lx = this%dof%Xh%lx
+    lxy = this%dof%Xh%lxy
+    this%nelv_2d = 0
+    do i = 1, this%msh%nelv
+       if (this%map_1d%el_lvl(i) .eq. 1) then
+          this%nelv_2d = this%nelv_2d + 1
+       end if
+    end do
+    this%glb_nelv_2d = 0
+    call MPI_Allreduce(this%nelv_2d, this%glb_nelv_2d, 1, &
+         MPI_INTEGER, MPI_SUM, NEKO_COMM, ierr)
+
+    this%offset_el_2d = 0
+    call MPI_Exscan(this%nelv_2d, this%offset_el_2d, 1, &
+         MPI_INTEGER, MPI_SUM, NEKO_COMM, ierr)
+    allocate(this%el_idx_2d(this%nelv_2d)) 
+    do i = 1, this%nelv_2d
+       this%el_idx_2d(i) = this%offset_el_2d + i
+    end do
+    this%n_2d = this%nelv_2d*lxy
+    allocate(this%idx_2d(this%n_2d))
+
+    j = 0
+    do e = 1, this%msh%nelv
+       if (this%map_1d%el_lvl(e) .eq. 1) then
+          if (this%map_1d%dir_el(e) .eq. 1) then
+             do i = 1, lxy
+                this%idx_2d(j*lxy+i) = linear_index(1,i,1,e,lx,lx,lx)
+             end do
+          end if
+          if (this%map_1d%dir_el(e) .eq. 2) then
+             do i = 1, lx
+                do k = 1, lx
+                   this%idx_2d(j*lxy+k+lx*(i-1)) = linear_index(k,1,i,e,lx,lx,lx)
+                end do
+             end do
+          end if
+          if (this%map_1d%dir_el(e) .eq. 3) then
+             do i = 1, lxy
+                this%idx_2d(j*lxy+i) = linear_index(i,1,1,e,lx,lx,lx)
+             end do
+          end if
+          j = j +1
+       end if
+    end do
+    do i = 1, this%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(this%msh%elements(i)%e%pts(1)%p%x-this%msh%elements(i)%e%pts(2)%p%x)
+       el_dim(2,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-this%msh%elements(i)%e%pts(3)%p%x)
+       el_dim(3,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-this%msh%elements(i)%e%pts(5)%p%x)
+       ! 1 corresponds to x, 2 to y, 3 to z
+       el_h = el_dim(this%map_1d%dir_el(i),dir)
+       this%el_heights%x(:,:,:,i) = el_h
+    end do
+    !Need to compute mapping for 3d to 2d
+    !Does order matter? Think its ok as long as values written in same order
+
+    call copy(this%u%x,this%el_heights%x,n)
+    call copy(this%old_u%x,this%el_heights%x,n)
+    call copy(this%avg_u%x,this%el_heights%x,n)
+    call perform_global_summation(this%u, this%avg_u, this%old_u, this%map_1d%n_el_lvls, &
+         this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, this%msh%nelv, lx)
+    this%domain_height = this%u%x(1,1,1,1)
+
+  end subroutine map_2d_init
+   
+  subroutine map_2d_init_char(this, coef,  dir, tol)
+    class(map_2d_t) :: this
+    type(coef_t), intent(inout), target :: coef
+    character(len=*), intent(in) :: dir
+    real(kind=rp), intent(in) :: tol
+    integer :: idir
+
+    if (trim(dir) .eq. 'x') then
+       idir = 1
+    else if (trim(dir) .eq. 'y') then
+       idir = 2
+    else if (trim(dir) .eq. 'z') then
+       idir = 3
+    else
+       call neko_error('Direction not supported, map_2d')
+    end if
+
+    call this%init(coef,idir,tol)
+
+  end subroutine map_2d_init_char
+
+
+  !> Computes average if field list in one direction and outputs 2D field
+  !! with averaged values. 
+  !! @param fld_data2D output 2D averages
+  !! @param field_list list of fields to be averaged
+  subroutine map_2d_average_field_list(this,fld_data2D,fld_data3D)
+    class(map_2d_t), intent(inout) :: this
+    type(fld_file_data_t), intent(inout) :: fld_data2D
+    type(field_list_t), intent(inout) :: fld_data3D
+    real(kind=rp), pointer, dimension(:,:,:,:) :: x_ptr, y_ptr
+    
+    type(vector_ptr_t), allocatable :: fields2d(:)
+    integer :: n_2d, n
+    integer :: i, j, lx, lxy, e
+
+    call fld_data2D%init(this%nelv_2d,this%offset_el_2d)
+    fld_data2D%gdim = 2
+    fld_data2D%lx = this%dof%Xh%lx
+    fld_data2D%ly = this%dof%Xh%ly
+    fld_data2D%lz = 1
+    lx = this%dof%Xh%lx
+    fld_data2D%glb_nelv = this%glb_nelv_2d
+    lxy = fld_data2D%lx*fld_data2D%ly
+    n_2d = lxy*this%nelv_2d
+    n = this%dof%size()
+    call fld_data2D%x%init(n_2d)
+    call fld_data2D%y%init(n_2d)
+    allocate(fld_data2D%idx(this%nelv_2d))
+
+    if (this%dir .eq. 1) then
+       x_ptr => this%dof%z
+       y_ptr => this%dof%y
+    end if
+    if (this%dir .eq. 2) then
+        x_ptr => this%dof%x
+        y_ptr => this%dof%z
+    end if
+    if (this%dir .eq. 3) then
+        x_ptr => this%dof%x
+        y_ptr => this%dof%y
+    end if
+    do j = 1, this%nelv_2d
+       fld_data2d%idx(j) = this%el_idx_2d(j)
+    end do
+    do j = 1, n_2d
+       fld_data2d%x%x(j) = x_ptr(this%idx_2d(j),1,1,1)
+       fld_data2d%y%x(j) = y_ptr(this%idx_2d(j),1,1,1)
+    end do
+    allocate(fields2D(fld_data3D%size()))
+
+    call fld_data2D%init_n_fields(fld_data3D%size(),n_2d)
+    this%u = 0.0_rp
+    this%old_u = 0.0_rp
+    this%avg_u = 0.0_rp
+    do i = 1, fld_data3D%size()
+       call copy(this%old_u%x,fld_data3D%items(i)%ptr%x,n)
+       call perform_local_summation(this%u,this%old_u, this%el_heights, this%domain_height, &
+            this%map_1d%dir_el, this%coef, this%msh%nelv, lx)
+       call copy(this%old_u%x,this%u%x,n)
+       call copy(this%avg_u%x,this%u%x,n)
+       call perform_global_summation(this%u, this%avg_u, this%old_u, this%map_1d%n_el_lvls, &
+            this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, this%msh%nelv, lx)
+       call copy(fld_data3D%items(i)%ptr%x,this%u%x,n)
+    end do
+    call fld_data2D%get_list(fields2d,fld_data2D%size())
+    do i = 1, fld_data3D%size()
+       do j = 1, n_2d
+          fields2d(i)%ptr%x(j) = fld_data3D%items(i)%ptr%x(this%idx_2d(j),1,1,1)
+       end do
+    end do
+
+  end subroutine map_2d_average_field_list
+ 
+  
+  !> Computes average if field list in one direction and outputs 2D field
+  !! with averaged values. 
+  !! @param fld_data2D output 2D averages
+  !! @param fld_data3D fld_file_data of fields to be averaged
+  subroutine map_2d_average(this,fld_data2D,fld_data3D)
+    class(map_2d_t), intent(inout) :: this
+    type(fld_file_data_t), intent(inout) :: fld_data2D
+    type(fld_file_data_t), intent(inout) :: fld_data3D
+    real(kind=rp), pointer, dimension(:,:,:,:) :: x_ptr, y_ptr
+    
+    type(vector_ptr_t), allocatable :: fields3d(:), fields2d(:)
+    integer :: n_2d, n
+    integer :: i, j, lx, lxy, e
+
+    call fld_data2D%init(this%nelv_2d,this%offset_el_2d)
+    fld_data2D%gdim = 2
+    fld_data2D%lx = this%dof%Xh%lx
+    fld_data2D%ly = this%dof%Xh%ly
+    fld_data2D%lz = 1
+    lx = this%dof%Xh%lx
+    fld_data2D%glb_nelv = this%glb_nelv_2d
+    lxy = fld_data2D%lx*fld_data2D%ly
+    n_2d = lxy*this%nelv_2d
+    n = this%dof%size()
+    call fld_data2D%x%init(n_2d)
+    call fld_data2D%y%init(n_2d)
+    allocate(fld_data2D%idx(n_2d))
+
+    if (this%dir .eq. 1) then
+       x_ptr => this%dof%z
+       y_ptr => this%dof%y
+    end if
+    if (this%dir .eq. 2) then
+        x_ptr => this%dof%x
+        y_ptr => this%dof%z
+    end if
+    if (this%dir .eq. 3) then
+        x_ptr => this%dof%x
+        y_ptr => this%dof%y
+    end if
+    do j = 1, n_2d
+       fld_data2d%idx(j) = this%idx_2d(j)
+       fld_data2d%x%x(j) = x_ptr(this%idx_2d(j),1,1,1)
+       fld_data2d%y%x(j) = y_ptr(this%idx_2d(j),1,1,1)
+    end do
+    allocate(fields3D(fld_data3D%size()))
+    allocate(fields2D(fld_data3D%size()))
+
+    call fld_data2D%init_n_fields(fld_data3D%size(),n_2d)
+    call fld_data3D%get_list(fields3D,fld_data3D%size())
+    
+    this%u = 0.0_rp
+    this%old_u = 0.0_rp
+    this%avg_u = 0.0_rp
+    do i = 1, fld_data3D%size()
+       call copy(this%old_u%x,fields3D(i)%ptr%x,n)
+       call perform_local_summation(this%u,this%old_u, this%el_heights, this%domain_height, &
+            this%map_1d%dir_el, this%coef, this%msh%nelv, lx)
+       call copy(this%old_u%x,this%u%x,n)
+       call copy(this%avg_u%x,this%u%x,n)
+       call perform_global_summation(this%u, this%avg_u, this%old_u, this%map_1d%n_el_lvls, &
+            this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, this%msh%nelv, lx)
+       call copy(fields3D(i)%ptr%x,this%u%x,n)
+    end do
+    call fld_data2D%get_list(fields2d,fld_data2D%size())
+    do i = 1, fld_data3D%size()
+       do j = 1, n_2d
+          fields2d(i)%ptr%x(j) = fields3D(i)%ptr%x(this%idx_2d(j))
+       end do
+    end do
+  end subroutine map_2d_average
+ 
+  subroutine perform_global_summation(u, avg_u, old_u, n_levels, hom_dir_el, gs_h, mult, nelv,lx)
+    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
+    type(c_ptr) :: ptr
+  
+    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=.true.)
+       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)
+       !Assumes sugarcube I think
+       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)
+    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
+
+end module map_2d

From 37978dba54ddcc93bc316e323c0a35827c937fb1 Mon Sep 17 00:00:00 2001
From: Martin Karp <makarp@n173-p15.eduroam.kth.se>
Date: Thu, 19 Sep 2024 18:02:46 +0200
Subject: [PATCH 3/9] deps

---
 src/.depends | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/.depends b/src/.depends
index fbf5f876676..7e40894fcf8 100644
--- a/src/.depends
+++ b/src/.depends
@@ -304,5 +304,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 

From d46be98b5bbf601c769e8db41157484a7b8ebeeb Mon Sep 17 00:00:00 2001
From: Martin Karp <makarp@n173-p15.eduroam.kth.se>
Date: Thu, 19 Sep 2024 18:14:15 +0200
Subject: [PATCH 4/9] fld data stuff

---
 src/io/fld_file_data.f90 | 70 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 70 insertions(+)

diff --git a/src/io/fld_file_data.f90 b/src/io/fld_file_data.f90
index 638402a31d4..a13472c4fd6 100644
--- a/src/io/fld_file_data.f90
+++ b/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,7 +224,9 @@ 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
     this%glb_nelv = 0
@@ -169,6 +238,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

From adf201d5457533aa630bf2f671580468a473cecf Mon Sep 17 00:00:00 2001
From: Martin Karp <33422901+MartinKarp@users.noreply.github.com>
Date: Fri, 20 Sep 2024 09:12:20 +0200
Subject: [PATCH 5/9] Update
 contrib/average_field_in_space/average_field_in_space.f90

Co-authored-by: Tim Felle Olsen <tife@dtu.dk>
---
 contrib/average_field_in_space/average_field_in_space.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/contrib/average_field_in_space/average_field_in_space.f90 b/contrib/average_field_in_space/average_field_in_space.f90
index e783fadede5..e01339ad6de 100644
--- a/contrib/average_field_in_space/average_field_in_space.f90
+++ b/contrib/average_field_in_space/average_field_in_space.f90
@@ -124,7 +124,7 @@ program average_field_in_space
      call neko_error('homogenous direction not supported')
   end if
   if (avg_to_1d) then
-     call map_1d%init(coef,  dir, 1e-7_rp)
+     call map_1d%init(coef, dir, 1e-7_rp)
   else
      call map_2d%init(coef,  dir, 1e-7_rp)
   end if

From 677bba6f6ea65abef0e6c8b168ab48c9266ce384 Mon Sep 17 00:00:00 2001
From: Martin Karp <33422901+MartinKarp@users.noreply.github.com>
Date: Fri, 20 Sep 2024 09:13:43 +0200
Subject: [PATCH 6/9] Apply suggestions from code review

Co-authored-by: Tim Felle Olsen <tife@dtu.dk>
---
 contrib/average_field_in_space/average_field_in_space.f90 | 2 +-
 contrib/calc_lift_from_field/calc_lift_from_field.f90     | 2 +-
 src/io/fld_file_data.f90                                  | 5 ++---
 3 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/contrib/average_field_in_space/average_field_in_space.f90 b/contrib/average_field_in_space/average_field_in_space.f90
index e01339ad6de..8fcd9865f68 100644
--- a/contrib/average_field_in_space/average_field_in_space.f90
+++ b/contrib/average_field_in_space/average_field_in_space.f90
@@ -126,7 +126,7 @@ program average_field_in_space
   if (avg_to_1d) then
      call map_1d%init(coef, dir, 1e-7_rp)
   else
-     call map_2d%init(coef,  dir, 1e-7_rp)
+     call map_2d%init(coef, dir, 1e-7_rp)
   end if
   
   !allocate array with pointers to all vectors in the file
diff --git a/contrib/calc_lift_from_field/calc_lift_from_field.f90 b/contrib/calc_lift_from_field/calc_lift_from_field.f90
index b62beefeb86..5628f9926a8 100644
--- a/contrib/calc_lift_from_field/calc_lift_from_field.f90
+++ b/contrib/calc_lift_from_field/calc_lift_from_field.f90
@@ -1,4 +1,4 @@
-	!> Program to calculate the force and acting on a single boundary zone as well as
+!> Program to calculate the force and acting on a single boundary zone as well as
 !! the torque around a point if one changes the value of center (defaults to 0,0,0).
 !! Outputs the x,y,z prjections of the pressure and viscous forces and torques and
 !! additionally saves the distribution of these quantities along a selected homogenous
diff --git a/src/io/fld_file_data.f90 b/src/io/fld_file_data.f90
index a13472c4fd6..4721de827a6 100644
--- a/src/io/fld_file_data.f90
+++ b/src/io/fld_file_data.f90
@@ -107,8 +107,8 @@ 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
@@ -226,7 +226,6 @@ subroutine fld_file_data_free(this)
        end do
        deallocate(this%s)
     end if
-   
     this%n_scalars = 0
     this%time = 0.0
     this%glb_nelv = 0

From 2f5feefbf9d8621ab8998b8abee6e0ea0bde1c66 Mon Sep 17 00:00:00 2001
From: Martin Karp <makarp@n173-p15.eduroam.kth.se>
Date: Fri, 20 Sep 2024 09:14:26 +0200
Subject: [PATCH 7/9] rm print

---
 contrib/average_field_in_space/average_field_in_space.f90 | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/contrib/average_field_in_space/average_field_in_space.f90 b/contrib/average_field_in_space/average_field_in_space.f90
index 8fcd9865f68..e1dd3be431a 100644
--- a/contrib/average_field_in_space/average_field_in_space.f90
+++ b/contrib/average_field_in_space/average_field_in_space.f90
@@ -133,7 +133,6 @@ program average_field_in_space
 
   output_file = file_t(trim(output_fname))
   do tstep = 0, field_data%meta_nsamples-1
-     print *,'whatup in loop', field_data%meta_nsamples
     if (pe_rank .eq. 0) write(*,*) 'Averaging field:', tstep
      if (tstep .gt. 0) call field_file%read(field_data)
      if (avg_to_1d) then
@@ -143,9 +142,7 @@ program average_field_in_space
         ! Should output a 2d field in principle
      else
         call map_2d%average(output_data,field_data)
-        print *,'lol'
         call output_file%write(output_data,field_data%time)
-        print *,'lol2'
      end if
   end do
   if (pe_rank .eq. 0) write(*,*) 'Done'

From ce93b250d6ae428e8bb012c7088aa1c02cffe396 Mon Sep 17 00:00:00 2001
From: Martin Karp <makarp@n173-p15.eduroam.kth.se>
Date: Fri, 20 Sep 2024 09:24:57 +0200
Subject: [PATCH 8/9] make lines shorter

---
 src/sem/map_1d.f90 | 27 ++++++++++++++++++---------
 src/sem/map_2d.f90 | 39 ++++++++++++++++++++++++++-------------
 2 files changed, 44 insertions(+), 22 deletions(-)

diff --git a/src/sem/map_1d.f90 b/src/sem/map_1d.f90
index 278cdd31fb6..95410d5ae08 100644
--- a/src/sem/map_1d.f90
+++ b/src/sem/map_1d.f90
@@ -21,9 +21,11 @@ module map_1d
   private
   !> Type that encapsulates a mapping from each gll point in the mesh
   !! to its corresponding (global) GLL point index in one direction.
-  !! @remark Could also be rather easily extended to say polar coordinates as well.
+  !! @remark Could also be rather easily extended to say polar coordinates
+  !! as well ( I think). Martin Karp
   type, public :: map_1d_t
-     !> Checks whether the specified direction is in the r,s, or t direction for each element.
+     !> Checks whether the specified direction is in the r,s, or t
+     !! direction for each element.
      integer, allocatable :: dir_el(:)
      !> Checks which level an element belongs to.
      integer, allocatable :: el_lvl(:)
@@ -80,7 +82,8 @@ subroutine map_1d_init(this, coef,  dir, tol)
 
     if (NEKO_BCKND_DEVICE .eq. 1) then
        if (pe_rank .eq. 0) then
-          call neko_warning('map_1d does not copy indices to device, but ok if used on cpu and for io')
+          call neko_warning('map_1d does not copy indices to device,'\\ &
+                            ' but ok if used on cpu and for io')
        end if
     end if
 
@@ -118,11 +121,14 @@ subroutine map_1d_init(this, coef,  dir, tol)
        !we assume elements are stacked on each other...
        ! Check which one of the normalized vectors are closest to dir
        ! If we want to incorporate other directions, we should look here
-       el_dim(1,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - this%msh%elements(i)%e%pts(2)%p%x)
+       el_dim(1,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
+                     this%msh%elements(i)%e%pts(2)%p%x)
        el_dim(1,:) = el_dim(1,:)/norm2(el_dim(1,:))
-       el_dim(2,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - this%msh%elements(i)%e%pts(3)%p%x)
+       el_dim(2,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
+                     this%msh%elements(i)%e%pts(3)%p%x)
        el_dim(2,:) = el_dim(2,:)/norm2(el_dim(2,:))
-       el_dim(3,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - this%msh%elements(i)%e%pts(5)%p%x)
+       el_dim(3,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
+                     this%msh%elements(i)%e%pts(5)%p%x)
        el_dim(3,:) = el_dim(3,:)/norm2(el_dim(3,:))
        ! Checks which directions in rst the xyz corresponds to
        ! 1 corresponds to r, 2 to s, 3 to t and are stored in dir_el
@@ -142,7 +148,8 @@ subroutine map_1d_init(this, coef,  dir, tol)
           if(this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
        end if
     end do
-    ! While loop where at each iteation the global maximum value propagates down one level.
+    ! While loop where at each iteation the global maximum value 
+    ! propagates down one level.
     ! When the minumum value has propagated to the highest level this stops.
     ! Only works when the bottom plate of the domain is flat.
     do while (.not. relcmp(glmax(min_vals,n), glb_min, this%tol))
@@ -313,7 +320,8 @@ subroutine map_1d_average_field_list(this, avg_planes, field_list)
     do j = 2, field_list%size() + 1
        do i = 1, n
           avg_planes%x(this%pt_lvl(i,1,1,1),j) = &
-          avg_planes%x(this%pt_lvl(i,1,1,1),j) + field_list%items(j-1)%ptr%x(i,1,1,1)*this%coef%B(i,1,1,1) &
+          avg_planes%x(this%pt_lvl(i,1,1,1),j) + &
+          field_list%items(j-1)%ptr%x(i,1,1,1)*this%coef%B(i,1,1,1) &
           /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
        end do
     end do
@@ -355,7 +363,8 @@ subroutine map_1d_average_vector_ptr(this, avg_planes, vector_ptr)
     do j = 2, size(vector_ptr)+1
        do i = 1, n
           avg_planes%x(this%pt_lvl(i,1,1,1),j) = &
-          avg_planes%x(this%pt_lvl(i,1,1,1),j) + vector_ptr(j-1)%ptr%x(i)*this%coef%B(i,1,1,1) &
+          avg_planes%x(this%pt_lvl(i,1,1,1),j) + &
+          vector_ptr(j-1)%ptr%x(i)*this%coef%B(i,1,1,1) &
           /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
        end do
     end do
diff --git a/src/sem/map_2d.f90 b/src/sem/map_2d.f90
index bb18fc7fb4b..8659899a9e3 100644
--- a/src/sem/map_2d.f90
+++ b/src/sem/map_2d.f90
@@ -17,7 +17,8 @@ module map_2d
   use vector, only: vector_ptr_t
   use logger, only: neko_log, LOG_SIZE
   use utils, only: neko_error, neko_warning
-  use math, only: glmax, glmin, glimax, glimin, relcmp, cmult, add2s1, col2, copy, rzero
+  use math, only: glmax, glmin, glimax, glimin, relcmp, cmult, &
+                  add2s1, col2, copy, rzero
   use neko_mpi_types
   use fld_file_data
   use field
@@ -123,9 +124,12 @@ subroutine map_2d_init(this,coef, dir, tol)
        !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(this%msh%elements(i)%e%pts(1)%p%x-this%msh%elements(i)%e%pts(2)%p%x)
-       el_dim(2,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-this%msh%elements(i)%e%pts(3)%p%x)
-       el_dim(3,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-this%msh%elements(i)%e%pts(5)%p%x)
+       el_dim(1,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-&
+                     this%msh%elements(i)%e%pts(2)%p%x)
+       el_dim(2,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-&
+                     this%msh%elements(i)%e%pts(3)%p%x)
+       el_dim(3,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-&
+                     this%msh%elements(i)%e%pts(5)%p%x)
        ! 1 corresponds to x, 2 to y, 3 to z
        el_h = el_dim(this%map_1d%dir_el(i),dir)
        this%el_heights%x(:,:,:,i) = el_h
@@ -136,7 +140,8 @@ subroutine map_2d_init(this,coef, dir, tol)
     call copy(this%u%x,this%el_heights%x,n)
     call copy(this%old_u%x,this%el_heights%x,n)
     call copy(this%avg_u%x,this%el_heights%x,n)
-    call perform_global_summation(this%u, this%avg_u, this%old_u, this%map_1d%n_el_lvls, &
+    call perform_global_summation(this%u, this%avg_u, this%old_u, &
+         this%map_1d%n_el_lvls, &
          this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, this%msh%nelv, lx)
     this%domain_height = this%u%x(1,1,1,1)
 
@@ -219,12 +224,15 @@ subroutine map_2d_average_field_list(this,fld_data2D,fld_data3D)
     this%avg_u = 0.0_rp
     do i = 1, fld_data3D%size()
        call copy(this%old_u%x,fld_data3D%items(i)%ptr%x,n)
-       call perform_local_summation(this%u,this%old_u, this%el_heights, this%domain_height, &
+       call perform_local_summation(this%u,this%old_u, &
+            this%el_heights, this%domain_height, &
             this%map_1d%dir_el, this%coef, this%msh%nelv, lx)
        call copy(this%old_u%x,this%u%x,n)
        call copy(this%avg_u%x,this%u%x,n)
-       call perform_global_summation(this%u, this%avg_u, this%old_u, this%map_1d%n_el_lvls, &
-            this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, this%msh%nelv, lx)
+       call perform_global_summation(this%u, this%avg_u, & 
+            this%old_u, this%map_1d%n_el_lvls, &
+            this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, &
+            this%msh%nelv, lx)
        call copy(fld_data3D%items(i)%ptr%x,this%u%x,n)
     end do
     call fld_data2D%get_list(fields2d,fld_data2D%size())
@@ -293,12 +301,15 @@ subroutine map_2d_average(this,fld_data2D,fld_data3D)
     this%avg_u = 0.0_rp
     do i = 1, fld_data3D%size()
        call copy(this%old_u%x,fields3D(i)%ptr%x,n)
-       call perform_local_summation(this%u,this%old_u, this%el_heights, this%domain_height, &
+       call perform_local_summation(this%u,this%old_u,&
+            this%el_heights, this%domain_height, &
             this%map_1d%dir_el, this%coef, this%msh%nelv, lx)
        call copy(this%old_u%x,this%u%x,n)
        call copy(this%avg_u%x,this%u%x,n)
-       call perform_global_summation(this%u, this%avg_u, this%old_u, this%map_1d%n_el_lvls, &
-            this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, this%msh%nelv, lx)
+       call perform_global_summation(this%u, this%avg_u, &
+            this%old_u, this%map_1d%n_el_lvls, &
+            this%map_1d%dir_el,this%coef%gs_h,&
+            this%coef%mult, this%msh%nelv, lx)
        call copy(fields3D(i)%ptr%x,this%u%x,n)
     end do
     call fld_data2D%get_list(fields2d,fld_data2D%size())
@@ -309,7 +320,8 @@ subroutine map_2d_average(this,fld_data2D,fld_data3D)
     end do
   end subroutine map_2d_average
  
-  subroutine perform_global_summation(u, avg_u, old_u, n_levels, hom_dir_el, gs_h, mult, nelv,lx)
+  subroutine perform_global_summation(u, avg_u, old_u, n_levels, &
+                                     hom_dir_el, gs_h, mult, nelv, lx)
     type(field_t), intent(inout) :: u, avg_u, old_u
     type(gs_t), intent(inout) :: gs_h
     integer, intent(in) :: n_levels, nelv, lx
@@ -364,7 +376,8 @@ subroutine perform_global_summation(u, avg_u, old_u, n_levels, hom_dir_el, gs_h,
     end do
   end subroutine perform_global_summation
   
-  subroutine perform_local_summation(u_out, u, el_heights,domain_height, hom_dir_el, coef, nelv,lx)
+  subroutine perform_local_summation(u_out, u, el_heights, domain_height,&
+                                     hom_dir_el, coef, nelv, lx)
     type(field_t), intent(inout) :: u, u_out, el_heights
     type(coef_t), intent(inout) :: coef
     integer, intent(in) :: nelv, lx

From e6b04914dd3ef1885559674e31775e77c77224cc Mon Sep 17 00:00:00 2001
From: Martin Karp <makarp@n173-p15.eduroam.kth.se>
Date: Fri, 20 Sep 2024 09:35:17 +0200
Subject: [PATCH 9/9] wrong slash

---
 src/sem/map_1d.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/sem/map_1d.f90 b/src/sem/map_1d.f90
index 95410d5ae08..c84708fcd69 100644
--- a/src/sem/map_1d.f90
+++ b/src/sem/map_1d.f90
@@ -82,7 +82,7 @@ subroutine map_1d_init(this, coef,  dir, tol)
 
     if (NEKO_BCKND_DEVICE .eq. 1) then
        if (pe_rank .eq. 0) then
-          call neko_warning('map_1d does not copy indices to device,'\\ &
+          call neko_warning('map_1d does not copy indices to device,'// &
                             ' but ok if used on cpu and for io')
        end if
     end if