FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_trcadvect3d_heve.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module FElib / Fluid dyn solver / Atmosphere / Tracer advection
3!!
4!! @par Description
5!! HEVE DGM scheme for tracer advection.
6!!
7!! To preserve nonnegativity, a limiter proposed by Light and Durran (2016, MWR) is used:
8!! we apply FCT for the lowest mode and the truncation and mass aware rescaling (TMAR).
9!!
10!! @par Reference
11!! - Light and Durran 2016:
12!! Preserving Nonnegativity in Discontinuous Galerkin Approximations to Scalar Transport via Truncation and Mass Aware Rescaling (TMAR).
13!! Monthly Weather Review, 144(12), 4771–4786.
14!!
15!! @author Yuta Kawai, Team SCALE
16!<
17!-------------------------------------------------------------------------------
18#include "scaleFElib.h"
20 !-----------------------------------------------------------------------------
21 !
22 !++ Used modules
23 !
24 use scale_precision
25 use scale_io
26 use scale_prc
27 use scale_prof
28
29 use scale_sparsemat, only: &
31 use scale_element_base, only: &
41
42 !-----------------------------------------------------------------------------
43 implicit none
44 private
45 !-----------------------------------------------------------------------------
46 !
47 !++ Public procedures
48 !
56
57 !-----------------------------------------------------------------------------
58 !
59 !++ Public parameters & variables
60 !
61
62 !-----------------------------------------------------------------------------
63 !
64 !++ Private procedures & variables
65 !
66 !-------------------
67
68 private :: atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn
69 private :: atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc
70 private :: atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc
71
72contains
73
74!OCL SERIAL
75 subroutine atm_dyn_dgm_trcadvect3d_heve_init( mesh, FaceIntMat )
77 implicit none
78 class(meshbase3d), intent(in), target :: mesh
79 type(sparsemat), intent(inout) :: faceintmat
80
81 class(localmesh3d), pointer :: lcmesh
82 class(elementbase3d), pointer :: elem
83 real(rp), allocatable :: intweight_lgl1dpts_h(:)
84 real(rp), allocatable :: intweight_lgl1dpts_v(:)
85 real(rp), allocatable :: intweight_h(:)
86 real(rp), allocatable :: intweight_v(:)
87
88 integer :: f
89 integer :: i, j, k, l
90 integer :: is, ie
91
92 real(rp), allocatable :: intweight(:,:)
93 !-------------------------------------------------------------
94
95 lcmesh => mesh%lcmesh_list(1)
96 elem => lcmesh%refElem3D
97 allocate( intweight(elem%Nfaces,elem%NfpTot) )
98 intweight(:,:) = 0.0_rp
99
100 allocate( intweight_lgl1dpts_h(elem%Nnode_h1D) )
101 allocate( intweight_lgl1dpts_v(elem%Nnode_v) )
102 allocate( intweight_h(elem%Nnode_h1D*elem%Nnode_v) )
103 allocate( intweight_v(elem%Nnode_h1D**2) )
104
105 intweight_lgl1dpts_h(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder_h)
106 intweight_lgl1dpts_v(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder_v)
107
108 do f=1, elem%Nfaces_h
109 do k=1, elem%Nnode_v
110 do i=1, elem%Nnode_h1D
111 l = i + (k-1)*elem%Nnode_h1D
112 intweight_h(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_v(k)
113 end do
114 end do
115
116 is = (f-1)*elem%Nfp_h + 1
117 ie = is + elem%Nfp_h - 1
118 intweight(f,is:ie) = intweight_h(:)
119 end do
120
121 do f=1, elem%Nfaces_v
122 do j=1, elem%Nnode_h1D
123 do i=1, elem%Nnode_h1D
124 l = i + (j-1)*elem%Nnode_h1D
125 intweight_v(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_h(j)
126 end do
127 end do
128
129 is = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v + 1
130 ie = is + elem%Nfp_v - 1
131 intweight(elem%Nfaces_h+f,is:ie) = intweight_v(:)
132 end do
133
134 call faceintmat%Init( intweight )
135
136 return
138
139!OCL SERIAL
141 implicit none
142 !--------------------------------------------
143
144 return
146
147
148!OCL SERIAL
150 QTRC_dt, & ! (out)
151 qtrc_, momx_, momy_, momz_, & ! (in)
152 alphdens_m, alphdens_p, fct_coef, & ! (in)
153 rhoq_tp, & ! (in)
154 element3d_operation, faceintmat, & ! (in)
155 lmesh, elem, lmesh2d, elem2d ) ! (in)
156
157 class(localmesh3d), intent(in) :: lmesh
158 class(elementbase3d), intent(in) :: elem
159 class(localmesh2d), intent(in) :: lmesh2d
160 class(elementbase2d), intent(in) :: elem2d
161
162 real(rp), intent(out) :: qtrc_dt(elem%np,lmesh%nea)
163 real(rp), intent(in) :: qtrc_(elem%np,lmesh%nea)
164 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
165 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
166 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
167 real(rp), intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
168 real(rp), intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
169 real(rp), intent(in) :: fct_coef(elem%np,lmesh%nea)
170 real(rp), intent(in) :: rhoq_tp(elem%np,lmesh%nea)
171 class(elementoperationbase3d), intent(in) :: element3d_operation
172 type(sparsemat), intent(in) :: faceintmat
173
174 real(rp) :: flux(elem%np,3), dflux(elem%np,4)
175 real(rp) :: del_flux(elem%nfptot,lmesh%ne)
176 real(rp) :: momwt_
177 real(rp) :: gsqrt_
178 real(rp) :: rgsqrt(elem%np), rgsqrtv(elem%np)
179
180 integer :: ke, ke2d
181
182 real(rp) :: q0, q1, vol
183 integer :: p
184 !---------------------------------------------------------------------------
185
186 call prof_rapstart('cal_trcadv_tend_bndflux', 3)
187 call atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc( &
188 del_flux, & ! (out)
189 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef, & ! (in)
190 lmesh%Gsqrt, lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), & ! (in)
191 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
192 lmesh%J(:,:), lmesh%Fscale(:,:), & ! (in)
193 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, faceintmat, & ! (in)
194 lmesh, elem, lmesh2d, elem2d ) ! (in)
195 call prof_rapend('cal_trcadv_tend_bndflux', 3)
196
197 call prof_rapstart('cal_trcadv_tend_interior', 3)
198 !$omp parallel do private( ke, ke2D, &
199 !$omp momwt_, Flux, DFlux, Gsqrt_, RGsqrt, RGsqrtV )
200 do ke=lmesh%NeS, lmesh%NeE
201 ke2d = lmesh%EMap3Dto2D(ke)
202
203 do p=1, elem%Np
204 rgsqrt(p) = 1.0_rp / lmesh%Gsqrt(p,ke)
205 rgsqrtv(p) = lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d) * rgsqrt(p)
206 end do
207
208 do p=1, elem%Np
209 gsqrt_ = lmesh%Gsqrt(p,ke)
210 flux(p,1) = gsqrt_ * momx_(p,ke) * qtrc_(p,ke)
211 flux(p,2) = gsqrt_ * momy_(p,ke) * qtrc_(p,ke)
212 flux(p,3) = gsqrt_ * ( &
213 momz_(p,ke) * rgsqrtv(p) &
214 + lmesh%GI3(p,ke,1) * momx_(p,ke) &
215 + lmesh%GI3(p,ke,2) * momy_(p,ke) ) * qtrc_(p,ke)
216 end do
217 call element3d_operation%Div( &
218 flux, del_flux(:,ke), & ! (in)
219 dflux ) ! (out)
220
221 qtrc_dt(:,ke) = - ( &
222 lmesh%Escale(:,ke,1,1) * dflux(:,1) &
223 + lmesh%Escale(:,ke,2,2) * dflux(:,2) &
224 + lmesh%Escale(:,ke,3,3) * dflux(:,3) &
225 + dflux(:,4) ) / lmesh%Gsqrt(:,ke) &
226 + rhoq_tp(:,ke)
227 end do
228 call prof_rapend('cal_trcadv_tend_interior', 3)
229
230 return
232
233!OCL SERIAL
235 fct_coef, & ! (out)
236 qtrc_, momx_, momy_, momz_, rhoq_tp_, alphdens_m, alphdens_p, & ! (in)
237 dens_hyd, ddens_, ddens0_, rk_c_ssm1, dt, & ! (in)
238 faceintmat, lmesh, elem, lmesh2d, elem2d, & ! (in)
239 disable_limiter ) ! (in)
240
241 implicit none
242 class(localmesh3d), intent(in) :: lmesh
243 class(elementbase3d), intent(in) :: elem
244 class(localmesh2d), intent(in) :: lmesh2d
245 class(elementbase2d), intent(in) :: elem2d
246 real(rp), intent(out) :: fct_coef(elem%np,lmesh%nea)
247 real(rp), intent(in) :: qtrc_(elem%np,lmesh%nea)
248 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
249 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
250 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
251 real(rp), intent(in) :: rhoq_tp_(elem%np,lmesh%nea)
252 real(rp), intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
253 real(rp), intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
254 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
255 real(rp), intent(in) :: ddens_ (elem%np,lmesh%nea)
256 real(rp), intent(in) :: ddens0_(elem%np,lmesh%nea)
257 real(rp), intent(in) :: rk_c_ssm1
258 real(rp), intent(in) :: dt
259 type(sparsemat), intent(in) :: faceintmat
260 logical, intent(in), optional :: disable_limiter
261
262 real(rp) :: netoutwardflux(lmesh%ne)
263 real(rp) :: momwt_(elem%np)
264
265 integer :: ke
266 real(rp) :: q
267 real(rp) :: dens_ssm1(elem%np)
268 !---------------------------------------------------------------------------
269
270 if ( present(disable_limiter) ) then
271 if ( disable_limiter ) then
272 !$omp parallel do
273 do ke=lmesh%NeS, lmesh%NeE
274 fct_coef(:,ke) = 1.0_rp
275 end do
276 return
277 end if
278 end if
279
280 call prof_rapstart('cal_trcadv_fct_coef_bndflux', 3)
281 call atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc( &
282 netoutwardflux, & ! (out)
283 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, & ! (in)
284 lmesh%Gsqrt, lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), & ! (in)
285 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
286 lmesh%J(:,:), lmesh%Fscale(:,:), lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, & ! (in)
287 faceintmat, & ! (in)
288 lmesh, elem, lmesh2d, elem2d ) ! (in)
289 call prof_rapend('cal_trcadv_fct_coef_bndflux', 3)
290
291 call prof_rapstart('cal_trcadv_fct_coef', 3)
292
293 !$omp parallel do private( ke, Q, dens_ssm1 )
294 do ke=lmesh%NeS, lmesh%NeE
295
296 dens_ssm1(:) = dens_hyd(:,ke) &
297 + ( 1.0_rp - rk_c_ssm1 ) * ddens0_(:,ke) + rk_c_ssm1 * ddens_(:,ke)
298 q = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * ( dens_ssm1(:) * qtrc_(:,ke) / dt + rhoq_tp_(:,ke) ) )
299
300 fct_coef(:,ke) = max( 0.0_rp, min( 1.0_rp, q / ( netoutwardflux(ke) + 1.0e-10_rp ) ) )
301 end do
302
303 call prof_rapend('cal_trcadv_fct_coef', 3)
304
305 return
307
308
309!OCL SERIAL
310!> Second Step of limiter in which nonlinear truncation and mass aware rescaling (TMAR)
311 subroutine atm_dyn_dgm_trcadvect3d_tmar( QTRC_, & ! (inout)
312 dens_hyd, ddens_, lmesh, elem, lmesh2d, elem2d ) ! (in)
313
314 implicit none
315 class(localmesh3d), intent(in) :: lmesh
316 class(elementbase3d), intent(in) :: elem
317 class(localmesh2d), intent(in) :: lmesh2d
318 class(elementbase2d), intent(in) :: elem2d
319 real(rp), intent(inout) :: qtrc_(elem%np,lmesh%nea)
320 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
321 real(rp), intent(in) :: ddens_ (elem%np,lmesh%nea)
322
323 integer :: ke
324 real(rp) :: q0, q1
325 real(rp) :: q(elem%np)
326 real(rp) :: dens(elem%np)
327 !--------------------------------------------------------
328
329 !$omp parallel do private( dens, Q, Q0, Q1 )
330 do ke=lmesh%NeS, lmesh%NeE
331 dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
332 q(:) = max( 0.0_rp, qtrc_(:,ke) )
333
334 q0 = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * dens(:) * qtrc_(:,ke) )
335 q1 = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * dens(:) * q(:) )
336 qtrc_(:,ke) = q0 / ( q1 + 1.0e-32_rp ) * q(:)
337 end do
338
339 return
340 end subroutine atm_dyn_dgm_trcadvect3d_tmar
341
342!OCL SERIAL
344 MFLX_x_tavg, MFLX_y_tavg, MFLX_z_tavg, alph_dens_M, alph_dens_P, &
345 DDENS, MOMX, MOMY, MOMZ, DPRES, DENS_hyd, PRES_hyd, &
346 Rtot, CVtot, CPtot, &
347 lmesh, elem, rkstage, tavg_weight_h, tavg_weight_v, is_hevi )
348
349 implicit none
350 class(localmesh3d), intent(in) :: lmesh
351 class(elementbase3d), intent(in) :: elem
352 real(rp), intent(inout) :: mflx_x_tavg(elem%np,lmesh%nea)
353 real(rp), intent(inout) :: mflx_y_tavg(elem%np,lmesh%nea)
354 real(rp), intent(inout) :: mflx_z_tavg(elem%np,lmesh%nea)
355 real(rp), intent(inout) :: alph_dens_m(elem%nfptot,lmesh%ne)
356 real(rp), intent(inout) :: alph_dens_p(elem%nfptot,lmesh%ne)
357 real(rp), intent(in) :: ddens(elem%np,lmesh%nea)
358 real(rp), intent(in) :: momx(elem%np,lmesh%nea)
359 real(rp), intent(in) :: momy(elem%np,lmesh%nea)
360 real(rp), intent(in) :: momz(elem%np,lmesh%nea)
361 real(rp), intent(in) :: dpres(elem%np,lmesh%nea)
362 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
363 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
364 real(rp), intent(in) :: rtot(elem%np,lmesh%nea)
365 real(rp), intent(in) :: cvtot(elem%np,lmesh%nea)
366 real(rp), intent(in) :: cptot(elem%np,lmesh%nea)
367 integer, intent(in) :: rkstage
368 real(rp), intent(in) :: tavg_weight_h
369 real(rp), intent(in) :: tavg_weight_v
370 logical, intent(in) :: is_hevi
371
372 integer :: ke
373 !--------------------------------------------------------------
374
375 !$omp parallel do private(ke)
376 do ke=lmesh%NeS, lmesh%NeE
377 if (rkstage == 1) then
378 mflx_x_tavg(:,ke) = tavg_weight_h * momx(:,ke)
379 mflx_y_tavg(:,ke) = tavg_weight_h * momy(:,ke)
380 mflx_z_tavg(:,ke) = tavg_weight_v * momz(:,ke)
381 alph_dens_m(:,ke) = 0.0_rp
382 alph_dens_p(:,ke) = 0.0_rp
383 else
384 mflx_x_tavg(:,ke) = mflx_x_tavg(:,ke) + tavg_weight_h * momx(:,ke)
385 mflx_y_tavg(:,ke) = mflx_y_tavg(:,ke) + tavg_weight_h * momy(:,ke)
386 mflx_z_tavg(:,ke) = mflx_z_tavg(:,ke) + tavg_weight_v * momz(:,ke)
387 end if
388 end do
389
390 call atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn( &
391 alph_dens_m, alph_dens_p, & ! (out)
392 ddens, momx, momy, momz, dpres, dens_hyd, pres_hyd, & ! (in)
393 rtot, cvtot, cptot, & ! (in)
394 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
395 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
396 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, & ! (in)
397 lmesh, lmesh%lcmesh2D, lmesh%refElem3D, lmesh%lcmesh2D%refElem2D, & ! (in)
398 tavg_weight_h, tavg_weight_v, is_hevi ) ! (in)
399
400 return
402
403!OCL SERIAL
404 subroutine atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_advtest( alph_dens_M, alph_dens_P, &
405 DDENS_, MOMX_, MOMY_, MOMZ_, DENS_hyd, &
406 Gsqrt, nx, ny, nz, vmapM, vmapP, lmesh, elem )
407
408 use scale_const, only: &
409 grav => const_grav, &
410 rdry => const_rdry, &
411 cpdry => const_cpdry, &
412 cvdry => const_cvdry, &
413 pres00 => const_pre00
414
415 implicit none
416
417 class(localmesh3d), intent(in) :: lmesh
418 class(elementbase3d), intent(in) :: elem
419 real(rp), intent(inout) :: alph_dens_m(elem%nfptot*lmesh%ne)
420 real(rp), intent(inout) :: alph_dens_p(elem%nfptot*lmesh%ne)
421 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
422 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
423 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
424 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
425 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
426 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
427 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
428 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
429 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
430 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
431 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
432
433 integer :: i, ip, im
434 real(rp) :: velp, velm, alpha, densm, densp
435 !------------------------------------------------------------------------
436
437
438 !$omp parallel do private( &
439 !$omp iM, iP, VelP, VelM, alpha, densM, densP )
440 do i=1, elem%NfpTot*lmesh%Ne
441 im = vmapm(i); ip = vmapp(i)
442
443 densm = ddens_(im) + dens_hyd(im)
444 densp = ddens_(ip) + dens_hyd(ip)
445
446 velm = ( momx_(im) * nx(i) + momy_(im) * ny(i) + momz_(im) * nz(i) ) / densm
447 velp = ( momx_(ip) * nx(i) + momy_(ip) * ny(i) + momz_(ip) * nz(i) ) / densp
448
449 alpha = max( abs(velm), abs(velp) )
450 alph_dens_m(i) = alpha * densm * gsqrt(im)
451 alph_dens_p(i) = alpha * densp * gsqrt(ip)
452 end do
453
454 return
456
457!-- private ---
458
459!OCL SERIAL
460 subroutine atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn( alph_dens_M, alph_dens_P, &
461 DDENS_, MOMX_, MOMY_, MOMZ_, DPRES_, DENS_hyd, PRES_hyd, &
462 Rtot, CVtot, CPtot, &
463 Gsqrt, G11, G12, G22, nx, ny, nz, vmapM, vmapP, iM2Dto3D, &
464 lmesh, lmesh2D, elem, elem2D, &
465 tavg_weight_h, tavg_weight_v, is_hevi )
466
467 use scale_const, only: &
468 grav => const_grav, &
469 rdry => const_rdry, &
470 cpdry => const_cpdry, &
471 cvdry => const_cvdry, &
472 pres00 => const_pre00
473
474 implicit none
475
476 class(localmesh3d), intent(in) :: lmesh
477 class(localmesh2d), intent(in) :: lmesh2d
478 class(elementbase3d), intent(in) :: elem
479 class(elementbase2d), intent(in) :: elem2d
480 real(rp), intent(inout) :: alph_dens_m(elem%nfptot,lmesh%ne)
481 real(rp), intent(inout) :: alph_dens_p(elem%nfptot,lmesh%ne)
482 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
483 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
484 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
485 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
486 real(rp), intent(in) :: dpres_(elem%np*lmesh%nea)
487 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
488 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
489 real(rp), intent(in) :: rtot(elem%np*lmesh%nea)
490 real(rp), intent(in) :: cvtot(elem%np*lmesh%nea)
491 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
492 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
493 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
494 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
495 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
496 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
497 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
498 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
499 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
500 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
501 integer, intent(in) :: im2dto3d(elem%nfptot)
502 real(rp), intent(in) :: tavg_weight_h
503 real(rp), intent(in) :: tavg_weight_v
504 logical, intent(in) :: is_hevi
505
506 integer :: ke, ip(elem%nfptot), im(elem%nfptot)
507 integer :: ke2d
508 real(rp) :: velp(elem%nfptot), velm(elem%nfptot), alpha(elem%nfptot)
509 real(rp) :: densm(elem%nfptot), densp(elem%nfptot)
510 real(rp) :: gamm, rgamm
511 real(rp) :: tavg_weight(elem%nfptot)
512
513 real(rp) :: gnn_m(elem%nfptot), gnn_p(elem%nfptot)
514 !------------------------------------------------------------------------
515
516 gamm = cpdry/cvdry
517 rgamm = cvdry/cpdry
518
519 !$omp parallel do private( &
520 !$omp ke, iM, iP, ke2D, VelP, VelM, alpha, densM, densP, tavg_weight, &
521 !$omp Gnn_M, Gnn_P )
522 do ke=lmesh%NeS, lmesh%NeE
523 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
524 ke2d = lmesh%EMap3Dto2D(ke)
525
526 densm(:) = ddens_(im) + dens_hyd(im)
527 densp(:) = ddens_(ip) + dens_hyd(ip)
528
529 if ( is_hevi ) then
530 gnn_m(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) )
531 gnn_p(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) )
532 else
533 gnn_m(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) ) &
534 + abs(nz(:,ke))
535 gnn_p(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) ) &
536 + abs(nz(:,ke))
537 end if
538
539 velm(:) = ( momx_(im)*nx(:,ke) + momy_(im)*ny(:,ke) + momz_(im)*nz(:,ke) ) / densm(:)
540 velp(:) = ( momx_(ip)*nx(:,ke) + momy_(ip)*ny(:,ke) + momz_(ip)*nz(:,ke) ) / densp(:)
541
542 alpha(:) = max( sqrt( gnn_m(:) * gamm * ( pres_hyd(im) + dpres_(im) ) / densm(:) ) + abs(velm(:)), &
543 sqrt( gnn_p(:) * gamm * ( pres_hyd(ip) + dpres_(ip) ) / densp(:) ) + abs(velp(:)) )
544 tavg_weight = tavg_weight_h * ( abs(nx(:,ke)) + abs(ny(:,ke)) ) + tavg_weight_v * abs(nz(:,ke))
545
546 alph_dens_m(:,ke) = alph_dens_m(:,ke) + tavg_weight * alpha(:) * densm(:) * gsqrt(im)
547 alph_dens_p(:,ke) = alph_dens_p(:,ke) + tavg_weight * alpha(:) * densp(:) * gsqrt(ip)
548 end do
549
550 return
551 end subroutine atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn
552
553!OCL SERIAL
554 subroutine atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc( &
555 del_flux, & ! (out)
556 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef, & ! (in)
557 gsqrt, gsqrth, g13, g23, nx, ny, nz, j, fscale, & ! (in)
558 vmapm, vmapp, im2dto3d, faceintmat, & ! (in)
559 lmesh, elem, lmesh2d, elem2d ) ! (in)
560
561 implicit none
562
563 class(localmesh3d), intent(in) :: lmesh
564 class(elementbase3d), intent(in) :: elem
565 class(localmesh2d), intent(in) :: lmesh2d
566 class(elementbase2d), intent(in) :: elem2d
567 real(rp), intent(out) :: del_flux(elem%nfptot,lmesh%ne)
568 real(rp), intent(in) :: qtrc_(elem%np*lmesh%nea)
569 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
570 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
571 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
572 real(rp), intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
573 real(rp), intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
574 real(rp), intent(in) :: fct_coef(elem%np*lmesh%nea)
575 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
576 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
577 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
578 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
579 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
580 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
581 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
582 real(rp), intent(in) :: j(elem%np*lmesh%ne)
583 real(rp), intent(in) :: fscale(elem%nfptot,lmesh%ne)
584 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
585 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
586 integer, intent(in) :: im2dto3d(elem%nfptot)
587 type(sparsemat), intent(in) :: faceintmat
588
589 integer :: ke, i, ip(elem%nfptot), im(elem%nfptot)
590 integer :: ke2d
591 real(rp) :: momflxp(elem%nfptot), momflxm(elem%nfptot), alpha(elem%nfptot)
592 real(rp) :: qtrc_p(elem%nfptot), qtrc_m(elem%nfptot)
593 real(rp) :: gsqrtmomx_p(elem%nfptot), gsqrtmomx_m(elem%nfptot)
594 real(rp) :: gsqrtmomy_p(elem%nfptot), gsqrtmomy_m(elem%nfptot)
595 real(rp) :: gsqrtmomz_p(elem%nfptot), gsqrtmomz_m(elem%nfptot)
596 real(rp) :: gsqrt_p(elem%nfptot), gsqrt_m(elem%nfptot)
597 real(rp) :: gsqrtv_p(elem%nfptot), gsqrtv_m(elem%nfptot)
598 real(rp) :: g13_m(elem%nfptot), g13_p(elem%nfptot)
599 real(rp) :: g23_m(elem%nfptot), g23_p(elem%nfptot)
600
601 real(rp) :: r_m(elem%nfptot), r_p(elem%nfptot)
602 real(rp) :: numflux(elem%nfptot)
603 real(rp) :: outward_flux_tmp(elem%nfaces)
604 integer :: f, p, fp, is, ie
605 !------------------------------------------------------------------------
606
607 !$omp parallel do private( &
608 !$omp ke, iM, iP, ke2D, &
609 !$omp alpha, MomFlxM, MomFlxP, R_M, R_P, numflux, outward_flux_tmp, &
610 !$omp f, p, fp, &
611 !$omp GsqrtMOMX_M, GsqrtMOMX_P, GsqrtMOMY_M, GsqrtMOMY_P, GsqrtMOMZ_M, GsqrtMOMZ_P, &
612 !$omp QTRC_M, QTRC_P, Gsqrt_P, Gsqrt_M, GsqrtV_P, GsqrtV_M, G13_P, G13_M, G23_P, G23_M )
613 do ke=lmesh%NeS, lmesh%NeE
614 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
615 ke2d = lmesh%EMap3Dto2D(ke)
616
617 gsqrt_m(:) = gsqrt(im)
618 gsqrt_p(:) = gsqrt(ip)
619 gsqrtv_m(:) = gsqrt_m(:) / gsqrth(im2dto3d(:),ke2d)
620 gsqrtv_p(:) = gsqrt_p(:) / gsqrth(im2dto3d(:),ke2d)
621
622 g13_m(:) = g13(im)
623 g13_p(:) = g13(ip)
624 g23_m(:) = g23(im)
625 g23_p(:) = g23(ip)
626
627 r_m(:) = fct_coef(im)
628 r_p(:) = fct_coef(ip)
629
630 qtrc_m(:) = qtrc_(im)
631 qtrc_p(:) = qtrc_(ip)
632 gsqrtmomx_m(:) = gsqrt_m(:) * momx_(im)
633 gsqrtmomx_p(:) = gsqrt_p(:) * momx_(ip)
634 gsqrtmomy_m(:) = gsqrt_m(:) * momy_(im)
635 gsqrtmomy_p(:) = gsqrt_p(:) * momy_(ip)
636 gsqrtmomz_m(:) = gsqrt_m(:) * momz_(im)
637 gsqrtmomz_p(:) = gsqrt_p(:) * momz_(ip)
638
639 momflxm(:) = ( gsqrtmomx_m(:) * nx(:,ke) + gsqrtmomy_m(:) * ny(:,ke) &
640 + ( ( gsqrtmomz_m(:) / gsqrtv_m(:) &
641 + g13_m(:) * gsqrtmomx_m(:) + g23_m(:) * gsqrtmomy_m(:) ) * nz(:,ke) ) &
642 )
643 momflxp(:) = ( gsqrtmomx_p(:) * nx(:,ke) + gsqrtmomy_p(:) * ny(:,ke) &
644 + ( ( gsqrtmomz_p(:) / gsqrtv_p(:) &
645 + g13_p(:) * gsqrtmomx_p(:) + g23_p(:) * gsqrtmomy_p(:) ) * nz(:,ke) ) &
646 )
647
648 numflux(:) = 0.5_rp * ( ( qtrc_p(:) * momflxp(:) + qtrc_m(:) * momflxm(:) ) &
649 - alphdens_p(:,ke) * qtrc_p(:) + alphdens_m(:,ke) * qtrc_m(:) )
650 ! alpha(:) = max( abs(MomFlxM(:)), abs(MomFlxP(:)) )
651 ! numflux(:) = 0.5_RP * ( ( QTRC_P(:) * MomFlxP(:) + QTRC_M(:) * MomFlxM(:) ) &
652 ! - alpha(:) * (QTRC_P(:) - QTRC_M(:) ) )
653
654 call sparsemat_matmul( faceintmat, j(im) * fscale(:,ke) * numflux(:), outward_flux_tmp )
655 do f=1, elem%Nfaces_h
656 do p=1, elem%Nfp_h
657 fp = p + (f-1)*elem%Nfp_h
658 del_flux(fp,ke) = lmesh%Fscale(fp,ke) * &
659 ( numflux(fp) * 0.5_rp * ( r_p(fp) + r_m(fp) - ( r_p(fp) - r_m(fp) ) * sign( 1.0_rp, outward_flux_tmp(f) ) ) &
660 - qtrc_m(fp) * momflxm(fp) )
661 end do
662 end do
663 do f=1, elem%Nfaces_v
664 do p=1, elem%Nfp_v
665 fp = p + (f-1)*elem%Nfp_v + elem%Nfaces_h * elem%Nfp_h
666 del_flux(fp,ke) = lmesh%Fscale(fp,ke) * &
667 ( numflux(fp) * 0.5_rp * ( r_p(fp) + r_m(fp) - ( r_p(fp) - r_m(fp) ) * sign( 1.0_rp, outward_flux_tmp(elem%Nfaces_h+f) ) ) &
668 - qtrc_m(fp) * momflxm(fp) )
669 end do
670 end do
671 end do
672
673 return
674 end subroutine atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc
675
676
677!OCL SERIAL
678 subroutine atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc( &
679 net_outward_flux, & ! (out)
680 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, & ! (in)
681 gsqrt, gsqrth, g13, g23, nx, ny, nz, j, fscale, & ! (in)
682 vmapm, vmapp, im2dto3d, faceintmat, & ! (in)
683 lmesh, elem, lmesh2d, elem2d ) ! (in)
684
685 implicit none
686
687 class(localmesh3d), intent(in) :: lmesh
688 class(elementbase3d), intent(in) :: elem
689 class(localmesh2d), intent(in) :: lmesh2d
690 class(elementbase2d), intent(in) :: elem2d
691 real(rp), intent(out) :: net_outward_flux(lmesh%ne)
692 real(rp), intent(in) :: qtrc_(elem%np*lmesh%nea)
693 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
694 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
695 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
696 real(rp), intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
697 real(rp), intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
698 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
699 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
700 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
701 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
702 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
703 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
704 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
705 real(rp), intent(in) :: j(elem%np*lmesh%ne)
706 real(rp), intent(in) :: fscale(elem%nfptot,lmesh%ne)
707 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
708 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
709 integer, intent(in) :: im2dto3d(elem%nfptot)
710 type(sparsemat), intent(in) :: faceintmat
711
712 integer :: ke, i, ip(elem%nfptot), im(elem%nfptot)
713 integer :: ke2d
714 real(rp) :: momflxp(elem%nfptot), momflxm(elem%nfptot), alpha(elem%nfptot)
715 real(rp) :: qtrc_p(elem%nfptot), qtrc_m(elem%nfptot)
716 real(rp) :: gsqrtmomx_p(elem%nfptot), gsqrtmomx_m(elem%nfptot)
717 real(rp) :: gsqrtmomy_p(elem%nfptot), gsqrtmomy_m(elem%nfptot)
718 real(rp) :: gsqrtmomz_p(elem%nfptot), gsqrtmomz_m(elem%nfptot)
719 real(rp) :: gsqrt_p(elem%nfptot), gsqrt_m(elem%nfptot)
720 real(rp) :: gsqrtv_p(elem%nfptot), gsqrtv_m(elem%nfptot)
721 real(rp) :: g13_m(elem%nfptot), g13_p(elem%nfptot)
722 real(rp) :: g23_m(elem%nfptot), g23_p(elem%nfptot)
723
724 real(rp) :: numflux(elem%nfptot)
725 real(rp) :: outward_flux_tmp(elem%nfaces)
726 !------------------------------------------------------------------------
727
728 !$omp parallel do private( &
729 !$omp ke, iM, iP, ke2D, &
730 !$omp alpha, MomFlxM, MomFlxP, numflux, outward_flux_tmp, &
731 !$omp GsqrtMOMX_M, GsqrtMOMX_P, GsqrtMOMY_M, GsqrtMOMY_P, GsqrtMOMZ_M, GsqrtMOMZ_P, &
732 !$omp QTRC_M, QTRC_P, Gsqrt_P, Gsqrt_M, GsqrtV_P, GsqrtV_M, G13_P, G13_M, G23_P, G23_M )
733 do ke=lmesh%NeS, lmesh%NeE
734 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
735 ke2d = lmesh%EMap3Dto2D(ke)
736
737 gsqrt_m(:) = gsqrt(im)
738 gsqrt_p(:) = gsqrt(ip)
739 gsqrtv_m(:) = gsqrt_m(:) / gsqrth(im2dto3d(:),ke2d)
740 gsqrtv_p(:) = gsqrt_p(:) / gsqrth(im2dto3d(:),ke2d)
741
742 g13_m(:) = g13(im)
743 g13_p(:) = g13(ip)
744 g23_m(:) = g23(im)
745 g23_p(:) = g23(ip)
746
747 qtrc_m(:) = qtrc_(im)
748 qtrc_p(:) = qtrc_(ip)
749 gsqrtmomx_m(:) = gsqrt_m(:) * momx_(im)
750 gsqrtmomx_p(:) = gsqrt_p(:) * momx_(ip)
751 gsqrtmomy_m(:) = gsqrt_m(:) * momy_(im)
752 gsqrtmomy_p(:) = gsqrt_p(:) * momy_(ip)
753 gsqrtmomz_m(:) = gsqrt_m(:) * momz_(im)
754 gsqrtmomz_p(:) = gsqrt_p(:) * momz_(ip)
755
756 momflxm(:) = ( gsqrtmomx_m(:) * nx(:,ke) + gsqrtmomy_m(:) * ny(:,ke) &
757 + ( ( gsqrtmomz_m(:) / gsqrtv_m(:) &
758 + g13_m(:) * gsqrtmomx_m(:) + g23_m(:) * gsqrtmomy_m(:) ) * nz(:,ke) ) &
759 )
760 momflxp(:) = ( gsqrtmomx_p(:) * nx(:,ke) + gsqrtmomy_p(:) * ny(:,ke) &
761 + ( ( gsqrtmomz_p(:) / gsqrtv_p(:) &
762 + g13_p(:) * gsqrtmomx_p(:) + g23_p(:) * gsqrtmomy_p(:) ) * nz(:,ke) ) &
763 )
764
765 alpha(:) = max( abs(momflxm(:)), abs(momflxp(:)) )
766
767 numflux(:) = 0.5_rp * ( ( qtrc_p(:) * momflxp(:) + qtrc_m(:) * momflxm(:) ) &
768 - alphdens_p(:,ke) * qtrc_p(:) + alphdens_m(:,ke) * qtrc_m(:) )
769 ! numflux(:) = 0.5_RP * ( ( QTRC_P(:) * MomFlxP(:) + QTRC_M(:) * MomFlxM(:) ) &
770 ! - alpha(:) * ( QTRC_P(:) - QTRC_M(:) ) )
771
772 call sparsemat_matmul( faceintmat, j(im) * fscale(:,ke) * numflux(:), outward_flux_tmp )
773 net_outward_flux(ke) = sum( max( 0.0_rp, outward_flux_tmp(:) ) )
774 end do
775
776 return
777 end subroutine atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc
778
module FElib / Fluid dyn solver / Atmosphere / Tracer advection
subroutine, public atm_dyn_dgm_trcadvect3d_tmar(qtrc_, dens_hyd, ddens_, lmesh, elem, lmesh2d, elem2d)
Second Step of limiter in which nonlinear truncation and mass aware rescaling (TMAR)
subroutine, public atm_dyn_dgm_trcadvect3d_save_massflux(mflx_x_tavg, mflx_y_tavg, mflx_z_tavg, alph_dens_m, alph_dens_p, ddens, momx, momy, momz, dpres, dens_hyd, pres_hyd, rtot, cvtot, cptot, lmesh, elem, rkstage, tavg_weight_h, tavg_weight_v, is_hevi)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef(fct_coef, qtrc_, momx_, momy_, momz_, rhoq_tp_, alphdens_m, alphdens_p, dens_hyd, ddens_, ddens0_, rk_c_ssm1, dt, faceintmat, lmesh, elem, lmesh2d, elem2d, disable_limiter)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_cal_tend(qtrc_dt, qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef, rhoq_tp, element3d_operation, faceintmat, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_advtest(alph_dens_m, alph_dens_p, ddens_, momx_, momy_, momz_, dens_hyd, gsqrt, nx, ny, nz, vmapm, vmapp, lmesh, elem)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_init(mesh, faceintmat)
module FElib / Element / Base
module FElib / Element / hexahedron
module FElib / Element / Operation / Base
module FElib / Mesh / Local 2D
module FElib / Mesh / Local 3D
module FElib / Mesh / Base 2D
module FElib / Mesh / Base 3D
module FElib / Data / base
Module common / Polynomial.
real(rp) function, dimension(nord+1), public polynomial_gengausslobattoptintweight(nord)
A function to calculate the Gauss-Lobbato weights.
Module common / sparsemat.
Derived type representing a 2D reference element.
Derived type representing a 3D reference element.
Derived type representing a hexahedral element.
Derived type representing a local mesh for 2D domain.
Derived type to manage a local 3D computational domain.
Derived type representing a field with 3D local mesh.
Derived type to manage a computational mesh (base type for 2D domain)
Derived type to manage a computational mesh (base type for 3D domain)
Derived type representing a field with 3D mesh.
Derived type to manage a sparse matrix.