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 !$acc enter data create( this )
169
170 allocate( this%velBC_ids(dom_bnd_num) )
171 this%velBC_ids(dombnd_btm_id) = bndtype_nametoid(btm_vel_bc)
172 this%velBC_ids(dombnd_top_id) = bndtype_nametoid(top_vel_bc)
173 this%velBC_ids(dombnd_north_id) = bndtype_nametoid(north_vel_bc)
174 this%velBC_ids(dombnd_south_id) = bndtype_nametoid(south_vel_bc)
175 this%velBC_ids(dombnd_east_id) = bndtype_nametoid(east_vel_bc)
176 this%velBC_ids(dombnd_west_id) = bndtype_nametoid(west_vel_bc)
177 !$acc enter data copyin( this%velBC_ids )
178
179 !--
180
181 allocate( this%thermalBC_ids(dom_bnd_num) )
182 this%thermalBC_ids(dombnd_btm_id) = bndtype_nametoid(btm_thermal_bc)
183 this%thermalBC_ids(dombnd_top_id) = bndtype_nametoid(top_thermal_bc)
184 this%thermalBC_ids(dombnd_north_id) = bndtype_nametoid(north_thermal_bc)
185 this%thermalBC_ids(dombnd_south_id) = bndtype_nametoid(south_thermal_bc)
186 this%thermalBC_ids(dombnd_east_id) = bndtype_nametoid(east_thermal_bc)
187 this%thermalBC_ids(dombnd_west_id) = bndtype_nametoid(west_thermal_bc)
188 !$acc enter data copyin( this%thermalBC_ids )
189
190 allocate( this%thermal_fixval(dom_bnd_num) )
191 this%thermal_fixval(dombnd_btm_id) = btm_thermal_fixval
192 this%thermal_fixval(dombnd_top_id) = top_thermal_fixval
193 this%thermal_fixval(dombnd_north_id) = north_thermal_fixval
194 this%thermal_fixval(dombnd_south_id) = south_thermal_fixval
195 this%thermal_fixval(dombnd_east_id) = east_thermal_fixval
196 this%thermal_fixval(dombnd_west_id) = west_thermal_fixval
197 !$acc enter data copyin( this%thermal_fixval )
198
199 !------
200
201 return
202 end subroutine atmos_dyn_bnd_setup
203
204!> Finalize an object to manage boundary conditions with dynamical process
205!OCL SERIAL
206 subroutine atmos_dyn_bnd_finalize( this )
207 implicit none
208 class(atmdynbnd), intent(inout) :: this
209
210 integer :: n
211 !--------------------------------------
212
213 if ( allocated(this%VelBC_list) ) then
214 do n=1, size(this%VelBC_list)
215 call this%VelBC_list(n)%Final()
216 call this%ThermalBC_list(n)%Final()
217 end do
218 deallocate( this%VelBC_list, this%ThermalBC_list )
219 !$acc exit data delete( this%VelBC_list, this%ThermalBC_list )
220
221 deallocate( this%velBC_ids, this%thermalBC_ids )
222 !$acc exit data delete( this%velBC_ids, this%thermalBC_ids )
223 end if
224 !$acc exit data delete( this )
225 return
226 end subroutine atmos_dyn_bnd_finalize
227
228!OCL SERIAL
229 subroutine atmos_dyn_bnd_setbcinfo( this, mesh )
230
231 implicit none
232
233 class(atmdynbnd), intent(inout) :: this
234 class(meshbase), target, intent(in) :: mesh
235
236 integer :: n
237 integer :: b
238 class(localmeshbase), pointer :: ptr_lcmesh
239 class(localmesh3d), pointer :: lcmesh3d
240 !--------------------------------------------------
241
242
243 allocate( this%VelBC_list(mesh%LOCAL_MESH_NUM) )
244 allocate( this%ThermalBC_list(mesh%LOCAL_MESH_NUM) )
245 !$acc enter data create( this%VelBC_list, this%ThermalBC_list )
246
247 nullify( lcmesh3d )
248
249 do n=1, mesh%LOCAL_MESH_NUM
250 call mesh%GetLocalMesh( n, ptr_lcmesh )
251 select type (ptr_lcmesh)
252 type is (localmesh3d)
253 lcmesh3d => ptr_lcmesh
254 end select
255
256 call bnd_init_lc( &
257 this%VelBC_list(n), this%ThermalBC_list(n), & ! (inout)
258 this%velBC_ids(:), this%thermalBC_ids(:), & ! (in)
259 this%thermal_fixval(:), & ! (in)
260 lcmesh3d%VMapB, mesh, lcmesh3d, lcmesh3d%refElem3D ) ! (in)
261
262 end do
263
264 return
265 end subroutine atmos_dyn_bnd_setbcinfo
266
267!> Set exterior values at element boundaries to evaluate inviscid numerical fluxes
268!!
269!OCL SERIAL
270 subroutine atmos_dyn_bnd_applybc_prgvars_lc( this, &
271 domID, & ! (in)
272 ddens, momx, momy, momz, therm, & ! (inout)
273 dens_hyd, pres_hyd, & ! (in)
274 gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
275 vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d ) ! (in)
276
277 use scale_mesh_bndinfo, only: &
278 bnd_type_slip_id, bnd_type_noslip_id, &
279 bnd_type_adiabat_id
280 use scale_prc
281 implicit none
282
283 class(atmdynbnd), intent(in) :: this
284 integer, intent(in) :: domid
285 class(localmesh3d), intent(in) :: lmesh
286 class(elementbase3d), intent(in) :: elem
287 class(localmesh2d), intent(in) :: lmesh2d
288 class(elementbase2d), intent(in) :: elem2d
289 real(rp), intent(inout) :: ddens(elem%np*lmesh%nea)
290 real(rp), intent(inout) :: momx(elem%np*lmesh%nea)
291 real(rp), intent(inout) :: momy(elem%np*lmesh%nea)
292 real(rp), intent(inout) :: momz(elem%np*lmesh%nea)
293 real(rp), intent(inout) :: therm(elem%np*lmesh%nea)
294 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
295 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
296 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
297 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
298 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
299 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
300 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
301 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
302 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
303 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
304 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
305 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
306 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
307 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
308 integer, intent(in) :: vmapb(:)
309
310 integer :: p, ke, ke2d
311 integer :: i, i_, im, ip
312 real(rp) :: mom_normal
313
314 real(rp) :: momw
315 real(rp) :: gsqrtv, g11_, g12_, g22_
316 real(rp) :: fac
317
318 integer :: indexh2dto3d_bnd(elem%nfptot)
319 !-----------------------------------------------
320
321 indexh2dto3d_bnd(:) = elem%IndexH2Dto3D_bnd(:)
322
323 !$omp parallel do collapse(2) private( &
324 !$omp ke, p, ke2D, i, i_, iM, iP, &
325 !$omp mom_normal, MOMW, GsqrtV, G11_, G12_, G22_, fac )
326 !$acc parallel loop gang &
327 !$acc present( DDENS, MOMX, MOMY, MOMZ, THERM, DENS_hyd, PRES_hyd, &
328 !$acc Gsqrt, GsqrtH, G11, G12, G22, G13, G23, nx, ny, nz, vmapM, vmapP, vmapB, lmesh,elem )
329 do ke=lmesh%NeS, lmesh%NeE
330 !$acc loop vector
331 do p=1, elem%NfpTot
332 i = p + (ke-1)*elem%NfpTot
333 ip = vmapp(i)
334 i_ = ip - elem%Np*lmesh%NeE
335
336 if (i_ > 0) then
337 im = vmapm(i)
338
339 select case( this%VelBC_list(domid)%list(i_) )
340 case ( bnd_type_slip_id)
341 ke2d = lmesh%EMap3Dto2D(ke)
342
343 gsqrtv = gsqrt(im) / gsqrth(indexh2dto3d_bnd(p),ke2d)
344 g11_ = g11(indexh2dto3d_bnd(p),ke2d)
345 g12_ = g12(indexh2dto3d_bnd(p),ke2d)
346 g22_ = g22(indexh2dto3d_bnd(p),ke2d)
347
348 momw = momz(im) / gsqrtv &
349 + g13(im) * momx(im) + g23(im) * momy(im)
350 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 )
351
352 mom_normal = momx(im) * nx(i) + momy(im) * ny(i) + momw * nz(i)
353 momx(ip) = momx(im) - 2.0_rp * mom_normal * ( nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) ) )
354 momy(ip) = momy(im) - 2.0_rp * mom_normal * ( ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) ) )
355 momz(ip) = momz(im) - 2.0_rp * mom_normal * fac / gsqrtv
356 case ( bnd_type_noslip_id )
357 momx(ip) = - momx(im)
358 momy(ip) = - momy(im)
359 momz(ip) = - momz(im)
360 end select
361 end if
362
363 end do
364 end do
365
366 return
367 end subroutine atmos_dyn_bnd_applybc_prgvars_lc
368
369!OCL SERIAL
370 subroutine atmos_dyn_bnd_applybc_numdiff_odd_lc( this, & ! (in)
371 gxvar, gyvar, gzvar, & ! (inout)
372 is_bound, & ! (out)
373 varid, domid, & ! (in)
374 nx, ny, nz, vmapm, vmapp, vmapb, lmesh, elem ) ! (in)
375
376 use scale_mesh_bndinfo, only: &
377 bnd_type_slip_id, bnd_type_noslip_id, &
378 bnd_type_adiabat_id
379
380 implicit none
381
382 class(atmdynbnd), intent(in) :: this
383 class(localmesh3d), intent(in) :: lmesh
384 class(elementbase3d), intent(in) :: elem
385 real(rp), intent(inout) :: gxvar(elem%np*lmesh%nea)
386 real(rp), intent(inout) :: gyvar(elem%np*lmesh%nea)
387 real(rp), intent(inout) :: gzvar(elem%np*lmesh%nea)
388 logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
389 integer, intent(in) :: varid
390 integer, intent(in) :: domid
391 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
392 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
393 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
394 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
395 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
396 integer, intent(in) :: vmapb(:)
397
398 integer :: i, i_, im, ip
399 real(rp) :: grad_normal
400 !-----------------------------------------------
401
402 do i=1, elem%NfpTot*lmesh%Ne
403 ip = vmapp(i)
404 i_ = ip - elem%Np*lmesh%NeE
405 is_bound(i) = .false.
406
407 if (i_ > 0) then
408 im = vmapm(i)
409 grad_normal = gxvar(im) * nx(i) + gyvar(im) * ny(i) + gzvar(im) * nz(i)
410
411 if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
412 select case(varid)
413 case(momx_vid)
414 gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
415 gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
416 case(momy_vid)
417 gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
418 gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
419 case(momz_vid)
420 gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
421 gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
422 end select
423 is_bound(i) = .true.
424
425 end if
426 if ( this%ThermalBC_list(domid)%list(i_) == bnd_type_adiabat_id ) then
427 select case(varid)
428 case(dens_vid, rhot_vid)
429 gxvar(ip) = gxvar(im) - 2.0_rp * grad_normal * nx(i)
430 gyvar(ip) = gyvar(im) - 2.0_rp * grad_normal * ny(i)
431 gzvar(ip) = gzvar(im) - 2.0_rp * grad_normal * nz(i)
432 end select
433 is_bound(i) = .true.
434
435 end if
436
437 end if
438
439 end do
440
441 return
442 end subroutine atmos_dyn_bnd_applybc_numdiff_odd_lc
443
444!OCL SERIAL
445 subroutine atmos_dyn_bnd_applybc_numdiff_even_lc( this, & ! (in)
446 var, & ! (inout)
447 is_bound, & ! (out)
448 varid, domid, & ! (in)
449 nx, ny, nz, vmapm, vmapp, vmapb, lmesh, elem ) ! (in)
450
451 use scale_mesh_bndinfo, only: &
452 bnd_type_slip_id, bnd_type_noslip_id, &
453 bnd_type_adiabat_id
454 implicit none
455
456 class(atmdynbnd), intent(in) :: this
457 class(localmesh3d), intent(in) :: lmesh
458 class(elementbase3d), intent(in) :: elem
459 real(rp), intent(inout) :: var(elem%np*lmesh%nea)
460 logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
461 integer, intent(in) :: varid
462 integer, intent(in) :: domid
463 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
464 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
465 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
466 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
467 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
468 integer, intent(in) :: vmapb(:)
469
470 integer :: i, i_, im, ip
471 real(rp) :: grad_normal
472 !-----------------------------------------------
473
474 do i=1, elem%NfpTot*lmesh%Ne
475 ip = vmapp(i)
476 i_ = ip - elem%Np*lmesh%NeE
477 is_bound(i) = .false.
478
479 if (i_ > 0) then
480 im = vmapm(i)
481
482 if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
483 select case(varid)
484 case(momx_vid)
485 grad_normal = var(im) * nx(i)
486 var(ip) = var(im) - 2.0_rp * grad_normal * nx(i)
487 case(momy_vid)
488 grad_normal = var(im) * ny(i)
489 var(ip) = var(im) - 2.0_rp * grad_normal * ny(i)
490 case(momz_vid)
491 grad_normal = var(im) * nz(i)
492 var(ip) = var(im) - 2.0_rp * grad_normal * nz(i)
493 end select
494 is_bound(i) = .true.
495
496 else if ( this%VelBC_list(domid)%list(i_) == bnd_type_noslip_id ) then
497 select case(varid)
498 case( momx_vid, momy_vid, momz_vid )
499 var(ip) = - var(im)
500 end select
501 is_bound(i) = .true.
502 end if
503
504 end if
505 end do
506
507 return
508 end subroutine atmos_dyn_bnd_applybc_numdiff_even_lc
509
510
511!> Set exterior values at boundaries for the turbulent schemes
512!!
513!OCL SERIAL
514 subroutine atmos_dyn_bnd_applybc_tbvars_lc( this, &
515 domID, & ! (in)
516 ddens, momx, momy, momz, pt, pres, & ! (inout)
517 dens_hyd, pres_hyd, rtot, cptot, & ! (in)
518 gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
519 vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d ) ! (in)
520
521 use scale_const, only: &
522 pres00 => const_pre00
523 use scale_mesh_bndinfo, only: &
524 bnd_type_slip_id, bnd_type_noslip_id, &
525 bnd_type_adiabat_id, bnd_type_fixval_id
526
527 implicit none
528
529 class(atmdynbnd), intent(in) :: this
530 integer, intent(in) :: domid
531 class(localmesh3d), intent(in) :: lmesh
532 class(elementbase3d), intent(in) :: elem
533 class(localmesh2d), intent(in) :: lmesh2d
534 class(elementbase2d), intent(in) :: elem2d
535 real(rp), intent(inout) :: ddens(elem%np*lmesh%nea)
536 real(rp), intent(inout) :: momx(elem%np*lmesh%nea)
537 real(rp), intent(inout) :: momy(elem%np*lmesh%nea)
538 real(rp), intent(inout) :: momz(elem%np*lmesh%nea)
539 real(rp), intent(inout) :: pt(elem%np*lmesh%nea)
540 real(rp), intent(inout) :: pres(elem%np*lmesh%nea)
541 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
542 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
543 real(rp), intent(in) :: rtot(elem%np*lmesh%nea)
544 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
545 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
546 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
547 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
548 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
549 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
550 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
551 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
552 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
553 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
554 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
555 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
556 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
557 integer, intent(in) :: vmapb(:)
558
559 integer :: p, ke, ke2d
560 integer :: i, i_, im, ip
561 real(rp) :: mom_normal
562
563 real(rp) :: momw
564 real(rp) :: temp_p, temp_b
565 real(rp) :: gsqrtv, g11_, g12_, g22_
566 real(rp) :: fac
567 !-----------------------------------------------
568
569 !$omp parallel do collapse(2) private( &
570 !$omp ke, p, ke2D, i, i_, iM, iP, &
571 !$omp mom_normal, MOMW, GsqrtV, G11_, G12_, G22_, fac, &
572 !$omp TEMP_P, TEMP_B )
573 do ke=lmesh%NeS, lmesh%NeE
574 do p=1, elem%NfpTot
575 i = p + (ke-1)*elem%NfpTot
576 ip = vmapp(i)
577 i_ = ip - elem%Np*lmesh%NeE
578
579 if (i_ > 0) then
580 im = vmapm(i)
581
582 select case( this%VelBC_list(domid)%list(i_) )
583 case ( bnd_type_slip_id)
584 ke2d = lmesh%EMap3Dto2D(ke)
585
586 gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d)
587 g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
588 g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
589 g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
590
591 momw = momz(im) / gsqrtv &
592 + g13(im) * momx(im) + g23(im) * momy(im)
593 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 )
594
595 mom_normal = momx(im) * nx(i) + momy(im) * ny(i) + momw * nz(i)
596 momx(ip) = momx(im) - 2.0_rp * mom_normal * ( nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) ) )
597 momy(ip) = momy(im) - 2.0_rp * mom_normal * ( ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) ) )
598 momz(ip) = momz(im) - 2.0_rp * mom_normal * fac / gsqrtv
599
600 case ( bnd_type_noslip_id )
601 momx(ip) = - momx(im)
602 momy(ip) = - momy(im)
603 momz(ip) = - momz(im)
604 end select
605
606 select case( this%ThermalBC_list(domid)%list(i_) )
607 case ( bnd_type_fixval_id )
608 temp_b = this%ThermalBC_list(domid)%val(i_)
609 temp_p = 2.0_rp * temp_b - pres(im) / ( ( dens_hyd(im) + ddens(im) ) * rtot(im) )
610 pt(ip) = temp_p * ( pres00 / pres(im) )**(rtot(im)/cptot(im))
611 ddens(ip) = pres(im) / ( rtot(im) * temp_p ) - dens_hyd(im)
612 end select
613 end if
614
615 end do
616 end do
617
618 return
619 end subroutine atmos_dyn_bnd_applybc_tbvars_lc
620
621!> Set exterior shear-stress tensor and heat flux at boundaries for the turbulent schemes
622!!
623!! Note: This codes are tentatively implementated.
624!! We need to check whether the formulations is valid for the case when general vertical coordinate is introduced.
625!OCL SERIAL
626 subroutine atmos_dyn_bnd_applybc_tbstress_lc( this, &
627 domID, & ! (in)
628 t11, t12, t13, t21, t22, t23, t31, t32, t33, & ! (inout)
629 df1, df2, df3, & ! (in)
630 gsqrt, gsqrth, g11, g12, g22, g13, g23, nx, ny, nz, & ! (in)
631 vmapm, vmapp, vmapb, lmesh, elem, lmesh2d, elem2d ) ! (in)
632
633 use scale_mesh_bndinfo, only: &
634 bnd_type_slip_id, bnd_type_noslip_id, &
635 bnd_type_adiabat_id, bnd_type_fixval_id
636
637 implicit none
638
639 class(atmdynbnd), intent(in) :: this
640 integer, intent(in) :: domid
641 class(localmesh3d), intent(in) :: lmesh
642 class(elementbase3d), intent(in) :: elem
643 class(localmesh2d), intent(in) :: lmesh2d
644 class(elementbase2d), intent(in) :: elem2d
645 real(rp), intent(inout) :: t11(elem%np*lmesh%nea), t12(elem%np*lmesh%nea), t13(elem%np*lmesh%nea)
646 real(rp), intent(inout) :: t21(elem%np*lmesh%nea), t22(elem%np*lmesh%nea), t23(elem%np*lmesh%nea)
647 real(rp), intent(inout) :: t31(elem%np*lmesh%nea), t32(elem%np*lmesh%nea), t33(elem%np*lmesh%nea)
648 real(rp), intent(inout) :: df1(elem%np*lmesh%nea), df2(elem%np*lmesh%nea), df3(elem%np*lmesh%nea)
649 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
650 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
651 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
652 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
653 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
654 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
655 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
656 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
657 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
658 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
659 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
660 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
661 integer, intent(in) :: vmapb(:)
662
663 integer :: p, ke, ke2d
664 integer :: i, i_, im, ip
665
666 real(rp) :: gsqrtv, g11_, g12_, g22_
667 real(rp) :: fac, normal_vec(3)
668
669 real(rp) :: stress_t_tmp, stress_t_tmpvec(3)
670 real(rp) :: vecw
671 real(rp) :: normal_flux
672 !-----------------------------------------------
673
674 !$omp parallel do collapse(2) private( &
675 !$omp ke, p, ke2D, i, i_, iM, iP, &
676 !$omp VECW, GsqrtV, G11_, G12_, G22_, fac, normal_vec, &
677 !$omp stress_t_tmp, stress_t_tmpvec, normal_flux )
678 do ke=lmesh%NeS, lmesh%NeE
679 do p=1, elem%NfpTot
680 i = p + (ke-1)*elem%NfpTot
681 ip = vmapp(i)
682 i_ = ip - elem%Np*lmesh%NeE
683
684 if (i_ > 0) then
685 ke2d = lmesh%EMap3Dto2D(ke)
686 im = vmapm(i)
687
688 gsqrtv = gsqrt(im) / gsqrth(elem%IndexH2Dto3D_bnd(p),ke2d)
689 g11_ = g11(elem%IndexH2Dto3D_bnd(p),ke2d)
690 g12_ = g12(elem%IndexH2Dto3D_bnd(p),ke2d)
691 g22_ = g22(elem%IndexH2Dto3D_bnd(p),ke2d)
692 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 )
693 normal_vec(1) = nx(i) + fac * ( g11_ * g13(im) + g12_ * g23(im) )
694 normal_vec(2) = ny(i) + fac * ( g12_ * g13(im) + g22_ * g23(im) )
695 normal_vec(3) = fac / gsqrtv
696
697 select case( this%VelBC_list(domid)%list(i_) )
698 case ( bnd_type_slip_id)
699 vecw = t13(im) / gsqrtv + g13(im) * t11(im) + g23(im) * t12(im)
700 stress_t_tmpvec(1) = t11(im) * nx(i) + t12(im) * ny(i) + vecw * nz(i)
701
702 vecw = t23(im) / gsqrtv + g13(im) * t21(im) + g23(im) * t22(im)
703 stress_t_tmpvec(2) = t21(im) * nx(i) + t22(im) * ny(i) + vecw * nz(i)
704
705 vecw = t33(im) / gsqrtv
706 stress_t_tmpvec(3) = t31(im) * nx(i) + t32(im) * ny(i) + vecw * nz(i)
707
708 stress_t_tmp = sum( stress_t_tmpvec(:) * normal_vec(:) )
709 stress_t_tmpvec(:) = stress_t_tmpvec(:) - stress_t_tmp * normal_vec(:)
710
711 ! T1j
712 t11(ip) = t11(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(1)
713 t12(ip) = t12(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(2)
714 t13(ip) = t13(im) - 2.0_rp * stress_t_tmpvec(1) * normal_vec(3)
715
716 ! T2j
717 t21(ip) = t21(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(1)
718 t22(ip) = t22(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(2)
719 t23(ip) = t23(im) - 2.0_rp * stress_t_tmpvec(2) * normal_vec(3)
720
721 ! T3j
722 t31(ip) = t31(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(1)
723 t32(ip) = t32(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(2)
724 t33(ip) = t33(im) - 2.0_rp * stress_t_tmpvec(3) * normal_vec(3)
725
726 case ( bnd_type_noslip_id )
727 end select
728
729 select case( this%ThermalBC_list(domid)%list(i_) )
730 case ( bnd_type_adiabat_id )
731 normal_flux = df1(im) * nx(i) + df2(im) * ny(i) &
732 + ( df3(im) / gsqrtv + g13(im) * df1(im) + g23(im) * df2(im) ) * nz(i)
733 df1(ip) = df1(im) - 2.0_rp * normal_flux * normal_vec(1)
734 df2(ip) = df2(im) - 2.0_rp * normal_flux * normal_vec(2)
735 df3(ip) = df3(im) - 2.0_rp * normal_flux * normal_vec(3)
736 end select
737 end if
738 end do
739 end do
740
741 return
742 end subroutine atmos_dyn_bnd_applybc_tbstress_lc
743
744!OCL SERIAL
745 subroutine atmos_dyn_bnd_inquire_bound_flag( this, & ! (in)
746 is_bound, & ! (out)
747 domid, vmapm, vmapp, vmapb, lmesh, elem ) ! (in)
748
749 use scale_mesh_bndinfo, only: &
750 bnd_type_slip_id, bnd_type_noslip_id, &
751 bnd_type_adiabat_id
752 implicit none
753
754 class(atmdynbnd), intent(in) :: this
755 class(localmesh3d), intent(in) :: lmesh
756 class(elementbase3d), intent(in) :: elem
757 logical, intent(out) :: is_bound(elem%nfptot*lmesh%ne)
758 integer, intent(in) :: domid
759 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
760 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
761 integer, intent(in) :: vmapb(:)
762
763 integer :: i, i_, im, ip
764 !-----------------------------------------------
765
766 do i=1, elem%NfpTot*lmesh%Ne
767 ip = vmapp(i)
768 i_ = ip - elem%Np*lmesh%NeE
769 is_bound(i) = .false.
770
771 if (i_ > 0) then
772 im = vmapm(i)
773 if ( this%VelBC_list(domid)%list(i_) == bnd_type_slip_id ) then
774 is_bound(i) = .true.
775 else if ( this%VelBC_list(domid)%list(i_) == bnd_type_noslip_id ) then
776 is_bound(i) = .true.
777 end if
778
779 end if
780 end do
781
782 return
783 end subroutine atmos_dyn_bnd_inquire_bound_flag
784
785 !-- private -------------------------------------------------------
786
787!OCL SERIAL
788 subroutine bnd_init_lc( &
789 velBCInfo, thermalBCInfo, & ! (inout)
790 velbc_ids, thermalbc_ids, & ! (in)
791 thermal_fixval, & ! (in)
792 vmapb, mesh, lmesh, elem ) ! (in)
793
794 use scale_mesh_bndinfo, only: bnd_type_nospec_id
795 implicit none
796
797 type(meshbndinfo), intent(inout) :: velbcinfo
798 type(meshbndinfo), intent(inout) :: thermalbcinfo
799 integer, intent(in) :: velbc_ids(dom_bnd_num)
800 integer, intent(in) :: thermalbc_ids(dom_bnd_num)
801 real(rp), intent(in) :: thermal_fixval(dom_bnd_num)
802 class(meshbase), intent(in) :: mesh
803 class(localmesh3d), intent(in) :: lmesh
804 class(elementbase3d), intent(in) :: elem
805 integer, intent(in) :: vmapb(:)
806
807 integer :: tileid
808 integer :: dom_bnd_sizes(dom_bnd_num)
809 integer :: bnd_buf_size
810 integer :: b, is_, ie_
811
812 !-----------------------------------------------
813
814 dom_bnd_sizes(:) = &
815 elem%Nfp_h*lmesh%NeZ*(/ lmesh%NeX, lmesh%NeY, lmesh%NeX, lmesh%NeY, 0, 0 /) &
816 + elem%Nfp_v*lmesh%NeX*lmesh%NeY*(/ 0, 0, 0, 0, 1, 1 /)
817 bnd_buf_size = sum(dom_bnd_sizes)
818
819 call velbcinfo%Init( bnd_buf_size )
820 call velbcinfo%Set(1, bnd_buf_size, bnd_type_nospec_id)
821
822 call thermalbcinfo%Init( bnd_buf_size )
823 call thermalbcinfo%Set(1, bnd_buf_size, bnd_type_nospec_id)
824
825 tileid = lmesh%tileID
826 is_ = 1
827 do b=1, dom_bnd_num
828 ie_ = is_ + dom_bnd_sizes(b) - 1
829 if ( mesh%tileID_globalMap(b,tileid) == tileid &
830 .and. mesh%tileFaceID_globalMap(b,tileid) == b ) then
831 call velbcinfo%Set( is_, ie_, velbc_ids(b) )
832 call thermalbcinfo%Set( is_, ie_, thermalbc_ids(b), thermal_fixval(b) )
833 end if
834 is_ = ie_ + 1
835 end do
836
837 return
838 end subroutine bnd_init_lc
839
840end 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:52
scale_mesh_bndinfo::bnd_type_nospec_name
character(len= *), parameter, public bnd_type_nospec_name
Definition scale_mesh_bndinfo.F90:47
scale_mesh_bndinfo::bnd_type_fixval_id
integer, parameter, public bnd_type_fixval_id
Definition scale_mesh_bndinfo.F90:58
scale_mesh_bndinfo::bndtype_nametoid
integer function, public bndtype_nametoid(bnd_type_name)
Definition scale_mesh_bndinfo.F90:156
scale_mesh_bndinfo::bnd_type_periodic_id
integer, parameter, public bnd_type_periodic_id
Definition scale_mesh_bndinfo.F90:50
scale_mesh_bndinfo::bnd_type_nospec_id
integer, parameter, public bnd_type_nospec_id
Definition scale_mesh_bndinfo.F90:48
scale_mesh_bndinfo::bnd_type_adiabat_id
integer, parameter, public bnd_type_adiabat_id
Definition scale_mesh_bndinfo.F90:56
scale_mesh_bndinfo::bnd_type_noslip_id
integer, parameter, public bnd_type_noslip_id
Definition scale_mesh_bndinfo.F90:54
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:111
scale_element_base::elementbase
Derived type representing an arbitrary finite element.
Definition scale_element_base.F90:30
scale_localmesh_2d::localmesh2d
Derived type representing a local mesh for 2D domain.
Definition scale_localmesh_2d.F90:30
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
Derived type to manage a computational mesh (base type for 2D domain)
Definition scale_mesh_base2d.F90:35
scale_mesh_base3d::meshbase3d
Derived type to manage a computational mesh (base type for 3D domain)
Definition scale_mesh_base3d.F90:37
scale_mesh_base::meshbase
Base type to manage a computational mesh.
Definition scale_mesh_base.F90:40
scale_mesh_bndinfo::meshbndinfo
Derived type to manage boundary information of computational domain.
Definition scale_mesh_bndinfo.F90:29
scale_meshfield_base::meshfield3d
Derived type representing a field with 3D mesh.
Definition scale_meshfield_base.F90:95
Generated by doxygen 1.14.0