FE-Project
Loading...
Searching...
No Matches
scale_element_operation_general.F90
Go to the documentation of this file.
1
2!-------------------------------------------------------------------------------
11#include "scaleFElib.h"
13
14 !-----------------------------------------------------------------------------
15 !
16 !++ used modules
17 !
18 use scale_precision
19
20 use scale_element_base, only: &
23
24 use scale_sparsemat, only: &
27
29
30 !-----------------------------------------------------------------------------
31 implicit none
32 private
33
34 !-----------------------------------------------------------------------------
35 !
36 !++ Public type & procedure
37 !
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 :: vfilterpm1 => element_operation_general_vfilterpm1
58 !-
59 procedure, public :: setup_modalfilter => element_operation_general_setup_modalfilter
60 procedure, public :: setup_modalfilter_tracer => element_operation_general_setup_modalfilter_tracer
61 procedure, public :: modalfilter_tracer => element_operation_general_modalfilter_tracer
62 procedure, public :: modalfilter_var5 => element_operation_general_modalfilter_var5
64
66 module procedure element_operation_general_generate_vpordm1
67 end interface
69
70contains
71
74!OCL SERIAL
75 subroutine element_operation_general_init( this, elem3D, &
76 Dx, Dy, Dz, Lift )
77 implicit none
78 class(ElementOperationGeneral), intent(inout) :: this
79 class(ElementBase3D), intent(in), target :: elem3D
80 type(SparseMat), intent(in), target :: Dx
81 type(SparseMat), intent(in), target :: Dy
82 type(SparseMat), intent(in), target :: Dz
83 type(SparseMat), intent(in), target :: Lift
84 !----------------------------------------------------------
85
86 this%elem3D => elem3d
87 this%Dx_sm => dx
88 this%Dy_sm => dy
89 this%Dz_sm => dz
90 this%Lift_sm => lift
91
92 !--
93 allocate( this%IntrpMat_VPOrdM1(elem3d%Np,elem3d%Np) )
94 call element_operation_general_generate_vpordm1( this%IntrpMat_VPOrdM1, &
95 elem3d )
96
97 return
98 end subroutine element_operation_general_init
99
102!OCL SERIAL
103 subroutine element_operation_general_generate_vpordm1( IntrpMat_VPOrdM1, &
104 elem3D )
105 implicit none
106 class(elementbase3d), intent(in), target :: elem3D
107 real(RP), intent(out) :: IntrpMat_VPOrdM1(elem3D%Np,elem3D%Np)
108
109 integer :: p1, p2, p_
110 real(RP) :: invV_VPOrdM1(elem3D%Np,elem3D%Np)
111 !----------------------------------------------------------
112
113 invv_vpordm1(:,:) = elem3d%invV(:,:)
114 do p2=1, elem3d%Nnode_h1D
115 do p1=1, elem3d%Nnode_h1D
116 p_ = p1 + (p2-1)*elem3d%Nnode_h1D + (elem3d%Nnode_v-1)*elem3d%Nnode_h1D**2
117 invv_vpordm1(p_,:) = 0.0_rp
118 end do
119 end do
120 intrpmat_vpordm1(:,:) = matmul(elem3d%V, invv_vpordm1)
121
122 return
123 end subroutine element_operation_general_generate_vpordm1
124
127!OCL SERIAL
128 subroutine element_operation_general_setup_modalfilter( this, &
129 MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
130 MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
131
132 implicit none
133 class(elementoperationgeneral), intent(inout) :: this
134 real(RP), intent(in) :: MF_ETAC_h
135 real(RP), intent(in) :: MF_ALPHA_h
136 integer, intent(in) :: MF_ORDER_h
137 real(RP), intent(in) :: MF_ETAC_v
138 real(RP), intent(in) :: MF_ALPHA_v
139 integer, intent(in) :: MF_ORDER_v
140 !--------------------------------------------------------
141
142 call setup_modalfilter( this%MFilter, &
143 mf_etac_h, mf_alpha_h, mf_order_h, &
144 mf_etac_v, mf_alpha_v, mf_order_v, &
145 this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
146
147 return
148 end subroutine element_operation_general_setup_modalfilter
149
152!OCL SERIAL
153 subroutine element_operation_general_setup_modalfilter_tracer( this, &
154 MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
155 MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
156 implicit none
157 class(elementoperationgeneral), intent(inout) :: this
158 real(RP), intent(in) :: MF_ETAC_h
159 real(RP), intent(in) :: MF_ALPHA_h
160 integer, intent(in) :: MF_ORDER_h
161 real(RP), intent(in) :: MF_ETAC_v
162 real(RP), intent(in) :: MF_ALPHA_v
163 integer, intent(in) :: MF_ORDER_v
164 !--------------------------------------------------------
165
166 call setup_modalfilter( this%MFilter_tracer, &
167 mf_etac_h, mf_alpha_h, mf_order_h, &
168 mf_etac_v, mf_alpha_v, mf_order_v, &
169 this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v )
170
171 return
172 end subroutine element_operation_general_setup_modalfilter_tracer
173
174
177 !OCL SERIAL
178 subroutine element_operation_general_final( this )
179 implicit none
180 class(elementoperationgeneral), intent(inout) :: this
181 !----------------------------------------------------------
182
183 nullify( this%elem3D )
184 nullify( this%Dx_sm, this%Dy_sm, this%Dz_sm, this%Lift_sm )
185
186 deallocate( this%IntrpMat_VPOrdM1 )
187
188 return
189 end subroutine element_operation_general_final
190
193!OCL SERIAL
194 subroutine element_operation_general_dx( this, vec_in, vec_out )
195 implicit none
196 class(elementoperationgeneral), intent(in) :: this
197 real(RP), intent(in) :: vec_in(this%elem3D%Np)
198 real(RP), intent(out) :: vec_out(this%elem3D%Np)
199 !----------------------------------------------------------
200 call sparsemat_matmul( this%Dx_sm, vec_in, vec_out )
201 return
202 end subroutine element_operation_general_dx
203
206!OCL SERIAL
207 subroutine element_operation_general_dy( this, vec_in, vec_out )
208 implicit none
209 class(elementoperationgeneral), intent(in) :: this
210 real(RP), intent(in) :: vec_in(this%elem3D%Np)
211 real(RP), intent(out) :: vec_out(this%elem3D%Np)
212 !----------------------------------------------------------
213 call sparsemat_matmul( this%Dy_sm, vec_in, vec_out )
214 return
215 end subroutine element_operation_general_dy
216
219!OCL SERIAL
220 subroutine element_operation_general_dz( this, vec_in, vec_out )
221 implicit none
222 class(elementoperationgeneral), intent(in) :: this
223 real(RP), intent(in) :: vec_in(this%elem3D%Np)
224 real(RP), intent(out) :: vec_out(this%elem3D%Np)
225 !----------------------------------------------------------
226 call sparsemat_matmul( this%Dz_sm, vec_in, vec_out )
227 return
228 end subroutine element_operation_general_dz
229
232!OCL SERIAL
233 subroutine element_operation_general_lift( this, vec_in, vec_out )
234 implicit none
235 class(elementoperationgeneral), intent(in) :: this
236 real(RP), intent(in) :: vec_in(this%elem3D%NfpTot)
237 real(RP), intent(out) :: vec_out(this%elem3D%Np)
238 !----------------------------------------------------------
239 call sparsemat_matmul( this%Lift_sm, vec_in, vec_out )
240 return
241 end subroutine element_operation_general_lift
242
245!OCL SERIAL
246 subroutine element_operation_general_dxdydzlift( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
247 implicit none
248 class(elementoperationgeneral), intent(in) :: this
249 real(RP), intent(in) :: vec_in(this%elem3D%Np)
250 real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
251 real(RP), intent(out) :: vec_out_dx(this%elem3D%Np)
252 real(RP), intent(out) :: vec_out_dy(this%elem3D%Np)
253 real(RP), intent(out) :: vec_out_dz(this%elem3D%Np)
254 real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
255 !----------------------------------------------------------
256
257 call sparsemat_matmul( this%Dx_sm, vec_in, vec_out_dx )
258 call sparsemat_matmul( this%Dy_sm, vec_in, vec_out_dy )
259 call sparsemat_matmul( this%Dz_sm, vec_in, vec_out_dz )
260 call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out_lift )
261 return
262 end subroutine element_operation_general_dxdydzlift
263
266!OCL SERIAL
267 subroutine element_operation_general_div( this, vec_in_x, vec_in_y, vec_in_z, vec_in_lift, &
268 vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
269 implicit none
270 class(elementoperationgeneral), intent(in) :: this
271 real(RP), intent(in) :: vec_in_x(this%elem3D%Np)
272 real(RP), intent(in) :: vec_in_y(this%elem3D%Np)
273 real(RP), intent(in) :: vec_in_z(this%elem3D%Np)
274 real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
275 real(RP), intent(out) :: vec_out_dx(this%elem3D%Np)
276 real(RP), intent(out) :: vec_out_dy(this%elem3D%Np)
277 real(RP), intent(out) :: vec_out_dz(this%elem3D%Np)
278 real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
279 !---------------------------------------------------------------
280
281 call sparsemat_matmul( this%Dx_sm, vec_in_x, vec_out_dx )
282 call sparsemat_matmul( this%Dy_sm, vec_in_y, vec_out_dy )
283 call sparsemat_matmul( this%Dz_sm, vec_in_z, vec_out_dz )
284 call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out_lift )
285 return
286 end subroutine element_operation_general_div
287
288
291!OCL SERIAL
292 subroutine element_operation_general_div_var5( this, vec_in, vec_in_lift, &
293 vec_out_d )
294 implicit none
295 class(elementoperationgeneral), intent(in) :: this
296 real(RP), intent(in) :: vec_in(this%elem3D%Np,3,5)
297 real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot,5)
298 real(RP), intent(out) :: vec_out_d(this%elem3D%Np,4,5)
299
300 integer :: iv
301 !---------------------------------------------------------------
302
303 do iv=1, 5
304 call sparsemat_matmul( this%Dx_sm, vec_in(:,1,iv), vec_out_d(:,1,iv) )
305 call sparsemat_matmul( this%Dy_sm, vec_in(:,2,iv), vec_out_d(:,2,iv) )
306 call sparsemat_matmul( this%Dz_sm, vec_in(:,3,iv), vec_out_d(:,3,iv) )
307 end do
308 do iv=1, 5
309 call sparsemat_matmul( this%Lift_sm, vec_in_lift(:,iv), vec_out_d(:,4,iv) )
310 end do
311 return
312 end subroutine element_operation_general_div_var5
313
314!OCL SERIAL
315 subroutine element_operation_general_vfilterpm1( this, vec_in, vec_out )
316 implicit none
317 class(elementoperationgeneral), intent(in) :: this
318 real(RP), intent(in) :: vec_in(this%elem3D%Np)
319 real(RP), intent(out) :: vec_out(this%elem3D%Np)
320 !---------------------------------------------------------------
321
322 call matmul_( this%IntrpMat_VPOrdM1, vec_in, this%elem3D%Np, &
323 vec_out )
324 return
325 end subroutine element_operation_general_vfilterpm1
326!--
327!OCL SERIAL
328 subroutine matmul_( IntrpMat_VPOrdM1, vec_in_, Np, vec_out_ )
329 implicit none
330 integer, intent(in) :: Np
331 real(RP), intent(in) :: IntrpMat_VPOrdM1(Np,Np)
332 real(RP), intent(in) :: vec_in_(Np)
333 real(RP), intent(out) :: vec_out_(Np)
334 !-------------------------------------------
335 vec_out_(:) = matmul( intrpmat_vpordm1(:,:), vec_in_(:) )
336 return
337 end subroutine matmul_
338
339!OCL SERIAL
340 subroutine element_operation_general_modalfilter_tracer( this, vec_in, vec_work, vec_out )
341 implicit none
342 class(elementoperationgeneral), intent(in) :: this
343 real(RP), intent(in) :: vec_in(this%elem3D%Np)
344 real(RP), intent(out) :: vec_work(this%elem3D%Np)
345 real(RP), intent(out) :: vec_out(this%elem3D%Np)
346
347 integer :: ii, kk, Np
348 real(RP) :: Mik
349 !---------------------------------------------
350
351 np = this%elem3D%Np
352 vec_out(:) = 0.0_rp
353
354 do ii=1, np
355 do kk=1, np
356 mik = this%MFilter_tracer%FilterMat(ii,kk)
357 vec_out(ii) = vec_out(ii) + mik * vec_in(kk)
358 end do
359 end do
360
361 return
362 end subroutine element_operation_general_modalfilter_tracer
363
364!OCL SERIAL
365 subroutine element_operation_general_modalfilter_var5( this, vec_in, vec_work, vec_out )
366 implicit none
367 class(elementoperationgeneral), intent(in) :: this
368 real(RP), intent(in) :: vec_in(this%elem3D%Np,5)
369 real(RP), intent(out) :: vec_work(this%elem3D%Np)
370 real(RP), intent(out) :: vec_out(this%elem3D%Np,5)
371
372 integer :: ii, kk, Np
373 real(RP) :: Mik
374 !---------------------------------------------
375
376 np = this%elem3D%Np
377 vec_out(:,:) = 0.0_rp
378
379 do ii=1, np
380 do kk=1, np
381 mik = this%MFilter%FilterMat(ii,kk)
382
383 vec_out(ii,1) = vec_out(ii,1) + mik * vec_in(kk,1)
384 vec_out(ii,2) = vec_out(ii,2) + mik * vec_in(kk,2)
385 vec_out(ii,3) = vec_out(ii,3) + mik * vec_in(kk,3)
386 vec_out(ii,4) = vec_out(ii,4) + mik * vec_in(kk,4)
387 vec_out(ii,5) = vec_out(ii,5) + mik * vec_in(kk,5)
388 end do
389 end do
390
391 return
392 end subroutine element_operation_general_modalfilter_var5
393
394!- private -
395
396!OCL SERIAL
397 subroutine setup_modalfilter( MFilter, &
398 MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
399 MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v, &
400 PolyOrder_h, PolyOrder_v )
401
403 implicit none
404
405 class(modalfilter), intent(inout) :: MFilter
406 real(RP), intent(in) :: MF_ETAC_h
407 real(RP), intent(in) :: MF_ALPHA_h
408 integer, intent(in) :: MF_ORDER_h
409 real(RP), intent(in) :: MF_ETAC_v
410 real(RP), intent(in) :: MF_ALPHA_v
411 integer, intent(in) :: MF_ORDER_v
412 integer, intent(in) :: PolyOrder_h
413 integer, intent(in) :: PolyOrder_v
414
415 type(hexahedralelement) :: elem3D
416 !--------------------------------------------------------
417
418 call elem3d%Init( polyorder_h, polyorder_v, .false. )
419
420 call mfilter%Init( &
421 elem3d, & ! (in)
422 mf_etac_h, mf_alpha_h, mf_order_h, & ! (in)
423 mf_etac_v, mf_alpha_v, mf_order_v ) ! (in)
424
425 call elem3d%Final()
426 return
427 end subroutine setup_modalfilter
428
429end module scale_element_operation_general
module FElib / Element / Base
subroutine, public elementbase3d_init(elem, lumpedmat_flag)
subroutine, public elementbase3d_final(elem)
module FElib / Element / hexahedron
module FElib / element/ ModalFilter
module FElib / Element / Operation / Base
module FElib / Element / Operation with arbitary elements
module common / sparsemat