25 use,
intrinsic :: iso_c_binding
30 integer :: nelv_2d = 0
31 integer :: glb_nelv_2d = 0
32 integer :: offset_el_2d = 0
35 integer,
allocatable :: idx_2d(:)
36 integer,
allocatable :: el_idx_2d(:)
40 type(
coef_t),
pointer :: coef => null()
46 real(kind=
rp) :: domain_height
50 generic :: init => init_int, init_char
53 generic :: average => average_list, average_file
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)
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)
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
86 call mpi_allreduce(this%nelv_2d, this%glb_nelv_2d, 1, &
90 call mpi_exscan(this%nelv_2d, this%offset_el_2d, 1, &
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
96 this%n_2d = this%nelv_2d*lxy
97 allocate(this%idx_2d(this%n_2d))
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
104 this%idx_2d(j*lxy+i) = linear_index(1,i,1,e,lx,lx,lx)
107 if (this%map_1d%dir_el(e) .eq. 2)
then
110 this%idx_2d(j*lxy+k+lx*(i-1)) = linear_index(k,1,i,e,lx,lx,lx)
114 if (this%map_1d%dir_el(e) .eq. 3)
then
116 this%idx_2d(j*lxy+i) = linear_index(i,1,1,e,lx,lx,lx)
122 do i = 1, this%msh%nelv
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)
134 el_h = el_dim(this%map_1d%dir_el(i),dir)
135 this%el_heights%x(:,:,:,i) = el_h
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)
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)
152 type(
coef_t),
intent(inout),
target :: coef
153 character(len=*),
intent(in) :: dir
154 real(kind=
rp),
intent(in) :: tol
157 if (trim(dir) .eq.
'x')
then
159 else if (trim(dir) .eq.
'y')
then
161 else if (trim(dir) .eq.
'z')
then
164 call neko_error(
'Direction not supported, map_2d')
167 call this%init(coef,idir,tol)
177 class(
map_2d_t),
intent(inout) :: this
180 real(kind=
rp),
pointer,
dimension(:,:,:,:) :: x_ptr, y_ptr
184 integer :: i, j, lx, lxy, e
186 call fld_data2d%init(this%nelv_2d,this%offset_el_2d)
188 fld_data2d%lx = this%dof%Xh%lx
189 fld_data2d%ly = this%dof%Xh%ly
192 fld_data2d%glb_nelv = this%glb_nelv_2d
193 lxy = fld_data2d%lx*fld_data2d%ly
194 n_2d = lxy*this%nelv_2d
196 call fld_data2d%x%init(n_2d)
197 call fld_data2d%y%init(n_2d)
198 allocate(fld_data2d%idx(this%nelv_2d))
200 if (this%dir .eq. 1)
then
204 if (this%dir .eq. 2)
then
208 if (this%dir .eq. 3)
then
212 do j = 1, this%nelv_2d
213 fld_data2d%idx(j) = this%el_idx_2d(j)
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)
219 allocate(fields2d(fld_data3d%size()))
221 call fld_data2d%init_n_fields(fld_data3d%size(),n_2d)
225 do i = 1, fld_data3d%size()
226 call copy(this%old_u%x,fld_data3d%items(i)%ptr%x,n)
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)
233 this%old_u, this%map_1d%n_el_lvls, &
234 this%map_1d%dir_el,this%coef%gs_h, this%coef%mult, &
236 call copy(fld_data3d%items(i)%ptr%x,this%u%x,n)
238 call fld_data2d%get_list(fields2d,fld_data2d%size())
239 do i = 1, fld_data3d%size()
241 fields2d(i)%ptr%x(j) = fld_data3d%items(i)%ptr%x(this%idx_2d(j),1,1,1)
253 class(
map_2d_t),
intent(inout) :: this
256 real(kind=
rp),
pointer,
dimension(:,:,:,:) :: x_ptr, y_ptr
258 type(
vector_ptr_t),
allocatable :: fields3d(:), fields2d(:)
260 integer :: i, j, lx, lxy, e
262 call fld_data2d%init(this%nelv_2d,this%offset_el_2d)
264 fld_data2d%lx = this%dof%Xh%lx
265 fld_data2d%ly = this%dof%Xh%ly
268 fld_data2d%glb_nelv = this%glb_nelv_2d
269 lxy = fld_data2d%lx*fld_data2d%ly
270 n_2d = lxy*this%nelv_2d
272 call fld_data2d%x%init(n_2d)
273 call fld_data2d%y%init(n_2d)
274 allocate(fld_data2d%idx(n_2d))
276 if (this%dir .eq. 1)
then
280 if (this%dir .eq. 2)
then
284 if (this%dir .eq. 3)
then
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)
293 allocate(fields3d(fld_data3d%size()))
294 allocate(fields2d(fld_data3d%size()))
296 call fld_data2d%init_n_fields(fld_data3d%size(),n_2d)
297 call fld_data3d%get_list(fields3d,fld_data3d%size())
302 do i = 1, fld_data3d%size()
303 call copy(this%old_u%x,fields3d(i)%ptr%x,n)
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)
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)
315 call fld_data2d%get_list(fields2d,fld_data2d%size())
316 do i = 1, fld_data3d%size()
318 fields2d(i)%ptr%x(j) = fields3d(i)%ptr%x(this%idx_2d(j))
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
338 if (neko_bcknd_device .eq. 1) &
341 call gs_h%op(u,gs_op_add)
342 if (neko_bcknd_device .eq. 1) &
344 call col2(u%x,mult,n)
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)
361 old_u%x(:,:,:,e) = u%x(:,:,:,e)
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)
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
387 integer :: n, i, j, e, k
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)
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)
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)
Copy data between host and device (or device and device)
type(mpi_comm) neko_comm
MPI communicator.
Device abstraction, common interface for various accelerators.
integer, parameter, public host_to_device
integer, parameter, public device_to_host
Defines a mapping of the degrees of freedom.
Simple module to handle fld file series. Provides an interface to the different fields sotred in a fl...
type(log_t), public neko_log
Global log stream.
integer, parameter, public log_size
Creates a 1d GLL point map along a specified direction based on the connectivity in the mesh.
Maps a 3D dofmap to a 2D spectral element grid.
subroutine perform_local_summation(u_out, u, el_heights, domain_height, hom_dir_el, coef, nelv, lx)
subroutine map_2d_average(this, fld_data2D, fld_data3D)
Computes average if field list in one direction and outputs 2D field with averaged values.
subroutine perform_global_summation(u, avg_u, old_u, n_levels, hom_dir_el, gs_h, mult, nelv, lx)
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.
subroutine map_2d_init(this, coef, dir, tol)
subroutine map_2d_init_char(this, coef, dir, tol)
subroutine, public cmult(a, c, n)
Multiplication by constant c .
subroutine, public add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
integer function, public glimin(a, n)
Min of an integer vector of length n.
subroutine, public col2(a, b, n)
Vector multiplication .
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
integer function, public glimax(a, n)
Max of an integer vector of length n.
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
integer, parameter, public rp
Global precision used in computations.
Defines a function space.
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
field_list_t, To be able to group fields together
Type that encapsulates a mapping from each gll point in the mesh to its corresponding (global) GLL po...
The function space for the SEM solution fields.