Skip to content

Commit 7bdae65

Browse files
hiromatsuihiromatsui
committed
Add test data for equatorial symmetry fields
1 parent 2058bde commit 7bdae65

35 files changed

+1190
-0
lines changed

‎src/Fortran_libraries/MHD_src/sph_MHD/cal_buoyancy_flux_rtp.f90‎

Lines changed: 181 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,181 @@
1+
!> @file cal_buoyancy_flux_rtp.f90
2+
!! module cal_buoyancy_flux_rtp
3+
!!
4+
!! @author H. Matsui
5+
!! @date Programmed in Oct., 2009
6+
!! @n Modified in Apr., 2013
7+
!
8+
!> @brief Evaluate energy fluxes for MHD dynamo in physical space
9+
!!
10+
!!@verbatim
11+
!! subroutine s_cal_buoyancy_flux_rtp &
12+
!! & (sph_rtp, fl_prop, ref_param_T, ref_param_C, &
13+
!! & bs_trns_base, bs_trns_scalar, fs_trns_eflux, &
14+
!! & trns_b_snap, trns_b_scl, trns_f_eflux)
15+
!! type(sph_rtp_grid), intent(in) :: sph_rtp
16+
!! type(fluid_property), intent(in) :: fl_prop
17+
!! type(reference_scalar_param), intent(in) :: ref_param_T
18+
!! type(reference_scalar_param), intent(in) :: ref_param_C
19+
!! type(base_field_address), intent(in) :: bs_trns_base
20+
!! type(base_field_address), intent(in) :: bs_trns_scalar
21+
!! type(energy_flux_address), intent(in) :: fs_trns_eflux
22+
!! type(spherical_transform_data), intent(in) :: trns_b_snap
23+
!! type(spherical_transform_data), intent(in) :: trns_b_scl
24+
!! type(spherical_transform_data), intent(inout) :: trns_f_eflux
25+
!!@endverbatim
26+
!
27+
module cal_buoyancy_flux_rtp
28+
!
29+
use m_precision
30+
use m_constants
31+
use m_machine_parameter
32+
!
33+
use t_phys_address
34+
use t_spheric_rtp_data
35+
use t_physical_property
36+
use t_reference_scalar_param
37+
use t_addresses_sph_transform
38+
use t_schmidt_poly_on_rtm
39+
!
40+
implicit none
41+
!
42+
private :: sel_buoyancy_flux_rtp, sel_pole_sph_buoyancy_flux
43+
!
44+
! -----------------------------------------------------------------------
45+
!
46+
contains
47+
!
48+
! -----------------------------------------------------------------------
49+
!
50+
subroutine s_cal_buoyancy_flux_rtp &
51+
& (sph_rtp, fl_prop, ref_param_T, ref_param_C, &
52+
& bs_trns_base, bs_trns_scalar, fs_trns_eflux, &
53+
& trns_b_snap, trns_b_scl, trns_f_eflux)
54+
!
55+
use cal_self_buoyancies_sph
56+
!
57+
type(sph_rtp_grid), intent(in) :: sph_rtp
58+
type(fluid_property), intent(in) :: fl_prop
59+
type(reference_scalar_param), intent(in) :: ref_param_T
60+
type(reference_scalar_param), intent(in) :: ref_param_C
61+
type(base_field_address), intent(in) :: bs_trns_base
62+
type(base_field_address), intent(in) :: bs_trns_scalar
63+
type(energy_flux_address), intent(in) :: fs_trns_eflux
64+
type(spherical_transform_data), intent(in) :: trns_b_snap
65+
type(spherical_transform_data), intent(in) :: trns_b_scl
66+
!
67+
type(spherical_transform_data), intent(inout) :: trns_f_eflux
68+
!
69+
integer(kind = kint) :: ibuo_temp, ibuo_comp
70+
!
71+
!
72+
call sel_field_address_for_buoyancies &
73+
& (bs_trns_scalar, ref_param_T, ref_param_C, &
74+
& ibuo_temp, ibuo_comp)
75+
!
76+
if(fs_trns_eflux%i_t_buo_gen .gt. 0) then
77+
call sel_buoyancy_flux_rtp(fl_prop%i_grav, sph_rtp, &
78+
& fl_prop%coef_buo, trns_b_scl%fld_rtp(1,ibuo_temp), &
79+
& trns_b_snap%fld_rtp(1,bs_trns_base%i_velo), &
80+
& trns_f_eflux%fld_rtp(1,fs_trns_eflux%i_t_buo_gen))
81+
call sel_pole_sph_buoyancy_flux &
82+
& (fl_prop%i_grav, sph_rtp%nnod_pole, sph_rtp%nidx_rtp(1), &
83+
& sph_rtp%radius_1d_rtp_r, fl_prop%coef_buo, &
84+
& trns_b_scl%fld_pole(1,ibuo_temp), &
85+
& trns_b_snap%fld_pole(1,bs_trns_base%i_velo), &
86+
& trns_f_eflux%fld_pole(1,fs_trns_eflux%i_t_buo_gen))
87+
end if
88+
!
89+
if(fs_trns_eflux%i_c_buo_gen .gt. 0) then
90+
call sel_buoyancy_flux_rtp(fl_prop%i_grav, sph_rtp, &
91+
& fl_prop%coef_comp_buo, trns_b_scl%fld_rtp(1,ibuo_comp), &
92+
& trns_b_snap%fld_rtp(1,bs_trns_base%i_velo), &
93+
& trns_f_eflux%fld_rtp(1,fs_trns_eflux%i_c_buo_gen))
94+
call sel_pole_sph_buoyancy_flux &
95+
& (fl_prop%i_grav, sph_rtp%nnod_pole, sph_rtp%nidx_rtp(1), &
96+
& sph_rtp%radius_1d_rtp_r, fl_prop%coef_comp_buo, &
97+
& trns_b_scl%fld_pole(1,ibuo_comp), &
98+
& trns_b_snap%fld_pole(1,bs_trns_base%i_velo), &
99+
& trns_f_eflux%fld_pole(1,fs_trns_eflux%i_c_buo_gen))
100+
end if
101+
!
102+
call cal_total_buoyancy_scalar &
103+
& (fs_trns_eflux%i_t_buo_gen, fs_trns_eflux%i_c_buo_gen, &
104+
& fs_trns_eflux%i_buoyancy_flux, sph_rtp%nnod_rtp, &
105+
& trns_f_eflux%ncomp, trns_f_eflux%fld_rtp)
106+
call cal_total_buoyancy_scalar &
107+
& (fs_trns_eflux%i_t_buo_gen, fs_trns_eflux%i_c_buo_gen, &
108+
& fs_trns_eflux%i_buoyancy_flux, sph_rtp%nnod_pole, &
109+
& trns_f_eflux%ncomp, trns_f_eflux%fld_pole)
110+
!
111+
end subroutine s_cal_buoyancy_flux_rtp
112+
!
113+
!-----------------------------------------------------------------------
114+
! -----------------------------------------------------------------------
115+
!
116+
subroutine sel_buoyancy_flux_rtp(i_grav, sph_rtp, coef, &
117+
& scalar, vr, prod)
118+
!
119+
use cal_sph_buoyancy_flux
120+
use cal_products_smp
121+
!
122+
type(sph_rtp_grid), intent(in) :: sph_rtp
123+
integer(kind = kint), intent(in) :: i_grav
124+
real(kind=kreal), intent(in) :: coef
125+
real(kind=kreal), intent(in) :: scalar(sph_rtp%nnod_rtp)
126+
real(kind=kreal), intent(in) :: vr(sph_rtp%nnod_rtp)
127+
!
128+
real(kind=kreal), intent(inout) :: prod(sph_rtp%nnod_rtp)
129+
!
130+
!
131+
if(i_grav .eq. iflag_radial_g) then
132+
call cal_scalar_product_w_coef(sph_rtp%nnod_rtp, coef, &
133+
& scalar, vr, prod)
134+
else
135+
if(sph_rtp%istep_rtp(1) .eq. 1) then
136+
call sph_self_buoyancy_flux_rin &
137+
& (sph_rtp%nnod_rtp, sph_rtp%nidx_rtp, &
138+
& sph_rtp%radius_1d_rtp_r, coef, scalar, vr, prod)
139+
else
140+
call sph_self_buoyancy_flux_pin &
141+
& (sph_rtp%nnod_rtp, sph_rtp%nidx_rtp, &
142+
& sph_rtp%radius_1d_rtp_r, coef, scalar, vr, prod)
143+
end if
144+
end if
145+
!
146+
end subroutine sel_buoyancy_flux_rtp
147+
!
148+
!-----------------------------------------------------------------------
149+
!-----------------------------------------------------------------------
150+
!
151+
subroutine sel_pole_sph_buoyancy_flux &
152+
& (i_grav, nnod_pole, nidx_rtp_r, radius, coef, &
153+
& t_pole, v_pole, d_pole)
154+
!
155+
use cal_sph_buoyancy_flux
156+
!
157+
integer(kind = kint), intent(in) :: i_grav
158+
integer(kind = kint), intent(in) :: nnod_pole
159+
integer(kind = kint), intent(in) :: nidx_rtp_r
160+
real(kind=kreal), intent(in) :: radius(nidx_rtp_r)
161+
!
162+
real(kind = kreal), intent(in) :: coef
163+
real(kind = kreal), intent(in) :: t_pole(nnod_pole)
164+
real(kind = kreal), intent(in) :: v_pole(nnod_pole,3)
165+
!
166+
real(kind = kreal), intent(inout) :: d_pole(nnod_pole)
167+
!
168+
!
169+
if(i_grav .eq. iflag_radial_g) then
170+
call pole_sph_const_buoyancy_flux(nnod_pole, nidx_rtp_r, coef, &
171+
& t_pole, v_pole, d_pole)
172+
else
173+
call pole_sph_self_buoyancy_flux(nnod_pole, nidx_rtp_r, radius, &
174+
& coef, t_pole, v_pole, d_pole)
175+
end if
176+
!
177+
end subroutine sel_pole_sph_buoyancy_flux
178+
!
179+
! -----------------------------------------------------------------------
180+
!
181+
end module cal_buoyancy_flux_rtp

