72 type(
coef_t),
intent(inout),
target :: coef
73 integer,
intent(in) :: dir
74 real(kind=
rp),
intent(in) :: tol
75 integer :: nelv, lx, n, i, e, lvl, ierr
76 real(kind=
rp),
contiguous,
pointer :: line(:,:,:,:)
77 real(kind=
rp),
allocatable :: min_vals(:,:,:,:)
78 real(kind=
rp),
allocatable :: min_temp(:,:,:,:)
79 type(c_ptr) :: min_vals_d = c_null_ptr
80 real(kind=
rp) :: el_dim(3,3), glb_min, glb_max, el_min
86 call neko_warning(
'map_1d does not copy indices to device,'// &
87 ' but ok if used on cpu and for io')
101 else if (dir .eq. 2)
then
103 else if(dir .eq. 3)
then
106 call neko_error(
'Invalid dir for geopmetric comm')
108 allocate(this%dir_el(nelv))
109 allocate(this%el_lvl(nelv))
110 allocate(this%pt_lvl(lx, lx, lx, nelv))
111 allocate(min_vals(lx, lx, lx, nelv))
112 allocate(min_temp(lx, lx, lx, nelv))
125 el_dim(1,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
126 this%msh%elements(i)%e%pts(2)%p%x)
127 el_dim(1,:) = el_dim(1,:)/norm2(el_dim(1,:))
128 el_dim(2,:) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
129 this%msh%elements(i)%e%pts(3)%p%x)
130 el_dim(2,:) = el_dim(2,:)/norm2(el_dim(2,:))
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 el_dim(3,:) = el_dim(3,:)/norm2(el_dim(3,:))
136 this%dir_el(i) = maxloc(el_dim(:,this%dir),dim=1)
138 glb_min =
glmin(line,n)
139 glb_max =
glmax(line,n)
145 el_min = minval(line(:,:,:,e))
146 min_vals(:,:,:,e) = el_min
148 if (
relcmp(el_min, glb_min, this%tol))
then
149 if(this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
156 do while (.not.
relcmp(
glmax(min_vals,n), glb_min, this%tol))
160 if (this%dir_el(e) .eq. 1)
then
161 if (line(1,1,1,e) .gt. line(lx,1,1,e))
then
162 min_vals(lx,:,:,e) = glb_max
164 min_vals(1,:,:,e) = glb_max
167 if (this%dir_el(e) .eq. 2)
then
168 if (line(1,1,1,e) .gt. line(1,lx,1,e))
then
169 min_vals(:,lx,:,e) = glb_max
171 min_vals(:,1,:,e) = glb_max
174 if (this%dir_el(e) .eq. 3)
then
175 if (line(1,1,1,e) .gt. line(1,1,lx,e))
then
176 min_vals(:,:,lx,e) = glb_max
178 min_vals(:,:,1,e) = glb_max
188 call coef%gs_h%op(min_vals, n, gs_op_add)
194 call col2(min_vals, coef%mult, n)
195 call cmult(min_temp, -1.0_rp, n)
196 call add2s1(min_vals, min_temp, 2.0_rp, n)
204 el_min = minval(min_vals(:,:,:,e))
205 min_vals(:,:,:,e) = el_min
206 if (
relcmp(el_min, glb_min, this%tol))
then
207 if (this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
211 this%n_el_lvls =
glimax(this%el_lvl,nelv)
212 this%n_gll_lvls = this%n_el_lvls*lx
217 lvl = lx * (this%el_lvl(e) - 1) + i
218 if (this%dir_el(e) .eq. 1)
then
219 if (line(1,1,1,e) .gt. line(lx,1,1,e))
then
220 this%pt_lvl(lx-i+1,:,:,e) = lvl
222 this%pt_lvl(i,:,:,e) = lvl
225 if (this%dir_el(e) .eq. 2)
then
226 if (line(1,1,1,e) .gt. line(1,lx,1,e))
then
227 this%pt_lvl(:,lx-i+1,:,e) = lvl
229 this%pt_lvl(:,i,:,e) = lvl
232 if (this%dir_el(e) .eq. 3)
then
233 if (line(1,1,1,e) .gt. line(1,1,lx,e))
then
234 this%pt_lvl(:,:,lx-i+1,e) = lvl
236 this%pt_lvl(:,:,i,e) = lvl
241 if(
allocated(min_vals))
deallocate(min_vals)
242 if(c_associated(min_vals_d))
call device_free(min_vals_d)
243 if(
allocated(min_temp))
deallocate(min_temp)
244 allocate(this%volume_per_gll_lvl(this%n_gll_lvls))
245 this%volume_per_gll_lvl =0.0_rp
247 this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1)) = &
248 this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1)) + coef%B(i,1,1,1)
250 call mpi_allreduce(mpi_in_place,this%volume_per_gll_lvl, this%n_gll_lvls, &
299 class(
map_1d_t),
intent(inout) :: this
301 type(
matrix_t),
intent(inout) :: avg_planes
302 integer :: n, ierr, j, i
303 real(kind=
rp) :: coord
304 call avg_planes%free()
305 call avg_planes%init(this%n_gll_lvls,
field_list%size()+1)
310 if (this%dir .eq. 1) coord = this%dof%x(i,1,1,1)
311 if (this%dir .eq. 2) coord = this%dof%y(i,1,1,1)
312 if (this%dir .eq. 3) coord = this%dof%z(i,1,1,1)
313 avg_planes%x(this%pt_lvl(i,1,1,1),1) = &
314 avg_planes%x(this%pt_lvl(i,1,1,1),1) + coord*this%coef%B(i,1,1,1) &
315 /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
319 avg_planes%x(this%pt_lvl(i,1,1,1),j) = &
320 avg_planes%x(this%pt_lvl(i,1,1,1),j) + &
321 field_list%items(j-1)%ptr%x(i,1,1,1)*this%coef%B(i,1,1,1) &
322 /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
326 call mpi_allreduce(mpi_in_place, avg_planes%x, (
field_list%size()+1)*this%n_gll_lvls, &
339 class(
map_1d_t),
intent(inout) :: this
342 type(
matrix_t),
intent(inout) :: avg_planes
343 integer :: n, ierr, j, i
344 real(kind=
rp) :: coord
346 call avg_planes%free()
347 call avg_planes%init(this%n_gll_lvls,
size(vector_ptr)+1)
353 if (this%dir .eq. 1) coord = this%dof%x(i,1,1,1)
354 if (this%dir .eq. 2) coord = this%dof%y(i,1,1,1)
355 if (this%dir .eq. 3) coord = this%dof%z(i,1,1,1)
356 avg_planes%x(this%pt_lvl(i,1,1,1),1) = &
357 avg_planes%x(this%pt_lvl(i,1,1,1),1) + coord*this%coef%B(i,1,1,1) &
358 /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
360 do j = 2,
size(vector_ptr)+1
362 avg_planes%x(this%pt_lvl(i,1,1,1),j) = &
363 avg_planes%x(this%pt_lvl(i,1,1,1),j) + &
364 vector_ptr(j-1)%ptr%x(i)*this%coef%B(i,1,1,1) &
365 /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
368 call mpi_allreduce(mpi_in_place,avg_planes%x, (
size(vector_ptr)+1)*this%n_gll_lvls, &