Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
map_1d.f90
Go to the documentation of this file.
1
3module map_1d
5 use num_types, only : rp
6 use space, only : space_t
7 use dofmap, only : dofmap_t
8 use gather_scatter, only : gs_op_add
9 use mesh, only : mesh_t
13 use coefs, only : coef_t
14 use field_list, only : field_list_t
15 use matrix, only : matrix_t
16 use vector, only : vector_ptr_t
17 use utils, only : neko_error, neko_warning
18 use math, only : glmax, glmin, glimax, relcmp, cmult, add2s1, col2
19 use mpi_f08, only : mpi_allreduce, mpi_sum, mpi_barrier, mpi_in_place
20 use, intrinsic :: iso_c_binding
21 implicit none
22 private
27 type, public :: map_1d_t
30 integer, allocatable :: dir_el(:)
32 integer, allocatable :: el_lvl(:)
34 integer, allocatable :: pt_lvl(:,:,:,:)
36 integer :: n_el_lvls
38 integer :: n_gll_lvls
41 type(dofmap_t), pointer :: dof => null()
43 type(coef_t), pointer :: coef => null()
45 type(mesh_t), pointer :: msh => null()
47 integer :: dir
49 real(kind=rp) :: tol = 1e-7
51 real(kind=rp), allocatable :: volume_per_gll_lvl(:)
52 contains
54 procedure, pass(this) :: init_int => map_1d_init
55
56 procedure, pass(this) :: init_char => map_1d_init_char
57 generic :: init => init_int, init_char
59 procedure, pass(this) :: free => map_1d_free
61 procedure, pass(this) :: average_planes_fld_lst => map_1d_average_field_list
62
63 procedure, pass(this) :: average_planes_vec_ptr => map_1d_average_vector_ptr
64 generic :: average_planes => average_planes_fld_lst, average_planes_vec_ptr
65 end type map_1d_t
66
67
68contains
69
70 subroutine map_1d_init(this, coef, dir, tol)
71 class(map_1d_t) :: this
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
81
82 call this%free()
83
84 if (neko_bcknd_device .eq. 1) then
85 if (pe_rank .eq. 0) then
86 call neko_warning('map_1d does not copy indices to device,'// &
87 ' but ok if used on cpu and for io')
88 end if
89 end if
90
91 this%dir = dir
92 this%dof => coef%dof
93 this%coef => coef
94 this%msh => coef%msh
95 nelv = this%msh%nelv
96 lx = this%dof%Xh%lx
97 n = this%dof%size()
98
99 if (dir .eq. 1) then
100 line => this%dof%x
101 else if (dir .eq. 2) then
102 line => this%dof%y
103 else if(dir .eq. 3) then
104 line => this%dof%z
105 else
106 call neko_error('Invalid dir for geopmetric comm')
107 end if
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))
113 call mpi_barrier(neko_comm)
114 if (neko_bcknd_device .eq. 1) then
115
116 call device_map(min_vals,min_vals_d,n)
117 end if
118 call mpi_barrier(neko_comm)
119
120 do i = 1, nelv
121 !store which direction r,s,t corresponds to speciifed direction, x,y,z
122 !we assume elements are stacked on each other...
123 ! Check which one of the normalized vectors are closest to dir
124 ! If we want to incorporate other directions, we should look here
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,:))
134 ! Checks which directions in rst the xyz corresponds to
135 ! 1 corresponds to r, 2 to s, 3 to t and are stored in dir_el
136 this%dir_el(i) = maxloc(el_dim(:,this%dir),dim=1)
137 end do
138 glb_min = glmin(line,n)
139 glb_max = glmax(line,n)
140
141 i = 1
142 this%el_lvl = -1
143 ! Check what the mimum value in each element and put in min_vals
144 do e = 1, nelv
145 el_min = minval(line(:,:,:,e))
146 min_vals(:,:,:,e) = el_min
147 ! Check if this element is on the bottom, in this case assign el_lvl = i = 1
148 if (relcmp(el_min, glb_min, this%tol)) then
149 if(this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
150 end if
151 end do
152 ! While loop where at each iteation the global maximum value
153 ! propagates down one level.
154 ! When the minumum value has propagated to the highest level this stops.
155 ! Only works when the bottom plate of the domain is flat.
156 do while (.not. relcmp(glmax(min_vals,n), glb_min, this%tol))
157 i = i + 1
158 do e = 1, nelv
159 !Sets the value at the bottom of each element to glb_max
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
163 else
164 min_vals(1,:,:,e) = glb_max
165 end if
166 end if
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
170 else
171 min_vals(:,1,:,e) = glb_max
172 end if
173 end if
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
177 else
178 min_vals(:,:,1,e) = glb_max
179 end if
180 end if
181 end do
182 !Make sketchy min as GS_OP_MIN is not supported with device mpi
183 min_temp = min_vals
184 if (neko_bcknd_device .eq. 1) &
185 call device_memcpy(min_vals, min_vals_d, n,&
186 host_to_device, sync=.false.)
187 !Propagates the minumum value along the element boundary.
188 call coef%gs_h%op(min_vals, n, gs_op_add)
189 if (neko_bcknd_device .eq. 1) &
190 call device_memcpy(min_vals, min_vals_d, n,&
191 device_to_host, sync=.true.)
192 !Obtain average along boundary
193
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)
197
198
199 !Checks the new minimum value on each element
200 !Assign this value to all points in this element in min_val
201 !If the element has not already been assinged a level,
202 !and it has obtained the minval, set el_lvl = i
203 do e = 1, nelv
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
208 end if
209 end do
210 end do
211 this%n_el_lvls = glimax(this%el_lvl,nelv)
212 this%n_gll_lvls = this%n_el_lvls*lx
213 !Numbers the points in each element based on the element level
214 !and its orientation
215 do e = 1, nelv
216 do i = 1, 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
221 else
222 this%pt_lvl(i,:,:,e) = lvl
223 end if
224 end if
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
228 else
229 this%pt_lvl(:,i,:,e) = lvl
230 end if
231 end if
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
235 else
236 this%pt_lvl(:,:,i,e) = lvl
237 end if
238 end if
239 end do
240 end do
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
246 do i = 1, n
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)
249 end do
250 call mpi_allreduce(mpi_in_place,this%volume_per_gll_lvl, this%n_gll_lvls, &
251 mpi_real_precision, mpi_sum, neko_comm, ierr)
252
253 end subroutine map_1d_init
254
255 subroutine map_1d_init_char(this, coef, dir, tol)
256 class(map_1d_t) :: this
257 type(coef_t), intent(inout), target :: coef
258 character(len=*), intent(in) :: dir
259 real(kind=rp), intent(in) :: tol
260 integer :: idir
261
262 if (trim(dir) .eq. 'yz' .or. trim(dir) .eq. 'zy') then
263 idir = 1
264 else if (trim(dir) .eq. 'xz' .or. trim(dir) .eq. 'zx') then
265 idir = 2
266 else if (trim(dir) .eq. 'xy' .or. trim(dir) .eq. 'yx') then
267 idir = 3
268 else
269 call neko_error('homogenous direction not supported')
270 end if
271
272 call this%init(coef,idir,tol)
273
274 end subroutine map_1d_init_char
275
276 subroutine map_1d_free(this)
277 class(map_1d_t) :: this
278
279 if(allocated(this%dir_el)) deallocate(this%dir_el)
280 if(allocated(this%el_lvl)) deallocate(this%el_lvl)
281 if(allocated(this%pt_lvl)) deallocate(this%pt_lvl)
282 if(associated(this%dof)) nullify(this%dof)
283 if(associated(this%msh)) nullify(this%msh)
284 if(associated(this%coef)) nullify(this%coef)
285 if(allocated(this%volume_per_gll_lvl)) deallocate(this%volume_per_gll_lvl)
286 this%dir = 0
287 this%n_el_lvls = 0
288 this%n_gll_lvls = 0
289
290 end subroutine map_1d_free
291
292
298 subroutine map_1d_average_field_list(this, avg_planes, field_list)
299 class(map_1d_t), intent(inout) :: this
300 type(field_list_t), intent(inout) :: field_list
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)
306 avg_planes = 0.0_rp
307 !ugly way of getting coordinates, computes average
308 n = this%dof%size()
309 do i = 1, n
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))
316 end do
317 do j = 2, field_list%size() + 1
318 do i = 1, n
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))
323 end do
324 end do
325 if (pe_size .gt. 1) then
326 call mpi_allreduce(mpi_in_place, avg_planes%x, (field_list%size()+1)*this%n_gll_lvls, &
327 mpi_real_precision, mpi_sum, neko_comm, ierr)
328 end if
329
330
331 end subroutine map_1d_average_field_list
332
338 subroutine map_1d_average_vector_ptr(this, avg_planes, vector_ptr)
339 class(map_1d_t), intent(inout) :: this
340 !Observe is an array...
341 type(vector_ptr_t), intent(inout) :: vector_ptr(:)
342 type(matrix_t), intent(inout) :: avg_planes
343 integer :: n, ierr, j, i
344 real(kind=rp) :: coord
345
346 call avg_planes%free()
347 call avg_planes%init(this%n_gll_lvls,size(vector_ptr)+1)
348 !ugly way of getting coordinates, computes average
349 avg_planes = 0.0_rp
350
351 n = this%dof%size()
352 do i = 1, n
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))
359 end do
360 do j = 2, size(vector_ptr)+1
361 do i = 1, n
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))
366 end do
367 end do
368 call mpi_allreduce(mpi_in_place,avg_planes%x, (size(vector_ptr)+1)*this%n_gll_lvls, &
369 mpi_real_precision, mpi_sum, neko_comm, ierr)
370
371
372 end subroutine map_1d_average_vector_ptr
373
374end module map_1d
Map a Fortran array to a device (allocate and associate)
Definition device.F90:72
Copy data between host and device (or device and device)
Definition device.F90:66
Coefficients.
Definition coef.f90:34
Definition comm.F90:1
type(mpi_datatype), public mpi_real_precision
MPI type for working precision of REAL types.
Definition comm.F90:50
integer, public pe_size
MPI size of communicator.
Definition comm.F90:58
integer, public pe_rank
MPI rank.
Definition comm.F90:55
type(mpi_comm), public neko_comm
MPI communicator.
Definition comm.F90:42
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:214
integer, parameter, public device_to_host
Definition device.F90:47
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Gather-scatter.
Creates a 1d GLL point map along a specified direction based on the connectivity in the mesh.
Definition map_1d.f90:3
subroutine map_1d_init(this, coef, dir, tol)
Definition map_1d.f90:71
subroutine map_1d_free(this)
Definition map_1d.f90:277
subroutine map_1d_average_field_list(this, avg_planes, field_list)
Computes average if field list in two directions and outputs matrix with averaged values avg_planes c...
Definition map_1d.f90:299
subroutine map_1d_init_char(this, coef, dir, tol)
Definition map_1d.f90:256
subroutine map_1d_average_vector_ptr(this, avg_planes, vector_ptr)
Computes average if vector_pt in two directions and outputs matrix with averaged values avg_planes co...
Definition map_1d.f90:339
Definition math.f90:60
subroutine, public cmult(a, c, n)
Multiplication by constant c .
Definition math.f90:417
subroutine, public add2s1(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
Definition math.f90:803
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:860
real(kind=rp) function, public glmax(a, n)
Max of a vector of length n.
Definition math.f90:522
integer function, public glimax(a, n)
Max of an integer vector of length n.
Definition math.f90:537
real(kind=rp) function, public glmin(a, n)
Min of a vector of length n.
Definition math.f90:552
Defines a matrix.
Definition matrix.f90:34
Defines a mesh.
Definition mesh.f90:34
Build configurations.
integer, parameter neko_bcknd_device
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:284
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:27
The function space for the SEM solution fields.
Definition space.f90:62