Skip to content

Commit 6a600f9

Browse files
hiromatsuihiromatsui
committed
Make new module for scalar diffusion with boundary condition at center
1 parent 8c76415 commit 6a600f9

16 files changed

+1075
-454
lines changed

‎src/Fortran_libraries/MHD_src/radial_FDM/sph_exp_fix_flx_diffuse_ICB.f90‎

Lines changed: 222 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,222 @@
1+
!>@file sph_exp_fix_flx_diffuse_ICB.f90
2+
!!@brief module sph_exp_fix_flx_diffuse_ICB
3+
!!
4+
!!@author H. Matsui
5+
!!@date Programmed in Jan, 2010
6+
!
7+
!>@brief Evaluate scalar fields using fixed flux condition
8+
!!
9+
!!@verbatim
10+
!! subroutine sph_in_fix_flux_scl_diffuse2(nnod_rj, jmax, g_sph_rj,&
11+
!! & kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
12+
!! & coef_d, scl_rj, dfs_rj)
13+
!! subroutine sph_in_fix_flux_val_diffuse2(nnod_rj, jmax, g_sph_rj,&
14+
!! & kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
15+
!! & coef_d, k_ratio, dk_dr, scl_rj, dfs_rj)
16+
!! integer (kind = kint), intent(in) :: nnod_rj, jmax
17+
!! integer(kind = kint), intent(in) :: kr_in
18+
!! real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
19+
!! real(kind = kreal), intent(in) :: coef_d
20+
!! real(kind = kreal), intent(in) :: k_ratio, dk_dr
21+
!! real(kind = kreal), intent(in) :: flux_ICB(jmax)
22+
!! real(kind = kreal), intent(in) :: r_ICB(0:2)
23+
!! real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
24+
!! real(kind = kreal), intent(in) :: scl_rj(nnod_rj)
25+
!! real(kind = kreal), intent(inout) :: dfs_rj(nnod_rj)
26+
!!
27+
!! subroutine adjust_sph_in_fix_flux(nnod_rj, jmax, &
28+
!! & kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
29+
!! & coef_d, coef_imp, dt, scl_rj)
30+
!! subroutine adjust_sph_in_fix_flx_v_diff(nnod_rj, jmax, &
31+
!! & kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
32+
!! & coef_d, k_ratio, dk_dr, coef_imp, dt, scl_rj)
33+
!! integer(kind = kint), intent(in) :: nnod_rj, jmax
34+
!! integer(kind = kint), intent(in) :: kr_in
35+
!! real(kind = kreal), intent(in) :: coef_imp, dt
36+
!! real(kind = kreal), intent(in) :: coef_d
37+
!! real(kind = kreal), intent(in) :: k_ratio, dk_dr
38+
!! real(kind = kreal), intent(in) :: flux_ICB(jmax)
39+
!! real(kind = kreal), intent(in) :: r_ICB(0:2)
40+
!! real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
41+
!! real (kind=kreal), intent(inout) :: scl_rj(nnod_rj)
42+
!!@endverbatim
43+
!!
44+
!!@n @param idx_rj_degree_zero Local address for degree 0
45+
!!@n @param n_point Number of points for spectrum data
46+
!!@n @param jmax Number of modes for spherical harmonics @f$L*(L+2)@f$
47+
!!@n @param g_sph_rj(jmax,13) Normalization coefficients
48+
!!@n @param j0 Local harmonics mode address for l = m = 0
49+
!!@n @param kr_in Radial ID for inner boundary
50+
!!@n @param r_ICB(0:2) Radius at ICB
51+
!!@n @param flux_ICB(jamx) Spectrum of fixed flux at ICB
52+
!!@n @param fdm2_fix_dr_ICB(-1:1,3)
53+
!! Matrix to evaluate field at ICB with fixed radial derivative
54+
!!
55+
!!@n @param coef_d Coefficient for diffusion term
56+
!!
57+
!!@n @param is_fld Address of spectrum data d_rj
58+
!! (poloidal component for vector)
59+
!!@n @param is_grd Address of gradient of spectrum data d_rj
60+
!! (poloidal component)
61+
!!@n @param is_div Address of divergence of spectrum data d_rj
62+
!!@n @param is_diffuse Address of divergence of spectrum data d_rj
63+
!! (poloidal component for vector)
64+
!!
65+
!!@n @param ntot_phys_rj Total number of components
66+
!!@n @param d_rj Spectrum data
67+
!
68+
module sph_exp_fix_flx_diffuse_ICB
69+
!
70+
use m_precision
71+
!
72+
use m_constants
73+
!
74+
implicit none
75+
!
76+
! -----------------------------------------------------------------------
77+
!
78+
contains
79+
!
80+
! -----------------------------------------------------------------------
81+
!
82+
subroutine sph_in_fix_flux_scl_diffuse2(nnod_rj, jmax, g_sph_rj, &
83+
& kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
84+
& coef_d, scl_rj, dfs_rj)
85+
!
86+
integer (kind = kint), intent(in) :: nnod_rj, jmax
87+
integer(kind = kint), intent(in) :: kr_in
88+
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
89+
real(kind = kreal), intent(in) :: coef_d
90+
real(kind = kreal), intent(in) :: flux_ICB(jmax)
91+
real(kind = kreal), intent(in) :: r_ICB(0:2)
92+
real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
93+
real(kind = kreal), intent(in) :: scl_rj(nnod_rj)
94+
!
95+
real(kind = kreal), intent(inout) :: dfs_rj(nnod_rj)
96+
!
97+
real(kind = kreal) :: d2t_dr2
98+
integer(kind = kint) :: inod, i_p1, j
99+
!
100+
!
101+
!$omp parallel do private(inod,i_p1,d2t_dr2)
102+
do j = 1, jmax
103+
inod = j + (kr_in-1) * jmax
104+
i_p1 = inod + jmax
105+
!
106+
d2t_dr2 = fdm2_fix_dr_ICB(-1,3) * flux_ICB(j) &
107+
& + fdm2_fix_dr_ICB( 0,3) * scl_rj(inod) &
108+
& + fdm2_fix_dr_ICB( 1,3) * scl_rj(i_p1)
109+
!
110+
dfs_rj(inod) = coef_d * (d2t_dr2 + two*r_ICB(1) * flux_ICB(j) &
111+
& - g_sph_rj(j,3)*r_ICB(2) * scl_rj(inod))
112+
!
113+
end do
114+
!$omp end parallel do
115+
!
116+
end subroutine sph_in_fix_flux_scl_diffuse2
117+
!
118+
! -----------------------------------------------------------------------
119+
!
120+
subroutine sph_in_fix_flux_val_diffuse2(nnod_rj, jmax, g_sph_rj, &
121+
& kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
122+
& coef_d, k_ratio, dk_dr, scl_rj, dfs_rj)
123+
!
124+
integer (kind = kint), intent(in) :: nnod_rj, jmax
125+
integer(kind = kint), intent(in) :: kr_in
126+
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
127+
real(kind = kreal), intent(in) :: coef_d
128+
real(kind = kreal), intent(in) :: k_ratio, dk_dr
129+
real(kind = kreal), intent(in) :: flux_ICB(jmax)
130+
real(kind = kreal), intent(in) :: r_ICB(0:2)
131+
real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
132+
real(kind = kreal), intent(in) :: scl_rj(nnod_rj)
133+
!
134+
real(kind = kreal), intent(inout) :: dfs_rj(nnod_rj)
135+
!
136+
real(kind = kreal) :: d2t_dr2
137+
integer(kind = kint) :: inod, i_p1, j
138+
!
139+
!
140+
!$omp parallel do private(inod,i_p1,d2t_dr2)
141+
do j = 1, jmax
142+
inod = j + (kr_in-1) * jmax
143+
i_p1 = inod + jmax
144+
!
145+
d2t_dr2 = fdm2_fix_dr_ICB(-1,3) * flux_ICB(j) &
146+
& + fdm2_fix_dr_ICB( 0,3) * scl_rj(inod) &
147+
& + fdm2_fix_dr_ICB( 1,3) * scl_rj(i_p1)
148+
!
149+
dfs_rj(inod) = coef_d * (d2t_dr2 + two*r_ICB(1) * flux_ICB(j) &
150+
& - g_sph_rj(j,3)*r_ICB(2) * scl_rj(inod))
151+
!
152+
end do
153+
!$omp end parallel do
154+
!
155+
end subroutine sph_in_fix_flux_val_diffuse2
156+
!
157+
! -----------------------------------------------------------------------
158+
! -----------------------------------------------------------------------
159+
!
160+
subroutine adjust_sph_in_fix_flux(nnod_rj, jmax, &
161+
& kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
162+
& coef_d, coef_imp, dt, scl_rj)
163+
!
164+
integer(kind = kint), intent(in) :: nnod_rj, jmax
165+
integer(kind = kint), intent(in) :: kr_in
166+
real(kind = kreal), intent(in) :: coef_d, coef_imp, dt
167+
real(kind = kreal), intent(in) :: flux_ICB(jmax)
168+
real(kind = kreal), intent(in) :: r_ICB(0:2)
169+
real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
170+
!
171+
real (kind=kreal), intent(inout) :: scl_rj(nnod_rj)
172+
!
173+
integer(kind = kint) :: inod, j
174+
!
175+
!
176+
!$omp parallel do private(inod)
177+
do j = 1, jmax
178+
inod = j + (kr_in-1) * jmax
179+
!
180+
scl_rj(inod) = scl_rj(inod) + dt * coef_imp * coef_d &
181+
& * (fdm2_fix_dr_ICB(-1,3) + two*r_ICB(1)) &
182+
& * flux_ICB(j)
183+
end do
184+
!$omp end parallel do
185+
!
186+
end subroutine adjust_sph_in_fix_flux
187+
!
188+
! -----------------------------------------------------------------------
189+
!
190+
subroutine adjust_sph_in_fix_flx_v_diff(nnod_rj, jmax, &
191+
& kr_in, r_ICB, fdm2_fix_dr_ICB, flux_ICB, &
192+
& coef_d, k_ratio, dk_dr, coef_imp, dt, scl_rj)
193+
!
194+
integer(kind = kint), intent(in) :: nnod_rj, jmax
195+
integer(kind = kint), intent(in) :: kr_in
196+
real(kind = kreal), intent(in) :: coef_imp, dt
197+
real(kind = kreal), intent(in) :: coef_d
198+
real(kind = kreal), intent(in) :: k_ratio, dk_dr
199+
real(kind = kreal), intent(in) :: flux_ICB(jmax)
200+
real(kind = kreal), intent(in) :: r_ICB(0:2)
201+
real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
202+
!
203+
real (kind=kreal), intent(inout) :: scl_rj(nnod_rj)
204+
!
205+
integer(kind = kint) :: inod, j
206+
!
207+
!
208+
!$omp parallel do private(inod)
209+
do j = 1, jmax
210+
inod = j + (kr_in-1) * jmax
211+
!
212+
scl_rj(inod) = scl_rj(inod) + dt * coef_imp * coef_d &
213+
& * (k_ratio * (fdm2_fix_dr_ICB(-1,3) &
214+
& + two*r_ICB(1)) + dk_dr) * flux_ICB(j)
215+
end do
216+
!$omp end parallel do
217+
!
218+
end subroutine adjust_sph_in_fix_flx_v_diff
219+
!
220+
! -----------------------------------------------------------------------
221+
!
222+
end module sph_exp_fix_flx_diffuse_ICB

