Neko  0.9.0
A portable framework for high-order spectral element flow simulations
map_2d.f90
Go to the documentation of this file.
1 
2 
3 module map_2d
4  use num_types, only: rp
5  use space, only: space_t
6  use dofmap, only: dofmap_t
7  use map_1d
9  use mesh, only: mesh_t
10  use device
11  use utils
12  use comm
13  use field_list
14  use coefs, only: coef_t
15  use field_list, only: field_list_t
16  use matrix, only: matrix_t
17  use vector, only: vector_ptr_t
18  use logger, only: neko_log, log_size
19  use utils, only: neko_error, neko_warning
20  use math, only: glmax, glmin, glimax, glimin, relcmp, cmult, &
22  use neko_mpi_types
23  use fld_file_data
24  use field
25  use, intrinsic :: iso_c_binding
26  implicit none
27  private
28 
29  type, public :: map_2d_t
30  integer :: nelv_2d = 0
31  integer :: glb_nelv_2d = 0
32  integer :: offset_el_2d = 0
33  integer :: lxy = 0
34  integer :: n_2d = 0
35  integer, allocatable :: idx_2d(:)
36  integer, allocatable :: el_idx_2d(:)
37  type(map_1d_t) :: map_1d
38  type(mesh_t), pointer :: msh
39  type(dofmap_t), pointer :: dof => null()
40  type(coef_t), pointer :: coef => null()
41  type(field_t) :: u
42  type(field_t) :: old_u
43  type(field_t) :: avg_u
44  type(field_t) :: el_heights
45  integer :: dir
46  real(kind=rp) :: domain_height
47  contains
48  procedure, pass(this) :: init_int => map_2d_init
49  procedure, pass(this) :: init_char => map_2d_init_char
50  generic :: init => init_int, init_char
51  procedure, pass(this) :: average_file => map_2d_average
52  procedure, pass(this) :: average_list => map_2d_average_field_list
53  generic :: average => average_list, average_file
54  end type map_2d_t
55 
56 contains
57 
58  subroutine map_2d_init(this,coef, dir, tol)
59  class(map_2d_t), intent(inout) :: this
60  type(coef_t), intent(inout), target :: coef
61  integer, intent(in) :: dir
62  real(kind=rp), intent(in) :: tol
63  real(kind=rp) :: el_dim(3,3), el_h
64  integer :: i, e, j, ierr, k, lx, lxy, n
65  call this%map_1d%init(coef,dir,tol)
66  this%msh => coef%msh
67  this%coef => coef
68  this%dof => coef%dof
69 
70  call this%u%init(this%dof)
71  call this%old_u%init(this%dof)
72  call this%avg_u%init(this%dof)
73  call this%el_heights%init(this%dof)
74  this%dir = dir
75 
76  n = this%dof%size()
77  lx = this%dof%Xh%lx
78  lxy = this%dof%Xh%lxy
79  this%nelv_2d = 0
80  do i = 1, this%msh%nelv
81  if (this%map_1d%el_lvl(i) .eq. 1) then
82  this%nelv_2d = this%nelv_2d + 1
83  end if
84  end do
85  this%glb_nelv_2d = 0
86  call mpi_allreduce(this%nelv_2d, this%glb_nelv_2d, 1, &
87  mpi_integer, mpi_sum, neko_comm, ierr)
88 
89  this%offset_el_2d = 0
90  call mpi_exscan(this%nelv_2d, this%offset_el_2d, 1, &
91  mpi_integer, mpi_sum, neko_comm, ierr)
92  allocate(this%el_idx_2d(this%nelv_2d))
93  do i = 1, this%nelv_2d
94  this%el_idx_2d(i) = this%offset_el_2d + i
95  end do
96  this%n_2d = this%nelv_2d*lxy
97  allocate(this%idx_2d(this%n_2d))
98 
99  j = 0
100  do e = 1, this%msh%nelv
101  if (this%map_1d%el_lvl(e) .eq. 1) then
102  if (this%map_1d%dir_el(e) .eq. 1) then
103  do i = 1, lxy
104  this%idx_2d(j*lxy+i) = linear_index(1,i,1,e,lx,lx,lx)
105  end do
106  end if
107  if (this%map_1d%dir_el(e) .eq. 2) then
108  do i = 1, lx
109  do k = 1, lx
110  this%idx_2d(j*lxy+k+lx*(i-1)) = linear_index(k,1,i,e,lx,lx,lx)
111  end do
112  end do
113  end if
114  if (this%map_1d%dir_el(e) .eq. 3) then
115  do i = 1, lxy
116  this%idx_2d(j*lxy+i) = linear_index(i,1,1,e,lx,lx,lx)
117  end do
118  end if
119  j = j +1
120  end if
121  end do
122  do i = 1, this%msh%nelv
123  !find height in hom-dir
124  !direction in local coords (r,s,t) that is hom is stored in map_1d%dir_el
125  !set element to height
126  !we assume elements are stacked on eachother...
127  el_dim(1,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-&
128  this%msh%elements(i)%e%pts(2)%p%x)
129  el_dim(2,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-&
130  this%msh%elements(i)%e%pts(3)%p%x)
131  el_dim(3,:) = abs(this%msh%elements(i)%e%pts(1)%p%x-&
132  this%msh%elements(i)%e%pts(5)%p%x)
133  ! 1 corresponds to x, 2 to y, 3 to z
134  el_h = el_dim(this%map_1d%dir_el(i),dir)
135  this%el_heights%x(:,:,:,i) = el_h
136  end do
137  !Need to compute mapping for 3d to 2d
138  !Does order matter? Think its ok as long as values written in same order
139 
140  call copy(this%u%x,this%el_heights%x,n)
141  call copy(this%old_u%x,this%el_heights%x,n)
142  call copy(this%avg_u%x,this%el_heights%x,n)
143  call perform_global_summation(this%u, this%avg_u, this%old_u, &
144  this%map_1d%n_el_lvls, &
145  this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, this%msh%nelv, lx)
146  this%domain_height = this%u%x(1,1,1,1)
147 
148  end subroutine map_2d_init
149 
150  subroutine map_2d_init_char(this, coef, dir, tol)
151  class(map_2d_t) :: this
152  type(coef_t), intent(inout), target :: coef
153  character(len=*), intent(in) :: dir
154  real(kind=rp), intent(in) :: tol
155  integer :: idir
156 
157  if (trim(dir) .eq. 'x') then
158  idir = 1
159  else if (trim(dir) .eq. 'y') then
160  idir = 2
161  else if (trim(dir) .eq. 'z') then
162  idir = 3
163  else
164  call neko_error('Direction not supported, map_2d')
165  end if
166 
167  call this%init(coef,idir,tol)
168 
169  end subroutine map_2d_init_char
170 
171 
176  subroutine map_2d_average_field_list(this,fld_data2D,fld_data3D)
177  class(map_2d_t), intent(inout) :: this
178  type(fld_file_data_t), intent(inout) :: fld_data2D
179  type(field_list_t), intent(inout) :: fld_data3D
180  real(kind=rp), pointer, dimension(:,:,:,:) :: x_ptr, y_ptr
181 
182  type(vector_ptr_t), allocatable :: fields2d(:)
183  integer :: n_2d, n
184  integer :: i, j, lx, lxy, e
185 
186  call fld_data2d%init(this%nelv_2d,this%offset_el_2d)
187  fld_data2d%gdim = 2
188  fld_data2d%lx = this%dof%Xh%lx
189  fld_data2d%ly = this%dof%Xh%ly
190  fld_data2d%lz = 1
191  lx = this%dof%Xh%lx
192  fld_data2d%glb_nelv = this%glb_nelv_2d
193  lxy = fld_data2d%lx*fld_data2d%ly
194  n_2d = lxy*this%nelv_2d
195  n = this%dof%size()
196  call fld_data2d%x%init(n_2d)
197  call fld_data2d%y%init(n_2d)
198  allocate(fld_data2d%idx(this%nelv_2d))
199 
200  if (this%dir .eq. 1) then
201  x_ptr => this%dof%z
202  y_ptr => this%dof%y
203  end if
204  if (this%dir .eq. 2) then
205  x_ptr => this%dof%x
206  y_ptr => this%dof%z
207  end if
208  if (this%dir .eq. 3) then
209  x_ptr => this%dof%x
210  y_ptr => this%dof%y
211  end if
212  do j = 1, this%nelv_2d
213  fld_data2d%idx(j) = this%el_idx_2d(j)
214  end do
215  do j = 1, n_2d
216  fld_data2d%x%x(j) = x_ptr(this%idx_2d(j),1,1,1)
217  fld_data2d%y%x(j) = y_ptr(this%idx_2d(j),1,1,1)
218  end do
219  allocate(fields2d(fld_data3d%size()))
220 
221  call fld_data2d%init_n_fields(fld_data3d%size(),n_2d)
222  this%u = 0.0_rp
223  this%old_u = 0.0_rp
224  this%avg_u = 0.0_rp
225  do i = 1, fld_data3d%size()
226  call copy(this%old_u%x,fld_data3d%items(i)%ptr%x,n)
227  call perform_local_summation(this%u,this%old_u, &
228  this%el_heights, this%domain_height, &
229  this%map_1d%dir_el, this%coef, this%msh%nelv, lx)
230  call copy(this%old_u%x,this%u%x,n)
231  call copy(this%avg_u%x,this%u%x,n)
232  call perform_global_summation(this%u, this%avg_u, &
233  this%old_u, this%map_1d%n_el_lvls, &
234  this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, &
235  this%msh%nelv, lx)
236  call copy(fld_data3d%items(i)%ptr%x,this%u%x,n)
237  end do
238  call fld_data2d%get_list(fields2d,fld_data2d%size())
239  do i = 1, fld_data3d%size()
240  do j = 1, n_2d
241  fields2d(i)%ptr%x(j) = fld_data3d%items(i)%ptr%x(this%idx_2d(j),1,1,1)
242  end do
243  end do
244 
245  end subroutine map_2d_average_field_list
246 
247 
252  subroutine map_2d_average(this,fld_data2D,fld_data3D)
253  class(map_2d_t), intent(inout) :: this
254  type(fld_file_data_t), intent(inout) :: fld_data2D
255  type(fld_file_data_t), intent(inout) :: fld_data3D
256  real(kind=rp), pointer, dimension(:,:,:,:) :: x_ptr, y_ptr
257 
258  type(vector_ptr_t), allocatable :: fields3d(:), fields2d(:)
259  integer :: n_2d, n
260  integer :: i, j, lx, lxy, e
261 
262  call fld_data2d%init(this%nelv_2d,this%offset_el_2d)
263  fld_data2d%gdim = 2
264  fld_data2d%lx = this%dof%Xh%lx
265  fld_data2d%ly = this%dof%Xh%ly
266  fld_data2d%lz = 1
267  lx = this%dof%Xh%lx
268  fld_data2d%glb_nelv = this%glb_nelv_2d
269  lxy = fld_data2d%lx*fld_data2d%ly
270  n_2d = lxy*this%nelv_2d
271  n = this%dof%size()
272  call fld_data2d%x%init(n_2d)
273  call fld_data2d%y%init(n_2d)
274  allocate(fld_data2d%idx(n_2d))
275 
276  if (this%dir .eq. 1) then
277  x_ptr => this%dof%z
278  y_ptr => this%dof%y
279  end if
280  if (this%dir .eq. 2) then
281  x_ptr => this%dof%x
282  y_ptr => this%dof%z
283  end if
284  if (this%dir .eq. 3) then
285  x_ptr => this%dof%x
286  y_ptr => this%dof%y
287  end if
288  do j = 1, n_2d
289  fld_data2d%idx(j) = this%idx_2d(j)
290  fld_data2d%x%x(j) = x_ptr(this%idx_2d(j),1,1,1)
291  fld_data2d%y%x(j) = y_ptr(this%idx_2d(j),1,1,1)
292  end do
293  allocate(fields3d(fld_data3d%size()))
294  allocate(fields2d(fld_data3d%size()))
295 
296  call fld_data2d%init_n_fields(fld_data3d%size(),n_2d)
297  call fld_data3d%get_list(fields3d,fld_data3d%size())
298 
299  this%u = 0.0_rp
300  this%old_u = 0.0_rp
301  this%avg_u = 0.0_rp
302  do i = 1, fld_data3d%size()
303  call copy(this%old_u%x,fields3d(i)%ptr%x,n)
304  call perform_local_summation(this%u,this%old_u,&
305  this%el_heights, this%domain_height, &
306  this%map_1d%dir_el, this%coef, this%msh%nelv, lx)
307  call copy(this%old_u%x,this%u%x,n)
308  call copy(this%avg_u%x,this%u%x,n)
309  call perform_global_summation(this%u, this%avg_u, &
310  this%old_u, this%map_1d%n_el_lvls, &
311  this%map_1d%dir_el,this%coef%gs_h,&
312  this%coef%mult, this%msh%nelv, lx)
313  call copy(fields3d(i)%ptr%x,this%u%x,n)
314  end do
315  call fld_data2d%get_list(fields2d,fld_data2d%size())
316  do i = 1, fld_data3d%size()
317  do j = 1, n_2d
318  fields2d(i)%ptr%x(j) = fields3d(i)%ptr%x(this%idx_2d(j))
319  end do
320  end do
321  end subroutine map_2d_average
322 
323  subroutine perform_global_summation(u, avg_u, old_u, n_levels, &
324  hom_dir_el, gs_h, mult, nelv, lx)
325  type(field_t), intent(inout) :: u, avg_u, old_u
326  type(gs_t), intent(inout) :: gs_h
327  integer, intent(in) :: n_levels, nelv, lx
328  integer, intent(in) :: hom_dir_el(nelv)
329  real(kind=rp), intent(in) :: mult(nelv*lx**3)
330  real(kind=rp) :: temp_el(lx,lx,lx)
331  integer :: n, i, j, e, k
332  type(c_ptr) :: ptr
333 
334  n = u%dof%size()
335 
336  do i = 1, n_levels-1
337  !compute average
338  if (neko_bcknd_device .eq. 1) &
339  call device_memcpy(u%x, u%x_d, n, &
340  host_to_device, sync=.true.)
341  call gs_h%op(u,gs_op_add)
342  if (neko_bcknd_device .eq. 1) &
343  call device_memcpy(u%x, u%x_d, n, device_to_host, sync=.true.)
344  call col2(u%x,mult,n)
345  !Assumes sugarcube I think
346  do e = 1, nelv
347  temp_el = 2.0*u%x(:,:,:,e)-old_u%x(:,:,:,e)
348  if (hom_dir_el(e) .eq. 1) then
349  u%x(1,:,:,e) = temp_el(lx,:,:)
350  avg_u%x(1,:,:,e) = avg_u%x(1,:,:,e)+temp_el(1,:,:)
351  u%x(lx,:,:,e) = temp_el(1,:,:)
352  else if (hom_dir_el(e) .eq. 2) then
353  u%x(:,1,:,e) = temp_el(:,lx,:)
354  avg_u%x(:,1,:,e) = avg_u%x(:,1,:,e)+temp_el(:,1,:)
355  u%x(:,lx,:,e) = temp_el(:,1,:)
356  else if (hom_dir_el(e) .eq. 3) then
357  u%x(:,:,1,e) = temp_el(:,:,lx)
358  avg_u%x(:,:,1,e) = avg_u%x(:,:,1,e)+temp_el(:,:,1)
359  u%x(:,:,lx,e) = temp_el(:,:,1)
360  end if
361  old_u%x(:,:,:,e) = u%x(:,:,:,e)
362  end do
363  end do
364  do e = 1, nelv
365  do i = 1,lx
366  do j = 1, lx
367  if (hom_dir_el(e) .eq. 1) then
368  u%x(:,i,j,e) = avg_u%x(1,i,j,e)
369  else if (hom_dir_el(e) .eq. 2) then
370  u%x(i,:,j,e) = avg_u%x(i,1,j,e)
371  else if (hom_dir_el(e) .eq. 3) then
372  u%x(i,j,:,e) = avg_u%x(i,j,1,e)
373  end if
374  end do
375  end do
376  end do
377  end subroutine perform_global_summation
378 
379  subroutine perform_local_summation(u_out, u, el_heights, domain_height,&
380  hom_dir_el, coef, nelv, lx)
381  type(field_t), intent(inout) :: u, u_out, el_heights
382  type(coef_t), intent(inout) :: coef
383  integer, intent(in) :: nelv, lx
384  integer, intent(in) :: hom_dir_el(nelv)
385  real(kind=rp), intent(in) :: domain_height
386  real(kind=rp) :: wt
387  integer :: n, i, j, e, k
388 
389  n = nelv*lx**3
390 
391  call col2(u%x,el_heights%x,n)
392  call cmult(u%x, 1.0_rp/(2.0*domain_height),n)
393  call rzero(u_out%x,n)
394 
395 
396  do e = 1, nelv
397  do i = 1,lx
398  do j = 1, lx
399  do k = 1, lx
400  wt = coef%Xh%wx(k)
401  if (hom_dir_el(e) .eq. 1) then
402  u_out%x(1,i,j,e) = u_out%x(1,i,j,e)+wt*u%x(k,i,j,e)
403  else if (hom_dir_el(e) .eq. 2) then
404  u_out%x(i,1,j,e) = u_out%x(i,1,j,e)+wt*u%x(i,k,j,e)
405  else if (hom_dir_el(e) .eq. 3) then
406  u_out%x(i,j,1,e) = u_out%x(i,j,1,e)+wt*u%x(i,j,k,e)
407  end if
408  end do
409  end do
410  end do
411  end do
412 
413  do e = 1, nelv
414  do i = 1,lx
415  do j = 1, lx
416  if (hom_dir_el(e) .eq. 1) then
417  u_out%x(:,i,j,e) = u_out%x(1,i,j,e)
418  else if (hom_dir_el(e) .eq. 2) then
419  u_out%x(i,:,j,e) = u_out%x(i,1,j,e)
420  else if (hom_dir_el(e) .eq. 3) then
421  u_out%x(i,j,:,e) = u_out%x(i,j,1,e)
422  end if
423  end do
424  end do
425  end do
426  end subroutine perform_local_summation
427 
428 end module map_2d
Copy data between host and device (or device and device)
Definition: device.F90:51
Coefficients.
Definition: coef.f90:34
Definition: comm.F90:1
type(mpi_comm) neko_comm
MPI communicator.
Definition: comm.F90:16
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
integer, parameter, public host_to_device
Definition: device.F90:47
integer, parameter, public device_to_host
Definition: device.F90:47
Defines a mapping of the degrees of freedom.
Definition: dofmap.f90:35
Defines a field.
Definition: field.f90:34
Simple module to handle fld file series. Provides an interface to the different fields sotred in a fl...
Gather-scatter.
Logging routines.
Definition: log.f90:34
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
integer, parameter, public log_size
Definition: log.f90:42
Creates a 1d GLL point map along a specified direction based on the connectivity in the mesh.
Definition: map_1d.f90:3
Maps a 3D dofmap to a 2D spectral element grid.
Definition: map_2d.f90:3
subroutine perform_local_summation(u_out, u, el_heights, domain_height, hom_dir_el, coef, nelv, lx)
Definition: map_2d.f90:381
subroutine map_2d_average(this, fld_data2D, fld_data3D)
Computes average if field list in one direction and outputs 2D field with averaged values.
Definition: map_2d.f90:253
subroutine perform_global_summation(u, avg_u, old_u, n_levels, hom_dir_el, gs_h, mult, nelv, lx)
Definition: map_2d.f90:325
subroutine map_2d_average_field_list(this, fld_data2D, fld_data3D)
Computes average if field list in one direction and outputs 2D field with averaged values.
Definition: map_2d.f90:177
subroutine map_2d_init(this, coef, dir, tol)
Definition: map_2d.f90:59
subroutine map_2d_init_char(this, coef, dir, tol)
Definition: map_2d.f90:151
Definition: math.f90:60
subroutine, public cmult(a, c, n)
Multiplication by constant c .
Definition: math.f90:311
subroutine, public add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
Definition: math.f90:658
integer function, public glimin(a, n)
Min of an integer vector of length n.
Definition: math.f90:422
subroutine, public col2(a, b, n)
Vector multiplication .
Definition: math.f90:729
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition: math.f90:377
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
subroutine, public rzero(a, n)
Zero a real vector.
Definition: math.f90:195
integer function, public glimax(a, n)
Max of an integer vector of length n.
Definition: math.f90:392
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition: math.f90:407
Defines a matrix.
Definition: matrix.f90:34
Defines a mesh.
Definition: mesh.f90:34
MPI derived types.
Definition: mpi_types.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Defines a function space.
Definition: space.f90:34
Utilities.
Definition: utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition: utils.f90:245
Defines a vector.
Definition: vector.f90:34
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
field_list_t, To be able to group fields together
Definition: field_list.f90:13
Type that encapsulates a mapping from each gll point in the mesh to its corresponding (global) GLL po...
Definition: map_1d.f90:26
The function space for the SEM solution fields.
Definition: space.f90:62