FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_bnd.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module FElib / Fluid dyn solver / Atmosphere / Boundary
3!!
4!! @par Description
5!! A module for setting halo data based on boundary conditions.
6!!
7!! @author Yuta Kawai, Team SCALE
8!<
9!-------------------------------------------------------------------------------
10#include "scaleFElib.h"
11module scale_atm_dyn_dgm_bnd
12 !-----------------------------------------------------------------------------
13 !
14 !++ Used modules
15 !
16 use scale_precision
17 use scale_io
18 use scale_prc
19 use scale_prof
20 use scale_const, only: &
21 grav => const_grav, &
22 rdry => const_rdry, &
23 cpdry => const_cpdry, &
24 cvdry => const_cvdry, &
25 pres00 => const_pre00
26
27 use scale_element_base, only: &
28 elementbase, elementbase2d, elementbase3d
29 use scale_mesh_base, only: meshbase
30 use scale_localmesh_base, only: localmeshbase
31
32 use scale_mesh_base3d, only: meshbase3d
33 use scale_mesh_base2d, only: meshbase2d
34 use scale_localmesh_3d, only: localmesh3d
35 use scale_localmesh_2d, only: localmesh2d
36
37 use scale_localmeshfield_base, only: localmeshfield3d
38 use scale_meshfield_base, only: meshfield3d
39
40 use scale_mesh_bndinfo, only: meshbndinfo
41
42 use scale_atm_dyn_dgm_nonhydro3d_common, only: &
43 dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
44 momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id, &
45 momz_vid => prgvar_momz_id, &
46 prgvar_num
47
48 !-----------------------------------------------------------------------------
49 implicit none
50 private
51 !-----------------------------------------------------------------------------
52 !
53 !++ Public parameter & type & procedures
54 !
55
56 !> A derived type useful for apply boundary conditions
57 type, public :: atmdynbnd
58 type(meshbndinfo), allocatable :: velbc_list(:)
59 type(meshbndinfo), allocatable :: thermalbc_list(:)
60
61 integer, allocatable:: velbc_ids(:)
62 integer, allocatable :: thermalbc_ids(:)
63 real(rp), allocatable :: thermal_fixval(:)
64 contains
65 procedure :: init => atmos_dyn_bnd_setup
66 procedure :: final => atmos_dyn_bnd_finalize
67 procedure :: setbcinfo => atmos_dyn_bnd_setbcinfo
68 procedure :: applybc_progvars_lc => atmos_dyn_bnd_applybc_prgvars_lc
69 procedure :: applybc_numdiff_odd_lc => atmos_dyn_bnd_applybc_numdiff_odd_lc
70 procedure :: applybc_numdiff_even_lc => atmos_dyn_bnd_applybc_numdiff_even_lc
71 procedure :: applybc_grad_tbvars_lc => atmos_dyn_bnd_applybc_tbvars_lc
72 procedure :: applybc_grad_tbstress_lc => atmos_dyn_bnd_applybc_tbstress_lc
73 procedure :: inquire_bound_flag => atmos_dyn_bnd_inquire_bound_flag
74 end type
75
76 !-----------------------------------------------------------------------------
77 !
78 !++ Private procedures & variables
79 !
80 !-------------------
81
82 private :: bnd_init_lc
83
84 integer, parameter :: dombnd_south_id = 1
85 integer, parameter :: dombnd_east_id = 2
86 integer, parameter :: dombnd_north_id = 3
87 integer, parameter :: dombnd_west_id = 4
88 integer, parameter :: dombnd_btm_id = 5
89 integer, parameter :: dombnd_top_id = 6
90 integer, parameter :: dom_bnd_num = 6
91
92contains
93
94!> Setup an object to manage boundary conditions with dynamical process
95!OCL SERIAL
96 subroutine atmos_dyn_bnd_setup( this )
97 use scale_const, only: &
98 undef8 => const_undef8
99 use scale_mesh_bndinfo, only: &
100 bnd_type_nospec_name, &
101 bnd_type_periodic_id, &
102 bndtype_nametoid
103 implicit none
104
105 class(atmdynbnd), intent(inout) :: this
106
107 character(len=H_SHORT) :: btm_vel_bc !< Velocity BC at bottom boundary
108 character(len=H_SHORT) :: top_vel_bc !< Velocity BC at top boundary
109 character(len=H_SHORT) :: north_vel_bc !< Velocity BC at northern boundary
110 character(len=H_SHORT) :: south_vel_bc !< Velocity BC at southern boundary
111 character(len=H_SHORT) :: east_vel_bc !< Velocity BC at east boundary
112 character(len=H_SHORT) :: west_vel_bc !< Velocity BC at west boundary
113
114 character(len=H_SHORT) :: btm_thermal_bc !< Thermal BC at bottom boundary
115 character(len=H_SHORT) :: top_thermal_bc !< Thermal BC at top boundary
116 character(len=H_SHORT) :: north_thermal_bc !< Thermal BC at northern boundary
117 character(len=H_SHORT) :: south_thermal_bc !< Thermal BC at southern boundary
118 character(len=H_SHORT) :: east_thermal_bc !< Thermal BC at east boundary
119 character(len=H_SHORT) :: west_thermal_bc !< Thermal BC at west boundary
120
121 real(rp) :: btm_thermal_fixval !< Fixed value with thermal BC at bottom boundary
122 real(rp) :: top_thermal_fixval !< Fixed value with thermal BC at top boundary
123 real(rp) :: north_thermal_fixval !< Fixed value with thermal BC at north boundary
124 real(rp) :: south_thermal_fixval !< Fixed value with thermal BC at south boundary
125 real(rp) :: east_thermal_fixval !< Fixed value with thermal BC at east boundary
126 real(rp) :: west_thermal_fixval !< Fixed value with thermal BC at west boundary
127
128 namelist /param_atmos_dyn_bnd/ &
129 btm_vel_bc, top_vel_bc, north_vel_bc, south_vel_bc, east_vel_bc, west_vel_bc, &
130 btm_thermal_bc, top_thermal_bc, north_thermal_bc, south_thermal_bc, east_thermal_bc, west_thermal_bc, &
131 btm_thermal_fixval, top_thermal_fixval, north_thermal_fixval, south_thermal_fixval, east_thermal_fixval, west_thermal_fixval
132
133 integer :: ierr
134 !-----------------------------------------------
135
136 btm_vel_bc = bnd_type_nospec_name
137 top_vel_bc = bnd_type_nospec_name
138 north_vel_bc = bnd_type_nospec_name
139 south_vel_bc = bnd_type_nospec_name
140 east_vel_bc = bnd_type_nospec_name
141 west_vel_bc = bnd_type_nospec_name
142
143 btm_thermal_bc = bnd_type_nospec_name
144 top_thermal_bc = bnd_type_nospec_name
145 north_thermal_bc = bnd_type_nospec_name
146 south_thermal_bc = bnd_type_nospec_name
147 east_thermal_bc = bnd_type_nospec_name
148 west_thermal_bc = bnd_type_nospec_name
149
150 btm_thermal_fixval = undef8
151 top_thermal_fixval = undef8
152 north_thermal_fixval = undef8
153 south_thermal_fixval = undef8
154 east_thermal_fixval = undef8
155 west_thermal_fixval = undef8
156
157 rewind(io_fid_conf)
158 read(io_fid_conf,nml=param_atmos_dyn_bnd,iostat=ierr)
159 if( ierr < 0 ) then !--- missing
160 log_info("ATMOS_dyn_bnd_setup",*) 'Not found namelist. Default used.'
161 elseif( ierr > 0 ) then !--- fatal error
162 log_error("ATMOS_dyn_bnd_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_BND. Check!'
163 call prc_abort
164 endif
165 log_nml(param_atmos_dyn_bnd)
166
167 !--
168 allocate( this%velBC_ids(dom_bnd_num) )
169 this%velBC_ids(dombnd_btm_id) = bndtype_nametoid(btm_vel_bc)
170 this%velBC_ids(dombnd_top_id) = bndtype_nametoid(top_vel_bc)
171 this%velBC_ids(dombnd_north_id) = bndtype_nametoid(north_vel_bc)
172 this%velBC_ids(dombnd_south_id) = bndtype_nametoid(south_vel_bc)
173 this%velBC_ids(dombnd_east_id) = bndtype_nametoid(east_vel_bc)
174 this%velBC_ids(dombnd_west_id) = bndtype_nametoid(west_vel_bc)
175
176 allocate( this%thermalBC_ids(dom_bnd_num) )
177 this%thermalBC_ids(dombnd_btm_id) = bndtype_nametoid(btm_thermal_bc)
178 this%thermalBC_ids(dombnd_top_id) = bndtype_nametoid(top_thermal_bc)
179 this%thermalBC_ids(dombnd_north_id) = bndtype_nametoid(north_thermal_bc)
180 this%thermalBC_ids(dombnd_south_id) = bndtype_nametoid(south_thermal_bc)
181 this%thermalBC_ids(dombnd_east_id) = bndtype_nametoid(east_thermal_bc)
182 this%thermalBC_ids(dombnd_west_id) = bndtype_nametoid(west_thermal_bc)
183
184 allocate( this%thermal_fixval(dom_bnd_num) )
185 this%thermal_fixval(dombnd_btm_id) = btm_thermal_fixval
186 this%thermal_fixval(dombnd_top_id) = top_thermal_fixval
187 this%thermal_fixval(dombnd_north_id) = north_thermal_fixval
188 this%thermal_fixval(dombnd_south_id) = south_thermal_fixval
189 this%thermal_fixval(dombnd_east_id) = east_thermal_fixval
190 this%thermal_fixval(dombnd_west_id) = west_thermal_fixval
191
192 !------
193
194 return
195 end subroutine atmos_dyn_bnd_setup
196
197!> Finalize an object to manage boundary conditions with dynamical process
198!OCL SERIAL
199 subroutine atmos_dyn_bnd_finalize( this )
200 implicit none
201 class(atmdynbnd), intent(inout) :: this
202
203 integer :: n
204 !--------------------------------------
205
206 if ( allocated(this%VelBC_list) ) then
207 do n=1, size(this%VelBC_list)
208 call this%VelBC_list(n)%Final()
209 call this%ThermalBC_list(n)%Final()
210 end do
211 deallocate( this%VelBC_list, this%ThermalBC_list )
212
213 deallocate( this%velBC_ids, this%thermalBC_ids )
214 end if
215
216 return
217 end subroutine atmos_dyn_bnd_finalize
218
219!OCL SERIAL
220 subroutine atmos_dyn_bnd_setbcinfo( this, mesh )
221
222 implicit none
223
224 class(atmdynbnd), intent(inout) :: this
225 class(meshbase), target, intent(in) :: mesh
226
227 integer :: n
228 integer :: b
229 class(localmeshbase), pointer :: ptr_lcmesh
230 class(localmesh3d), pointer :: lcmesh3d
231 !--------------------------------------------------
232
233
234 allocate( this%VelBC_list(mesh%LOCAL_MESH_NUM) )
235 allocate( this%ThermalBC_list(mesh%LOCAL_MESH_NUM) )
236
237 nullify( lcmesh3d )
238
239 do n=1, mesh%LOCAL_MESH_NUM
240 call mesh%GetLocalMesh( n, ptr_lcmesh )
241 select type (ptr_lcmesh)
242 type is (localmesh3d)
243 lcmesh3d => ptr_lcmesh
244 end select
245
246 call bnd_init_lc( &
247 this%VelBC_list(n), this%ThermalBC_list(n), & ! (inout)
248 this%velBC_ids(:), this%thermalBC_ids(:), & ! (in)
249 this%thermal_fixval(:), & ! (in)
250 lcmesh3d%VMapB, mesh, lcmesh3d, lcmesh3d%refElem3D ) ! (in)
251 end do
252
253 return
254 end subroutine atmos_dyn_bnd_setbcinfo
255
256!> Set exterior values at element boundaries to evaluate inviscid numerical fluxes
257!!
258!OCL SERIAL
259 subroutine atmos_dyn_bnd_applybc_prgvars_lc( this, &
260 domID, & ! (in)
261 ddens, momx, momy, momz, therm, & ! (inout)
262 dens_hyd, pres_hyd, & ! (in)
263 gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
264 vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d ) ! (in)
265
266 use scale_mesh_bndinfo, only: &
267 bnd_type_slip_id, bnd_type_noslip_id, &
268 bnd_type_adiabat_id
269
270 implicit none
271
272 class(atmdynbnd), intent(in) :: this
273 integer, intent(in) :: domid
274 class(localmesh3d), intent(in) :: lmesh
275 class(elementbase3d), intent(in) :: elem
276 class(localmesh2d), intent(in) :: lmesh2d
277 class(elementbase2d), intent(in) :: elem2d
278 real(rp), intent(inout) :: ddens(elem%np*lmesh%nea)
279 real(rp), intent(inout) :: momx(elem%np*lmesh%nea)
280 real(rp), intent(inout) :: momy(elem%np*lmesh%nea)
281 real(rp), intent(inout) :: momz(elem%np*lmesh%nea)
282 real(rp), intent(inout) :: therm(elem%np*lmesh%nea)
283 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
284 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
285 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
286 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
287 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
288 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
289 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
290 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
291 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
292 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
293 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
294 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
295 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
296 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
297 integer, intent(in) :: vmapb(:)
298
299 integer :: p, ke, ke2d
300 integer :: i, i_, im, ip
301 real(rp) :: mom_normal
302
303 real(rp) :: momw
304 real(rp) :: gsqrtv, g11_, g12_, g22_
305 real(rp) :: fac
306 !-----------------------------------------------
307
308 !$omp parallel do collapse(2) private( &
309 !$omp ke, p, ke2D, i, i_, iM, iP, &
310 !$omp mom_normal, MOMW, GsqrtV, G11_, G12_, G22_, fac )
311 do ke=lmesh%NeS, lmesh%NeE
312 do p=1, elem%NfpTot
313 i = p + (ke-1)*elem%NfpTot
314 ip = vmapp(i)
315 i_ = ip - elem%Np*lmesh%NeE
316
317 if (i_ > 0) then
318 im = vmapm(i)
319
320 select case( this%VelBC_list(domid)%list(i_) )
321 case ( bnd_type_slip_id)
322 ke2d = lmesh%EMap3Dto2D(ke)
323
324 gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d)
325 g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
326 g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
327 g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
328
329 momw = momz(im) / gsqrtv &
330 + g13(im) * momx(im) + g23(im) * momy(im)
331 fac = nz(i) * gsqrtv**2 / ( 1.0_rp + g11_ * ( gsqrtv * g13(im) )**2 + 2.0_rp * g12_ * ( gsqrtv**2 * g13(im) * g23(im) ) + g22_ * ( gsqrtv * g23(im) )**2 )
332
333 mom_normal = momx(im) * nx(i) + momy(im) * ny(i) + momw * nz(i)
334 momx(ip) = momx(im) - 2.0_rp * mom_normal * ( nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) ) )
335 momy(ip) = momy(im) - 2.0_rp * mom_normal * ( ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) ) )
336 momz(ip) = momz(im) - 2.0_rp * mom_normal * fac / gsqrtv
337
338 case ( bnd_type_noslip_id )
339 momx(ip) = - momx(im)
340 momy(ip) = - momy(im)
341 momz(ip) = - momz(im)
342 end select
343 end if
344
345 end do
346 end do
347
348 return
349 end subroutine atmos_dyn_bnd_applybc_prgvars_lc
350
351!OCL SERIAL
352 subroutine atmos_dyn_bnd_applybc_numdiff_odd_lc( this, & ! (in)
353 gxvar, gyvar, gzvar, & ! (inout)
354 is_bound, & ! (out)
355 varid, domid, & ! (in)
356 nx, ny, nz, vmapm, vmapp, vmapb, lmesh, elem ) ! (in)
357
358 use scale_mesh_bndinfo, only: &
359 bnd_type_slip_id, bnd_type_noslip_id, &
360 bnd_type_adiabat_id
361
362 implicit none
363
364 class(atmdynbnd), intent(in) :: this
365 class(localmesh3d), intent(in) :: lmesh
366 class(elementbase3d), intent(in) :: elem
367 real(rp), intent(inout) :: gxvar(elem%np*lmesh%nea)
368 real(rp), intent(inout) :: gyvar(elem%np*lmesh%nea)
369 real(rp), intent(inout) :: gzvar(elem%np*lmesh%nea)
370 logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
371 integer, intent(in) :: varid
372 integer, intent(in) :: domid
373 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
374 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
375 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
376 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
377 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
378 integer, intent(in) :: vmapb(:)
379
380 integer :: i, i_, im, ip
381 real(rp) :: grad_normal
382 !-----------------------------------------------
383
384 do i=1, elem%NfpTot*lmesh%Ne
385 ip = vmapp(i)
386 i_ = ip - elem%Np*lmesh%NeE
387 is_bound(i) = .false.
388
389 if (i_ > 0) then
390 im = vmapm(i)
391 grad_normal = gxvar(im) * nx(i) + gyvar(im) * ny(i) + gzvar(im) * nz(i)
392
393 if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
394 select case(varid)
395 case(momx_vid)
396 gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
397 gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
398 case(momy_vid)
399 gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
400 gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
401 case(momz_vid)
402 gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
403 gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
404 end select
405 is_bound(i) = .true.
406
407 end if
408 if ( this%ThermalBC_list(domid)%list(i_) == bnd_type_adiabat_id ) then
409 select case(varid)
410 case(dens_vid, rhot_vid)
411 gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
412 gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
413 gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
414 end select
415 is_bound(i) = .true.
416
417 end if
418
419 end if
420
421 end do
422
423 return
424 end subroutine atmos_dyn_bnd_applybc_numdiff_odd_lc
425
426!OCL SERIAL
427 subroutine atmos_dyn_bnd_applybc_numdiff_even_lc( this, & ! (in)
428 var, & ! (inout)
429 is_bound, & ! (out)
430 varid, domid, & ! (in)
431 nx, ny, nz, vmapm, vmapp, vmapb, lmesh, elem ) ! (in)
432
433 use scale_mesh_bndinfo, only: &
434 bnd_type_slip_id, bnd_type_noslip_id, &
435 bnd_type_adiabat_id
436 implicit none
437
438 class(atmdynbnd), intent(in) :: this
439 class(localmesh3d), intent(in) :: lmesh
440 class(elementbase3d), intent(in) :: elem
441 real(rp), intent(inout) :: var(elem%np*lmesh%nea)
442 logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
443 integer, intent(in) :: varid
444 integer, intent(in) :: domid
445 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
446 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
447 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
448 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
449 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
450 integer, intent(in) :: vmapb(:)
451
452 integer :: i, i_, im, ip
453 real(rp) :: grad_normal
454 !-----------------------------------------------
455
456 do i=1, elem%NfpTot*lmesh%Ne
457 ip = vmapp(i)
458 i_ = ip - elem%Np*lmesh%NeE
459 is_bound(i) = .false.
460
461 if (i_ > 0) then
462 im = vmapm(i)
463
464 if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
465 select case(varid)
466 case(momx_vid)
467 grad_normal = var(im) * nx(i)
468 var(ip) = var(im) - 2.0_rp * grad_normal * nx(i)
469 case(momy_vid)
470 grad_normal = var(im) * ny(i)
471 var(ip) = var(im) - 2.0_rp * grad_normal * ny(i)
472 case(momz_vid)
473 grad_normal = var(im) * nz(i)
474 var(ip) = var(im) - 2.0_rp * grad_normal * nz(i)
475 end select
476 is_bound(i) = .true.
477
478 else if ( this%VelBC_list(domid)%list(i_) == bnd_type_noslip_id ) then
479 select case(varid)
480 case( momx_vid, momy_vid, momz_vid )
481 var(ip) = - var(im)
482 end select
483 is_bound(i) = .true.
484 end if
485
486 end if
487 end do
488
489 return
490 end subroutine atmos_dyn_bnd_applybc_numdiff_even_lc
491
492
493!> Set exterior values at boundaries for the turbulent schemes
494!!
495!OCL SERIAL
496 subroutine atmos_dyn_bnd_applybc_tbvars_lc( this, &
497 domID, & ! (in)
498 ddens, momx, momy, momz, pt, pres, & ! (inout)
499 dens_hyd, pres_hyd, rtot, cptot, & ! (in)
500 gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
501 vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d ) ! (in)
502
503 use scale_const, only: &
504 pres00 => const_pre00
505 use scale_mesh_bndinfo, only: &
506 bnd_type_slip_id, bnd_type_noslip_id, &
507 bnd_type_adiabat_id, bnd_type_fixval_id
508
509 implicit none
510
511 class(atmdynbnd), intent(in) :: this
512 integer, intent(in) :: domid
513 class(localmesh3d), intent(in) :: lmesh
514 class(elementbase3d), intent(in) :: elem
515 class(localmesh2d), intent(in) :: lmesh2d
516 class(elementbase2d), intent(in) :: elem2d
517 real(rp), intent(inout) :: ddens(elem%np*lmesh%nea)
518 real(rp), intent(inout) :: momx(elem%np*lmesh%nea)
519 real(rp), intent(inout) :: momy(elem%np*lmesh%nea)
520 real(rp), intent(inout) :: momz(elem%np*lmesh%nea)
521 real(rp), intent(inout) :: pt(elem%np*lmesh%nea)
522 real(rp), intent(inout) :: pres(elem%np*lmesh%nea)
523 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
524 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
525 real(rp), intent(in) :: rtot(elem%np*lmesh%nea)
526 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
527 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
528 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
529 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
530 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
531 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
532 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
533 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
534 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
535 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
536 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
537 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
538 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
539 integer, intent(in) :: vmapb(:)
540
541 integer :: p, ke, ke2d
542 integer :: i, i_, im, ip
543 real(rp) :: mom_normal
544
545 real(rp) :: momw
546 real(rp) :: temp_p, temp_b
547 real(rp) :: gsqrtv, g11_, g12_, g22_
548 real(rp) :: fac
549 !-----------------------------------------------
550
551 !$omp parallel do collapse(2) private( &
552 !$omp ke, p, ke2D, i, i_, iM, iP, &
553 !$omp mom_normal, MOMW, GsqrtV, G11_, G12_, G22_, fac, &
554 !$omp TEMP_P, TEMP_B )
555 do ke=lmesh%NeS, lmesh%NeE
556 do p=1, elem%NfpTot
557 i = p + (ke-1)*elem%NfpTot
558 ip = vmapp(i)
559 i_ = ip - elem%Np*lmesh%NeE
560
561 if (i_ > 0) then
562 im = vmapm(i)
563
564 select case( this%VelBC_list(domid)%list(i_) )
565 case ( bnd_type_slip_id)
566 ke2d = lmesh%EMap3Dto2D(ke)
567
568 gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d)
569 g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
570 g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
571 g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
572
573 momw = momz(im) / gsqrtv &
574 + g13(im) * momx(im) + g23(im) * momy(im)
575 fac = nz(i) * gsqrtv**2 / ( 1.0_rp + g11_ * ( gsqrtv * g13(im) )**2 + 2.0_rp * g12_ * ( gsqrtv**2 * g13(im) * g23(im) ) + g22_ * ( gsqrtv * g23(im) )**2 )
576
577 mom_normal = momx(im) * nx(i) + momy(im) * ny(i) + momw * nz(i)
578 momx(ip) = momx(im) - 2.0_rp * mom_normal * ( nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) ) )
579 momy(ip) = momy(im) - 2.0_rp * mom_normal * ( ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) ) )
580 momz(ip) = momz(im) - 2.0_rp * mom_normal * fac / gsqrtv
581
582 case ( bnd_type_noslip_id )
583 momx(ip) = - momx(im)
584 momy(ip) = - momy(im)
585 momz(ip) = - momz(im)
586 end select
587
588 select case( this%ThermalBC_list(domid)%list(i_) )
589 case ( bnd_type_fixval_id )
590 temp_b = this%ThermalBC_list(domid)%val(i_)
591 temp_p = 2.0_rp * temp_b - pres(im) / ( ( dens_hyd(im) + ddens(im) ) * rtot(im) )
592 pt(ip) = temp_p * ( pres00 / pres(im) )**(rtot(im)/cptot(im))
593 ddens(ip) = pres(im) / ( rtot(im) * temp_p ) - dens_hyd(im)
594 end select
595 end if
596
597 end do
598 end do
599
600 return
601 end subroutine atmos_dyn_bnd_applybc_tbvars_lc
602
603!> Set exterior shear-stress tensor and heat flux at boundaries for the turbulent schemes
604!!
605!! Note: This codes are tentatively implementated.
606!! We need to check whether the formulations is valid for the case when general vertical coordinate is introduced.
607!OCL SERIAL
608 subroutine atmos_dyn_bnd_applybc_tbstress_lc( this, &
609 domID, & ! (in)
610 t11, t12, t13, t21, t22, t23, t31, t32, t33, & ! (inout)
611 df1, df2, df3, & ! (in)
612 gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
613 vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d ) ! (in)
614
615 use scale_mesh_bndinfo, only: &
616 bnd_type_slip_id, bnd_type_noslip_id, &
617 bnd_type_adiabat_id, bnd_type_fixval_id
618
619 implicit none
620
621 class(atmdynbnd), intent(in) :: this
622 integer, intent(in) :: domid
623 class(localmesh3d), intent(in) :: lmesh
624 class(elementbase3d), intent(in) :: elem
625 class(localmesh2d), intent(in) :: lmesh2d
626 class(elementbase2d), intent(in) :: elem2d
627 real(rp), intent(inout) :: t11(elem%np*lmesh%nea), t12(elem%np*lmesh%nea), t13(elem%np*lmesh%nea)
628 real(rp), intent(inout) :: t21(elem%np*lmesh%nea), t22(elem%np*lmesh%nea), t23(elem%np*lmesh%nea)
629 real(rp), intent(inout) :: t31(elem%np*lmesh%nea), t32(elem%np*lmesh%nea), t33(elem%np*lmesh%nea)
630 real(rp), intent(inout) :: df1(elem%np*lmesh%nea), df2(elem%np*lmesh%nea), df3(elem%np*lmesh%nea)
631 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
632 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
633 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
634 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
635 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
636 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
637 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
638 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
639 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
640 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
641 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
642 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
643 integer, intent(in) :: vmapb(:)
644
645 integer :: p, ke, ke2d
646 integer :: i, i_, im, ip
647
648 real(rp) :: gsqrtv, g11_, g12_, g22_
649 real(rp) :: fac, normal_vec(3)
650
651 real(rp) :: stress_t_tmp, stress_t_tmpvec(3)
652 real(rp) :: vecw
653 real(rp) :: normal_flux
654 !-----------------------------------------------
655
656 !$omp parallel do collapse(2) private( &
657 !$omp ke, p, ke2D, i, i_, iM, iP, &
658 !$omp VECW, GsqrtV, G11_, G12_, G22_, fac, normal_vec, &
659 !$omp stress_t_tmp, stress_t_tmpvec, normal_flux )
660 do ke=lmesh%NeS, lmesh%NeE
661 do p=1, elem%NfpTot
662 i = p + (ke-1)*elem%NfpTot
663 ip = vmapp(i)
664 i_ = ip - elem%Np*lmesh%NeE
665
666 if (i_ > 0) then
667 ke2d = lmesh%EMap3Dto2D(ke)
668 im = vmapm(i)
669
670 gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d)
671 g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
672 g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
673 g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
674 fac = nz(i) * gsqrtv**2 / ( 1.0_rp + g11_ * ( gsqrtv * g13(im) )**2 + 2.0_rp * g12_ * ( gsqrtv**2 * g13(im) * g23(im) ) + g22_ * ( gsqrtv * g23(im) )**2 )
675 normal_vec(1) = nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) )
676 normal_vec(2) = ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) )
677 normal_vec(3) = fac / gsqrtv
678
679 select case( this%VelBC_list(domid)%list(i_) )
680 case ( bnd_type_slip_id)
681 vecw = t13(im) / gsqrtv + g13(im) * t11(im) + g23(im) * t12(im)
682 stress_t_tmpvec(1) = t11(im) * nx(i) + t12(im) * ny(i) + vecw * nz(i)
683
684 vecw = t23(im) / gsqrtv + g13(im) * t21(im) + g23(im) * t22(im)
685 stress_t_tmpvec(2) = t21(im) * nx(i) + t22(im) * ny(i) + vecw * nz(i)
686
687 vecw = t33(im) / gsqrtv
688 stress_t_tmpvec(3) = t31(im) * nx(i) + t32(im) * ny(i) + vecw * nz(i)
689
690 stress_t_tmp = sum( stress_t_tmpvec(:) * normal_vec(:) )
691 stress_t_tmpvec(:) = stress_t_tmpvec(:) - stress_t_tmp * normal_vec(:)
692
693 ! T1j
694 t11(ip) = t11(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(1)
695 t12(ip) = t12(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(2)
696 t13(ip) = t13(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(3)
697
698 ! T2j
699 t21(ip) = t21(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(1)
700 t22(ip) = t22(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(2)
701 t23(ip) = t23(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(3)
702
703 ! T3j
704 t31(ip) = t31(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(1)
705 t32(ip) = t32(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(2)
706 t33(ip) = t33(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(3)
707
708 case ( bnd_type_noslip_id )
709 end select
710
711 select case( this%ThermalBC_list(domid)%list(i_) )
712 case ( bnd_type_adiabat_id )
713 normal_flux = df1(im) * nx(i) + df2(im) * ny(i) &
714 + ( df3(im) / gsqrtv + g13(im) * df1(im) + g23(im) * df2(im) ) * nz(i)
715 df1(ip) = df1(im) - 2.0_rp * normal_flux * normal_vec(1)
716 df2(ip) = df2(im) - 2.0_rp * normal_flux * normal_vec(2)
717 df3(ip) = df3(im) - 2.0_rp * normal_flux * normal_vec(3)
718 end select
719 end if
720 end do
721 end do
722
723 return
724 end subroutine atmos_dyn_bnd_applybc_tbstress_lc
725
726!OCL SERIAL
727 subroutine atmos_dyn_bnd_inquire_bound_flag( this, & ! (in)
728 is_bound, & ! (out)
729 domid, vmapm, vmapp, vmapb, lmesh, elem ) ! (in)
730
731 use scale_mesh_bndinfo, only: &
732 bnd_type_slip_id, bnd_type_noslip_id, &
733 bnd_type_adiabat_id
734 implicit none
735
736 class(atmdynbnd), intent(in) :: this
737 class(localmesh3d), intent(in) :: lmesh
738 class(elementbase3d), intent(in) :: elem
739 logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
740 integer, intent(in) :: domid
741 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
742 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
743 integer, intent(in) :: vmapb(:)
744
745 integer :: i, i_, im, ip
746 !-----------------------------------------------
747
748 do i=1, elem%NfpTot*lmesh%Ne
749 ip = vmapp(i)
750 i_ = ip - elem%Np*lmesh%NeE
751 is_bound(i) = .false.
752
753 if (i_ > 0) then
754 im = vmapm(i)
755 if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
756 is_bound(i) = .true.
757 else if ( this%VelBC_list(domid)%list(i_) == bnd_type_noslip_id ) then
758 is_bound(i) = .true.
759 end if
760
761 end if
762 end do
763
764 return
765 end subroutine atmos_dyn_bnd_inquire_bound_flag
766
767 !-- private -------------------------------------------------------
768
769!OCL SERIAL
770 subroutine bnd_init_lc( &
771 velBCInfo, thermalBCInfo, & ! (inout)
772 velbc_ids, thermalbc_ids, & ! (in)
773 thermal_fixval, & ! (in)
774 vmapb, mesh, lmesh, elem ) ! (in)
775
776 use scale_mesh_bndinfo, only: bnd_type_nospec_id
777 implicit none
778
779 type(meshbndinfo), intent(inout) :: velbcinfo
780 type(meshbndinfo), intent(inout) :: thermalbcinfo
781 integer, intent(in) :: velbc_ids(dom_bnd_num)
782 integer, intent(in) :: thermalbc_ids(dom_bnd_num)
783 real(rp), intent(in) :: thermal_fixval(dom_bnd_num)
784 class(meshbase), intent(in) :: mesh
785 class(localmesh3d), intent(in) :: lmesh
786 class(elementbase3d), intent(in) :: elem
787 integer, intent(in) :: vmapb(:)
788
789 integer :: tileid
790 integer :: dom_bnd_sizes(dom_bnd_num)
791 integer :: bnd_buf_size
792 integer :: b, is_, ie_
793
794 !-----------------------------------------------
795
796 dom_bnd_sizes(:) = &
797 elem%Nfp_h*lmesh%NeZ*(/ lmesh%NeX, lmesh%NeY, lmesh%NeX, lmesh%NeY, 0, 0 /) &
798 + elem%Nfp_v*lmesh%NeX*lmesh%NeY*(/ 0, 0, 0, 0, 1, 1 /)
799 bnd_buf_size = sum(dom_bnd_sizes)
800
801 call velbcinfo%Init( bnd_buf_size )
802 call velbcinfo%Set(1, bnd_buf_size, bnd_type_nospec_id)
803
804 call thermalbcinfo%Init( bnd_buf_size )
805 call thermalbcinfo%Set(1, bnd_buf_size, bnd_type_nospec_id)
806
807 tileid = lmesh%tileID
808 is_ = 1
809 do b=1, dom_bnd_num
810 ie_ = is_ + dom_bnd_sizes(b) - 1
811 if ( mesh%tileID_globalMap(b,tileid) == tileid &
812 .and. mesh%tileFaceID_globalMap(b,tileid) == b ) then
813 call velbcinfo%Set( is_, ie_, velbc_ids(b) )
814 call thermalbcinfo%Set( is_, ie_, thermalbc_ids(b), thermal_fixval(b) )
815 end if
816 is_ = ie_ + 1
817 end do
818
819 return
820 end subroutine bnd_init_lc
821
822end module scale_atm_dyn_dgm_bnd
scale_atm_dyn_dgm_bnd
module FElib / Fluid dyn solver / Atmosphere / Boundary
Definition scale_atm_dyn_dgm_bnd.F90:11
scale_atm_dyn_dgm_bnd::dombnd_south_id
integer, parameter dombnd_south_id
Definition scale_atm_dyn_dgm_bnd.F90:84
scale_atm_dyn_dgm_nonhydro3d_common
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / Common
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:11
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_momy_id
integer, parameter, public prgvar_momy_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:74
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_ddens_id
integer, parameter, public prgvar_ddens_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:68
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_momz_id
integer, parameter, public prgvar_momz_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:72
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_momx_id
integer, parameter, public prgvar_momx_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:73
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_drhot_id
integer, parameter, public prgvar_drhot_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:70
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_num
integer, parameter, public prgvar_num
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:77
scale_element_base
module FElib / Element / Base
Definition scale_element_base.F90:11
scale_localmesh_2d
module FElib / Mesh / Local 2D
Definition scale_localmesh_2d.F90:10
scale_localmesh_3d
module FElib / Mesh / Local 3D
Definition scale_localmesh_3d.F90:10
scale_localmesh_base
module FElib / Mesh / Local, Base
Definition scale_localmesh_base.F90:10
scale_localmeshfield_base
module FElib / Data / base
Definition scale_localmeshfield_base.F90:10
scale_mesh_base2d
module FElib / Mesh / Base 2D
Definition scale_mesh_base2d.F90:10
scale_mesh_base3d
module FElib / Mesh / Base 3D
Definition scale_mesh_base3d.F90:10
scale_mesh_base
module FElib / Mesh / Base
Definition scale_mesh_base.F90:10
scale_mesh_bndinfo
module FElib / Mesh / Boundary information
Definition scale_mesh_bndinfo.F90:11
scale_mesh_bndinfo::bnd_type_slip_id
integer, parameter, public bnd_type_slip_id
Definition scale_mesh_bndinfo.F90:51
scale_mesh_bndinfo::bnd_type_nospec_name
character(len= *), parameter, public bnd_type_nospec_name
Definition scale_mesh_bndinfo.F90:46
scale_mesh_bndinfo::bnd_type_fixval_id
integer, parameter, public bnd_type_fixval_id
Definition scale_mesh_bndinfo.F90:57
scale_mesh_bndinfo::bndtype_nametoid
integer function, public bndtype_nametoid(bnd_type_name)
Definition scale_mesh_bndinfo.F90:145
scale_mesh_bndinfo::bnd_type_periodic_id
integer, parameter, public bnd_type_periodic_id
Definition scale_mesh_bndinfo.F90:49
scale_mesh_bndinfo::bnd_type_nospec_id
integer, parameter, public bnd_type_nospec_id
Definition scale_mesh_bndinfo.F90:47
scale_mesh_bndinfo::bnd_type_adiabat_id
integer, parameter, public bnd_type_adiabat_id
Definition scale_mesh_bndinfo.F90:55
scale_mesh_bndinfo::bnd_type_noslip_id
integer, parameter, public bnd_type_noslip_id
Definition scale_mesh_bndinfo.F90:53
scale_meshfield_base
module FElib / Data / base
Definition scale_meshfield_base.F90:10
scale_atm_dyn_dgm_bnd::atmdynbnd
A derived type useful for apply boundary conditions.
Definition scale_atm_dyn_dgm_bnd.F90:57
scale_element_base::elementbase2d
Derived type representing a 2D reference element.
Definition scale_element_base.F90:71
scale_element_base::elementbase3d
Derived type representing a 3D reference element.
Definition scale_element_base.F90:107
scale_element_base::elementbase
Derived type representing an arbitrary finite element.
Definition scale_element_base.F90:30
scale_localmesh_2d::localmesh2d
Definition scale_localmesh_2d.F90:29
scale_localmesh_3d::localmesh3d
Derived type to manage a local 3D computational domain.
Definition scale_localmesh_3d.F90:32
scale_localmesh_base::localmeshbase
Derived type to manage a local computational domain (base type)
Definition scale_localmesh_base.F90:31
scale_localmeshfield_base::localmeshfield3d
Derived type representing a field with 3D local mesh.
Definition scale_localmeshfield_base.F90:69
scale_mesh_base2d::meshbase2d
Definition scale_mesh_base2d.F90:34
scale_mesh_base3d::meshbase3d
Definition scale_mesh_base3d.F90:36
scale_mesh_base::meshbase
Definition scale_mesh_base.F90:38
scale_mesh_bndinfo::meshbndinfo
Definition scale_mesh_bndinfo.F90:28
scale_meshfield_base::meshfield3d
Derived type representing a field with 3D mesh.
Definition scale_meshfield_base.F90:95
Generated by doxygen 1.14.0