‎src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_buoyancy_flux.f90‎

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
!> @file cal_sph_buoyancy_flux.f90
2+
!! module cal_sph_buoyancy_flux
3+
!!
4+
!! @author H. Matsui
5+
!! @date Programmed in Oct., 2009
6+
!! @n Modified in Apr., 2013
7+
!
8+
!> @brief Evaluate energy fluxes for MHD dynamo in physical space
9+
!!
10+
!!@verbatim
11+
!! subroutine sph_self_buoyancy_flux_rin(nnod_rtp, nidx_rtp, &
12+
!! & radius, coef, scalar, vr, prod)
13+
!! subroutine sph_self_buoyancy_flux_pin(nnod_rtp, nidx_rtp, &
14+
!! & radius, coef, scalar, vr, prod)
15+
!! integer(kind = kint), intent(in) :: nnod_rtp
16+
!! integer(kind = kint), intent(in) :: nidx_rtp(3)
17+
!! real(kind=kreal), intent(in) :: coef
18+
!! real(kind=kreal), intent(in) :: scalar(nnod_rtp), vr(nnod_rtp)
19+
!! real(kind=kreal), intent(in) :: radius(nidx_rtp(1))
20+
!! real(kind=kreal), intent(inout) :: prod(nnod_rtp)
21+
!!
22+
!! subroutine pole_sph_self_buoyancy_flux(nnod_pole, nidx_rtp_r, &
23+
!! & radius, coef, t_pole, v_pole, d_pole)
24+
!! subroutine pole_sph_const_buoyancy_flux(nnod_pole, nidx_rtp_r, &
25+
!! & coef, t_pole, v_pole, d_pole)
26+
!! integer(kind = kint), intent(in) :: nnod_pole
27+
!! integer(kind = kint), intent(in) :: nidx_rtp_r
28+
!! real(kind=kreal), intent(in) :: radius(nidx_rtp_r)
29+
!! real(kind = kreal), intent(in) :: coef
30+
!! real(kind = kreal), intent(in) :: t_pole(nnod_pole)
31+
!! real(kind = kreal), intent(in) :: v_pole(nnod_pole,3)
32+
!! real(kind = kreal), intent(inout) :: d_pole(nnod_pole)
33+
!!@endverbatim
34+
!
35+
module cal_sph_buoyancy_flux
36+
!
37+
use m_precision
38+
use m_constants
39+
use m_machine_parameter
40+
!
41+
implicit none
42+
!
43+
! -----------------------------------------------------------------------
44+
!
45+
contains
46+
!
47+
! -----------------------------------------------------------------------
48+
!
49+
subroutine sph_self_buoyancy_flux_rin(nnod_rtp, nidx_rtp, &
50+
& radius, coef, scalar, vr, prod)
51+
!
52+
integer(kind = kint), intent(in) :: nnod_rtp
53+
integer(kind = kint), intent(in) :: nidx_rtp(3)
54+
real(kind=kreal), intent(in) :: coef
55+
real(kind=kreal), intent(in) :: scalar(nnod_rtp), vr(nnod_rtp)
56+
real(kind=kreal), intent(in) :: radius(nidx_rtp(1))
57+
!
58+
real(kind=kreal), intent(inout) :: prod(nnod_rtp)
59+
!
60+
integer (kind=kint) :: inod, k, ml
61+
!
62+
!
63+
!$omp parallel do private(ml,k,inod)
64+
do ml = 1, nidx_rtp(2)*nidx_rtp(3)
65+
do k = 1, nidx_rtp(1)
66+
inod = k + (ml-1) * nidx_rtp(1)
67+
prod(inod) = coef*scalar(inod)*vr(inod)*radius(k)
68+
end do
69+
end do
70+
!$omp end parallel do
71+
!
72+
end subroutine sph_self_buoyancy_flux_rin
73+
!
74+
!-----------------------------------------------------------------------
75+
!
76+
subroutine sph_self_buoyancy_flux_pin(nnod_rtp, nidx_rtp, &
77+
& radius, coef, scalar, vr, prod)
78+
!
79+
integer(kind = kint), intent(in) :: nnod_rtp
80+
integer(kind = kint), intent(in) :: nidx_rtp(3)
81+
real(kind=kreal), intent(in) :: coef
82+
real(kind=kreal), intent(in) :: scalar(nnod_rtp), vr(nnod_rtp)
83+
real(kind=kreal), intent(in) :: radius(nidx_rtp(1))
84+
!
85+
real(kind=kreal), intent(inout) :: prod(nnod_rtp)
86+
!
87+
integer (kind=kint) :: inod, k, l, m
88+
!
89+
!
90+
!$omp parallel private(l,k)
91+
do l = 1, nidx_rtp(2)
92+
do k = 1, nidx_rtp(1)
93+
!$omp do private(m,inod)
94+
do m = 1, nidx_rtp(3)
95+
inod = m + (k-1) * nidx_rtp(3) &
96+
& + (l-1) * nidx_rtp(3) * nidx_rtp(1)
97+
prod(inod) = coef*scalar(inod)*vr(inod)*radius(k)
98+
end do
99+
!$omp end do nowait
100+
end do
101+
end do
102+
!$omp end parallel
103+
!
104+
end subroutine sph_self_buoyancy_flux_pin
105+
!
106+
!-----------------------------------------------------------------------
107+
! -----------------------------------------------------------------------
108+
!
109+
subroutine pole_sph_self_buoyancy_flux(nnod_pole, nidx_rtp_r, &
110+
& radius, coef, t_pole, v_pole, d_pole)
111+
!
112+
integer(kind = kint), intent(in) :: nnod_pole
113+
integer(kind = kint), intent(in) :: nidx_rtp_r
114+
real(kind=kreal), intent(in) :: radius(nidx_rtp_r)
115+
!
116+
real(kind = kreal), intent(in) :: coef
117+
real(kind = kreal), intent(in) :: t_pole(nnod_pole)
118+
real(kind = kreal), intent(in) :: v_pole(nnod_pole,3)
119+
!
120+
real(kind = kreal), intent(inout) :: d_pole(nnod_pole)
121+
!
122+
integer(kind = kint) :: inod, kr
123+
!
124+
!
125+
! field for north pole (kr) and south pole (inod)
126+
!$omp parallel do private(kr,inod)
127+
do kr = 1, nidx_rtp_r
128+
inod = kr + nidx_rtp_r
129+
d_pole(kr) = coef*t_pole(kr) * v_pole(kr,3) * radius(kr)
130+
d_pole(inod) = -coef*t_pole(inod)*v_pole(inod,3)*radius(kr)
131+
end do
132+
!$omp end parallel do
133+
!
134+
d_pole(2*nidx_rtp_r+1) = zero
135+
!
136+
end subroutine pole_sph_self_buoyancy_flux
137+
!
138+
! -----------------------------------------------------------------------
139+
!
140+
subroutine pole_sph_const_buoyancy_flux(nnod_pole, nidx_rtp_r, &
141+
& coef, t_pole, v_pole, d_pole)
142+
!
143+
integer(kind = kint), intent(in) :: nnod_pole
144+
integer(kind = kint), intent(in) :: nidx_rtp_r
145+
!
146+
real(kind = kreal), intent(in) :: coef
147+
real(kind = kreal), intent(in) :: t_pole(nnod_pole)
148+
real(kind = kreal), intent(in) :: v_pole(nnod_pole,3)
149+
!
150+
real(kind = kreal), intent(inout) :: d_pole(nnod_pole)
151+
!
152+
integer(kind = kint) :: inod, kr
153+
!
154+
!
155+
! field for north pole (kr) and south pole (inod)
156+
!$omp parallel do private(kr,inod)
157+
do kr = 1, nidx_rtp_r
158+
inod = kr + nidx_rtp_r
159+
d_pole(kr) &
160+
& = coef*t_pole(kr)*v_pole(kr,3)
161+
d_pole(inod) &
162+
& = -coef*t_pole(inod)*v_pole(inod,3)
163+
end do
164+
!$omp end parallel do
165+
!
166+
d_pole(2*nidx_rtp_r+1) = zero
167+
!
168+
end subroutine pole_sph_const_buoyancy_flux
169+
!
170+
! -----------------------------------------------------------------------
171+
!
172+
end module cal_sph_buoyancy_flux

0 commit comments

Comments
 (0)