‎src/Fortran_libraries/MHD_src/radial_FDM/sph_exp_fix_scalar_ICB.f90‎

Lines changed: 0 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,6 @@
1313
!! subroutine dsdr_sph_lm0_fix_scalar_in_2(idx_rj_degree_zero, &
1414
!! & jmax, kr_in, r_ICB, fdm2_fix_fld_ICB, fix_ICB, &
1515
!! & is_fld, is_grd, n_point, ntot_phys_rj, d_rj)
16-
!! subroutine cal_sph_fix_scalar_in_diffuse2(jmax, g_sph_rj, &
17-
!! & kr_in, r_ICB, fdm2_fix_fld_ICB, fix_ICB, coef_d, &
18-
!! & is_fld, is_diffuse, n_point, ntot_phys_rj, d_rj)
1916
!! subroutine cal_dsdr_sph_no_bc_in_2 &
2017
!! & (jmax, kr_in, fdm2_fix_fld_ICB, is_fld, is_grd, &
2118
!! & n_point, ntot_phys_rj, d_rj)
@@ -34,11 +31,9 @@
3431
!! Matrix to evaluate radial derivative at ICB with fixed field
3532
!!
3633
!!@n @param fix_ICB(jmax) Spectr data for fixed scalar at ICB
37-
!!@n @param coef_d Coefficient for diffusion term
3834
!!
3935
!!@n @param is_fld Field address of input field
4036
!!@n @param is_grd Field address of radial gradient of field
41-
!!@n @param is_diffuse Field address for diffusion of field
4237
!!
4338
!!@n @param ntot_phys_rj Total number of components
4439
!!@n @param d_rj Spectrum data
@@ -129,50 +124,6 @@ subroutine dsdr_sph_lm0_fix_scalar_in_2(idx_rj_degree_zero, &
129124
end subroutine dsdr_sph_lm0_fix_scalar_in_2
130125
!
131126
! -----------------------------------------------------------------------
132-
!
133-
subroutine cal_sph_fix_scalar_in_diffuse2(jmax, g_sph_rj, &
134-
& kr_in, r_ICB, fdm2_fix_fld_ICB, fix_ICB, coef_d, &
135-
& is_fld, is_diffuse, n_point, ntot_phys_rj, d_rj)
136-
!
137-
integer(kind = kint), intent(in) :: is_fld, is_diffuse
138-
integer(kind = kint), intent(in) :: jmax, kr_in
139-
integer(kind = kint), intent(in) :: n_point, ntot_phys_rj
140-
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
141-
real(kind = kreal), intent(in) :: fix_ICB(jmax)
142-
real(kind = kreal), intent(in) :: r_ICB(0:2)
143-
real(kind = kreal), intent(in) :: fdm2_fix_fld_ICB(0:2,3)
144-
real(kind = kreal), intent(in) :: coef_d
145-
!
146-
real(kind = kreal), intent(inout) :: d_rj(n_point,ntot_phys_rj)
147-
!
148-
real(kind = kreal) :: d1t_dr1, d2t_dr2
149-
integer(kind = kint) :: inod, i_p1, i_p2, j
150-
!
151-
!
152-
!$omp parallel do private(inod,i_p1,i_p2,d1t_dr1,d2t_dr2)
153-
do j = 1, jmax
154-
inod = j + (kr_in-1) * jmax
155-
i_p1 = inod + jmax
156-
i_p2 = i_p1 + jmax
157-
!
158-
d1t_dr1 = fdm2_fix_fld_ICB( 0,2) * fix_ICB(j) &
159-
& + fdm2_fix_fld_ICB( 1,2) * d_rj(i_p1,is_fld) &
160-
& + fdm2_fix_fld_ICB( 2,2) * d_rj(i_p2,is_fld)
161-
d2t_dr2 = fdm2_fix_fld_ICB( 0,3) * fix_ICB(j) &
162-
& + fdm2_fix_fld_ICB( 1,3) * d_rj(i_p1,is_fld) &
163-
& + fdm2_fix_fld_ICB( 2,3) * d_rj(i_p2,is_fld)
164-
!
165-
d_rj(inod,is_fld) = fix_ICB(j)
166-
d_rj(inod,is_diffuse) &
167-
& = coef_d * (d2t_dr2 + two*r_ICB(1)*d1t_dr1 &
168-
& - g_sph_rj(j,3)*r_ICB(2) * d_rj(inod,is_fld) )
169-
!
170-
end do
171-
!$omp end parallel do
172-
!
173-
end subroutine cal_sph_fix_scalar_in_diffuse2
174-
!
175-
! -----------------------------------------------------------------------
176127
!
177128
subroutine cal_dsdr_sph_no_bc_in_2 &
178129
& (jmax, kr_in, fdm2_fix_fld_ICB, is_fld, is_grd, &

‎src/Fortran_libraries/MHD_src/radial_FDM/sph_exp_fix_scl_diffuse_CMB.f90‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -67,25 +67,25 @@ subroutine sph_out_fix_scalar_diffuse2(nnod_rj, jmax, g_sph_rj, &
6767
real(kind = kreal), intent(inout) :: scl_rj(nnod_rj)
6868
real(kind = kreal), intent(inout) :: dfs_rj(nnod_rj)
6969
!
70-
real(kind = kreal) :: d1t_dr1, d2t_dr2
70+
real(kind = kreal) :: d1s_dr1, d2s_dr2
7171
integer(kind = kint) :: inod, i_n1, i_n2, j
7272
!
7373
!
74-
!$omp parallel do private(inod,i_n1,i_n2,d2t_dr2,d1t_dr1)
74+
!$omp parallel do private(inod,i_n1,i_n2,d2s_dr2,d1s_dr1)
7575
do j = 1, jmax
7676
inod = j + (kr_out-1) * jmax
7777
i_n1 = inod - jmax
7878
i_n2 = i_n1 - jmax
7979
!
80-
d1t_dr1 = fdm2_fix_fld_CMB(-2,2) * scl_rj(i_n2) &
80+
d1s_dr1 = fdm2_fix_fld_CMB(-2,2) * scl_rj(i_n2) &
8181
& + fdm2_fix_fld_CMB(-1,2) * scl_rj(i_n1) &
8282
& + fdm2_fix_fld_CMB( 0,2) * fix_CMB(j)
83-
d2t_dr2 = fdm2_fix_fld_CMB(-2,3) * scl_rj(i_n2) &
83+
d2s_dr2 = fdm2_fix_fld_CMB(-2,3) * scl_rj(i_n2) &
8484
& + fdm2_fix_fld_CMB(-1,3) * scl_rj(i_n1) &
8585
& + fdm2_fix_fld_CMB( 0,3) * fix_CMB(j)
8686
!
8787
scl_rj(inod) = fix_CMB(j)
88-
dfs_rj(inod) = coef_d * (d2t_dr2 + two*r_CMB(1) * d1t_dr1 &
88+
dfs_rj(inod) = coef_d * (d2s_dr2 + two*r_CMB(1) * d1s_dr1 &
8989
& - g_sph_rj(j,3)*r_CMB(2) * scl_rj(inod))
9090
!
9191
end do

0 commit comments

Comments
 (0)