Neko 0.9.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
map_2d.f90
Go to the documentation of this file.
1
2
3module 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, &
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(:)
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
56contains
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
428end 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 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_init(this, coef, dir, tol)
Definition map_2d.f90:59
subroutine map_2d_init_char(this, coef, dir, tol)
Definition map_2d.f90:151
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_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
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
pure integer function, public linear_index(i, j, k, l, lx, ly, lz)
Compute the address of a (i,j,k,l) array with sizes (1:lx, 1:ly, 1:lz, :)
Definition utils.f90:174
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
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