FE-Project
Loading...
Searching...
No Matches
scale_element_operation_general.F90
Go to the documentation of this file.
1
2!-------------------------------------------------------------------------------
3!> module FElib / Element / Operation with arbitary elements
4!!
5!! @par Description
6!! A module for providing mathematical operations with arbitary elements using a module for SpMV
7!!
8!! @author Yuta Kawai, Xuanzhengbo Ren, and Team SCALE
9!!
10!<
11#include "scaleFElib.h"
12module scale_element_operation_general
13
14 !-----------------------------------------------------------------------------
15 !
16 !++ used modules
17 !
18 use scale_precision
19
20 use scale_element_base, only: &
21 elementbase3d, &
22 elementbase3d_init, elementbase3d_final
23
24 use scale_sparsemat, only: &
25 sparsemat, sparsemat_matmul
26 use scale_element_modalfilter, only: modalfilter
27
28 use scale_element_operation_base, only: elementoperationbase3d
29
30 !-----------------------------------------------------------------------------
31 implicit none
32 private
33
34 !-----------------------------------------------------------------------------
35 !
36 !++ Public type & procedure
37 !
38 type, public, extends(elementoperationbase3d) :: elementoperationgeneral
39 type(sparsemat), pointer :: dx_sm
40 type(sparsemat), pointer :: dy_sm
41 type(sparsemat), pointer :: dz_sm
42 type(sparsemat), pointer :: lift_sm
43 real(rp), allocatable :: intrpmat_vpordm1(:,:)
44
45 type(modalfilter) :: mfilter
46 type(modalfilter) :: mfilter_tracer
47 contains
48 procedure, public :: init => element_operation_general_init
49 procedure, public :: final => element_operation_general_final
50 procedure, public :: dx => element_operation_general_dx
51 procedure, public :: dy => element_operation_general_dy
52 procedure, public :: dz => element_operation_general_dz
53 procedure, public :: lift => element_operation_general_lift
54 procedure, public :: dxdydzlift => element_operation_general_dxdydzlift
55 procedure, public :: div => element_operation_general_div
56 procedure, public :: div_var5 => element_operation_general_div_var5
57 procedure, public :: div_var5_2 => element_operation_general_div_var5_2
58 procedure, public :: lift_var5 => element_operation_general_lift_var5
59 procedure, public :: vfilterpm1 => element_operation_general_vfilterpm1
60 !-
61 procedure, public :: setup_modalfilter => element_operation_general_setup_modalfilter
62 procedure, public :: setup_modalfilter_tracer => element_operation_general_setup_modalfilter_tracer
63 procedure, public :: modalfilter_tracer => element_operation_general_modalfilter_tracer
64 procedure, public :: modalfilter_var5 => element_operation_general_modalfilter_var5
65 end type elementoperationgeneral
66
67 interface elementoperationgeneral_generate_vpordm1
68 module procedure element_operation_general_generate_vpordm1
69 end interface
70 public :: elementoperationgeneral_generate_vpordm1
71
72contains
73
74 !> Initialization
75 !!
76!OCL SERIAL
77 subroutine element_operation_general_init( this, elem3D, &
78 Dx, Dy, Dz, Lift )
79 implicit none
80 class(ElementOperationGeneral), intent(inout) :: this
81 class(ElementBase3D), intent(in), target :: elem3D
82 type(SparseMat), intent(in), target :: Dx
83 type(SparseMat), intent(in), target :: Dy
84 type(SparseMat), intent(in), target :: Dz
85 type(SparseMat), intent(in), target :: Lift
86 !----------------------------------------------------------
87
88 this%elem3D => elem3d
89 this%Dx_sm => dx
90 this%Dy_sm => dy
91 this%Dz_sm => dz
92 this%Lift_sm => lift
93
94 !--
95 allocate( this%IntrpMat_VPOrdM1(elem3d%Np,elem3d%Np) )
96 call element_operation_general_generate_vpordm1( this%IntrpMat_VPOrdM1, &
97 elem3d )
98
99 return
100 end subroutine element_operation_general_init
101
102 !> Generate a vertical filter matrix to remove the highest mode
103 !!
104!OCL SERIAL
105 subroutine element_operation_general_generate_vpordm1( IntrpMat_VPOrdM1, &
106 elem3D )
107 implicit none
108 class(elementbase3d), intent(in), target :: elem3D
109 real(RP), intent(out) :: IntrpMat_VPOrdM1(elem3D%Np,elem3D%Np)
110
111 integer :: p1, p2, p_
112 real(RP) :: invV_VPOrdM1(elem3D%Np,elem3D%Np)
113 !----------------------------------------------------------
114
115 invv_vpordm1(:,:) = elem3d%invV(:,:)
116 do p2=1, elem3d%Nnode_h1D
117 do p1=1, elem3d%Nnode_h1D
118 p_ = p1 + (p2-1)*elem3d%Nnode_h1D + (elem3d%Nnode_v-1)*elem3d%Nnode_h1D**2
119 invv_vpordm1(p_,:) = 0.0_rp
120 end do
121 end do
122 intrpmat_vpordm1(:,:) = matmul(elem3d%V, invv_vpordm1)
123
124 return
125 end subroutine element_operation_general_generate_vpordm1
126
127 !> Setup modal filter
128 !!
129!OCL SERIAL
130 subroutine element_operation_general_setup_modalfilter( this, &
131 MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
132 MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
133
134 implicit none
135 class(elementoperationgeneral), intent(inout) :: this
136 real(RP), intent(in) :: MF_ETAC_h
137 real(RP), intent(in) :: MF_ALPHA_h
138 integer, intent(in) :: MF_ORDER_h
139 real(RP), intent(in) :: MF_ETAC_v
140 real(RP), intent(in) :: MF_ALPHA_v
141 integer, intent(in) :: MF_ORDER_v
142 !--------------------------------------------------------
143
144 call setup_modalfilter( this%MFilter, &
145 mf_etac_h, mf_alpha_h, mf_order_h, &
146 mf_etac_v, mf_alpha_v, mf_order_v, &
147 this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
148
149 return
150 end subroutine element_operation_general_setup_modalfilter
151
152 !> Setup modal filter for tracer
153 !!
154!OCL SERIAL
155 subroutine element_operation_general_setup_modalfilter_tracer( this, &
156 MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
157 MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
158 implicit none
159 class(elementoperationgeneral), intent(inout) :: this
160 real(RP), intent(in) :: MF_ETAC_h
161 real(RP), intent(in) :: MF_ALPHA_h
162 integer, intent(in) :: MF_ORDER_h
163 real(RP), intent(in) :: MF_ETAC_v
164 real(RP), intent(in) :: MF_ALPHA_v
165 integer, intent(in) :: MF_ORDER_v
166 !--------------------------------------------------------
167
168 call setup_modalfilter( this%MFilter_tracer, &
169 mf_etac_h, mf_alpha_h, mf_order_h, &
170 mf_etac_v, mf_alpha_v, mf_order_v, &
171 this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
172
173 return
174 end subroutine element_operation_general_setup_modalfilter_tracer
175
176
177 !> Finalization
178 !!
179 !OCL SERIAL
180 subroutine element_operation_general_final( this )
181 implicit none
182 class(elementoperationgeneral), intent(inout) :: this
183 !----------------------------------------------------------
184
185 nullify( this%elem3D )
186 nullify( this%Dx_sm, this%Dy_sm, this%Dz_sm, this%Lift_sm )
187
188 deallocate( this%IntrpMat_VPOrdM1 )
189
190 return
191 end subroutine element_operation_general_final
192
193!> Calculate the differential in x-direction
194!!
195!OCL SERIAL
196 subroutine element_operation_general_dx( this, vec_in, vec_out )
197 implicit none
198 class(elementoperationgeneral), intent(in) :: this
199 real(RP), intent(in) :: vec_in(this%elem3D%Np)
200 real(RP), intent(out) :: vec_out(this%elem3D%Np)
201 !----------------------------------------------------------
202 call sparsemat_matmul( this%Dx_sm, vec_in, vec_out )
203 return
204 end subroutine element_operation_general_dx
205
206!> Calculate the differential in y-direction
207!!
208!OCL SERIAL
209 subroutine element_operation_general_dy( this, vec_in, vec_out )
210 implicit none
211 class(elementoperationgeneral), intent(in) :: this
212 real(RP), intent(in) :: vec_in(this%elem3D%Np)
213 real(RP), intent(out) :: vec_out(this%elem3D%Np)
214 !----------------------------------------------------------
215 call sparsemat_matmul( this%Dy_sm, vec_in, vec_out )
216 return
217 end subroutine element_operation_general_dy
218
219!> Calculate the differential in z-direction
220!!
221!OCL SERIAL
222 subroutine element_operation_general_dz( this, vec_in, vec_out )
223 implicit none
224 class(elementoperationgeneral), intent(in) :: this
225 real(RP), intent(in) :: vec_in(this%elem3D%Np)
226 real(RP), intent(out) :: vec_out(this%elem3D%Np)
227 !----------------------------------------------------------
228 call sparsemat_matmul( this%Dz_sm, vec_in, vec_out )
229 return
230 end subroutine element_operation_general_dz
231
232!> Calculate the differential in z-direction
233!!
234!OCL SERIAL
235 subroutine element_operation_general_lift( this, vec_in, vec_out )
236 implicit none
237 class(elementoperationgeneral), intent(in) :: this
238 real(RP), intent(in) :: vec_in(this%elem3D%NfpTot)
239 real(RP), intent(out) :: vec_out(this%elem3D%Np)
240 !----------------------------------------------------------
241 call sparsemat_matmul( this%Lift_sm, vec_in, vec_out )
242 return
243 end subroutine element_operation_general_lift
244
245!> Calculate the 3D gradient
246!!
247!OCL SERIAL
248 subroutine element_operation_general_dxdydzlift( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
249 implicit none
250 class(elementoperationgeneral), intent(in) :: this
251 real(RP), intent(in) :: vec_in(this%elem3D%Np)
252 real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
253 real(RP), intent(out) :: vec_out_dx(this%elem3D%Np)
254 real(RP), intent(out) :: vec_out_dy(this%elem3D%Np)
255 real(RP), intent(out) :: vec_out_dz(this%elem3D%Np)
256 real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
257 !----------------------------------------------------------
258
259 call sparsemat_matmul( this%Dx_sm, vec_in, vec_out_dx )
260 call sparsemat_matmul( this%Dy_sm, vec_in, vec_out_dy )
261 call sparsemat_matmul( this%Dz_sm, vec_in, vec_out_dz )
262 call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out_lift )
263 return
264 end subroutine element_operation_general_dxdydzlift
265
266!> Calculate the 3D gradient
267!!
268!OCL SERIAL
269 subroutine element_operation_general_div( this, vec_in, vec_in_lift, &
270 vec_out )
271 implicit none
272 class(elementoperationgeneral), intent(in) :: this
273 real(RP), intent(in) :: vec_in(this%elem3D%Np,3)
274 real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
275 real(RP), intent(out) :: vec_out(this%elem3D%Np,4)
276 !---------------------------------------------------------------
277
278 call sparsemat_matmul( this%Dx_sm, vec_in(:,1), vec_out(:,1) )
279 call sparsemat_matmul( this%Dy_sm, vec_in(:,2), vec_out(:,2) )
280 call sparsemat_matmul( this%Dz_sm, vec_in(:,3), vec_out(:,3) )
281 call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out(:,4) )
282 return
283 end subroutine element_operation_general_div
284
285
286!> Calculate the 3D gradient
287!!
288!OCL SERIAL
289 subroutine element_operation_general_div_var5( this, vec_in, vec_in_lift, &
290 vec_out_d )
291 implicit none
292 class(elementoperationgeneral), intent(in) :: this
293 real(RP), intent(in) :: vec_in(this%elem3D%Np,3,5)
294 real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot,5)
295 real(RP), intent(out) :: vec_out_d(this%elem3D%Np,4,5)
296
297 integer :: iv
298 !---------------------------------------------------------------
299
300 do iv=1, 5
301 call sparsemat_matmul( this%Dx_sm, vec_in(:,1,iv), vec_out_d(:,1,iv) )
302 call sparsemat_matmul( this%Dy_sm, vec_in(:,2,iv), vec_out_d(:,2,iv) )
303 call sparsemat_matmul( this%Dz_sm, vec_in(:,3,iv), vec_out_d(:,3,iv) )
304 end do
305 do iv=1, 5
306 call sparsemat_matmul( this%Lift_sm, vec_in_lift(:,iv), vec_out_d(:,4,iv) )
307 end do
308 return
309 end subroutine element_operation_general_div_var5
310
311!> Calculate the 3D gradient
312!!
313!OCL SERIAL
314 subroutine element_operation_general_div_var5_2( this, vec_in, &
315 vec_out_d )
316 implicit none
317 class(elementoperationgeneral), intent(in) :: this
318 real(RP), intent(in) :: vec_in(this%elem3D%Np,3,5)
319 real(RP), intent(out) :: vec_out_d(this%elem3D%Np,3,5)
320
321 integer :: iv
322 !---------------------------------------------------------------
323
324 do iv=1, 5
325 call sparsemat_matmul( this%Dx_sm, vec_in(:,1,iv), vec_out_d(:,1,iv) )
326 call sparsemat_matmul( this%Dy_sm, vec_in(:,2,iv), vec_out_d(:,2,iv) )
327 call sparsemat_matmul( this%Dz_sm, vec_in(:,3,iv), vec_out_d(:,3,iv) )
328 end do
329 return
330 end subroutine element_operation_general_div_var5_2
331!> Calculate the differential in z-direction
332!!
333!OCL SERIAL
334 subroutine element_operation_general_lift_var5( this, vec_in, vec_out )
335 implicit none
336 class(elementoperationgeneral), intent(in) :: this
337 real(RP), intent(in) :: vec_in(this%elem3D%NfpTot,5)
338 real(RP), intent(out) :: vec_out(this%elem3D%Np,5)
339
340 integer :: iv
341 !----------------------------------------------------------
342 do iv=1,5
343 call sparsemat_matmul( this%Lift_sm, vec_in(:,iv), vec_out(:,iv) )
344 end do
345 return
346 end subroutine element_operation_general_lift_var5
347
348!OCL SERIAL
349 subroutine element_operation_general_vfilterpm1( this, vec_in, vec_out )
350 implicit none
351 class(elementoperationgeneral), intent(in) :: this
352 real(RP), intent(in) :: vec_in(this%elem3D%Np)
353 real(RP), intent(out) :: vec_out(this%elem3D%Np)
354 !---------------------------------------------------------------
355
356 call matmul_( this%IntrpMat_VPOrdM1, vec_in, this%elem3D%Np, &
357 vec_out )
358 return
359 end subroutine element_operation_general_vfilterpm1
360!--
361!OCL SERIAL
362 subroutine matmul_( IntrpMat_VPOrdM1, vec_in_, Np, vec_out_ )
363 implicit none
364 integer, intent(in) :: Np
365 real(RP), intent(in) :: IntrpMat_VPOrdM1(Np,Np)
366 real(RP), intent(in) :: vec_in_(Np)
367 real(RP), intent(out) :: vec_out_(Np)
368 !-------------------------------------------
369 vec_out_(:) = matmul( intrpmat_vpordm1(:,:), vec_in_(:) )
370 return
371 end subroutine matmul_
372
373!OCL SERIAL
374 subroutine element_operation_general_modalfilter_tracer( this, vec_in, vec_work, vec_out )
375 implicit none
376 class(elementoperationgeneral), intent(in) :: this
377 real(RP), intent(in) :: vec_in(this%elem3D%Np)
378 real(RP), intent(out) :: vec_work(this%elem3D%Np)
379 real(RP), intent(out) :: vec_out(this%elem3D%Np)
380
381 integer :: ii, kk, Np
382 real(RP) :: Mik
383 !---------------------------------------------
384
385 np = this%elem3D%Np
386 vec_out(:) = 0.0_rp
387
388 do ii=1, np
389 do kk=1, np
390 mik = this%MFilter_tracer%FilterMat(ii,kk)
391 vec_out(ii) = vec_out(ii) + mik * vec_in(kk)
392 end do
393 end do
394
395 return
396 end subroutine element_operation_general_modalfilter_tracer
397
398!OCL SERIAL
399 subroutine element_operation_general_modalfilter_var5( this, vec_in, vec_work, vec_out )
400 implicit none
401 class(elementoperationgeneral), intent(in) :: this
402 real(RP), intent(in) :: vec_in(this%elem3D%Np,5)
403 real(RP), intent(out) :: vec_work(this%elem3D%Np)
404 real(RP), intent(out) :: vec_out(this%elem3D%Np,5)
405
406 integer :: ii, kk, Np
407 real(RP) :: Mik
408 !---------------------------------------------
409
410 np = this%elem3D%Np
411 vec_out(:,:) = 0.0_rp
412
413 do ii=1, np
414 do kk=1, np
415 mik = this%MFilter%FilterMat(ii,kk)
416
417 vec_out(ii,1) = vec_out(ii,1) + mik * vec_in(kk,1)
418 vec_out(ii,2) = vec_out(ii,2) + mik * vec_in(kk,2)
419 vec_out(ii,3) = vec_out(ii,3) + mik * vec_in(kk,3)
420 vec_out(ii,4) = vec_out(ii,4) + mik * vec_in(kk,4)
421 vec_out(ii,5) = vec_out(ii,5) + mik * vec_in(kk,5)
422 end do
423 end do
424
425 return
426 end subroutine element_operation_general_modalfilter_var5
427
428!- private -
429
430!OCL SERIAL
431 subroutine setup_modalfilter( MFilter, &
432 MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
433 MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v, &
434 PolyOrder_h, PolyOrder_v )
435
436 use scale_element_hexahedral, only: hexahedralelement
437 implicit none
438
439 class(modalfilter), intent(inout) :: MFilter
440 real(RP), intent(in) :: MF_ETAC_h
441 real(RP), intent(in) :: MF_ALPHA_h
442 integer, intent(in) :: MF_ORDER_h
443 real(RP), intent(in) :: MF_ETAC_v
444 real(RP), intent(in) :: MF_ALPHA_v
445 integer, intent(in) :: MF_ORDER_v
446 integer, intent(in) :: PolyOrder_h
447 integer, intent(in) :: PolyOrder_v
448
449 type(hexahedralelement) :: elem3D
450 !--------------------------------------------------------
451
452 call elem3d%Init( polyorder_h, polyorder_v, .false. )
453
454 call mfilter%Init( &
455 elem3d, & ! (in)
456 mf_etac_h, mf_alpha_h, mf_order_h, & ! (in)
457 mf_etac_v, mf_alpha_v, mf_order_v ) ! (in)
458
459 call elem3d%Final()
460 return
461 end subroutine setup_modalfilter
462
463end module scale_element_operation_general
scale_element_operation_general::elementoperationgeneral_generate_vpordm1
Definition scale_element_operation_general.F90:67
scale_sparsemat::sparsemat_matmul
Definition scale_sparsemat.F90:57
scale_element_base
module FElib / Element / Base
Definition scale_element_base.F90:11
scale_element_base::elementbase3d_init
subroutine, public elementbase3d_init(elem, lumpedmat_flag)
Initialize an object to manage a 3D reference element.
Definition scale_element_base.F90:372
scale_element_base::elementbase3d_final
subroutine, public elementbase3d_final(elem)
Finalize an object to manage a 3D reference element.
Definition scale_element_base.F90:396
scale_element_hexahedral
module FElib / Element / hexahedron
Definition scale_element_hexahedral.F90:10
scale_element_modalfilter
module FElib / Element/ ModalFilter
Definition scale_element_modalfilter.F90:11
scale_element_operation_base
module FElib / Element / Operation / Base
Definition scale_element_operation_base.F90:12
scale_element_operation_general
module FElib / Element / Operation with arbitary elements
Definition scale_element_operation_general.F90:12
scale_sparsemat
Module common / sparsemat.
Definition scale_sparsemat.F90:11
scale_element_base::elementbase3d
Derived type representing a 3D reference element.
Definition scale_element_base.F90:107
scale_element_hexahedral::hexahedralelement
Derived type representing a hexahedral element.
Definition scale_element_hexahedral.F90:31
scale_element_modalfilter::modalfilter
Derived type representing a modal filter.
Definition scale_element_modalfilter.F90:42
scale_element_operation_base::elementoperationbase3d
Definition scale_element_operation_base.F90:31
scale_element_operation_general::elementoperationgeneral
Definition scale_element_operation_general.F90:38
scale_sparsemat::sparsemat
Derived type to manage a sparse matrix.
Definition scale_sparsemat.F90:33
Generated by doxygen 1.14.0