Skip to content

Commit cf08b8f

Browse files
hiromatsuihiromatsui
committed
Make new module for scalar advection at center
1 parent 6a600f9 commit cf08b8f

File tree

3 files changed

+143
-109
lines changed

3 files changed

+143
-109
lines changed

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

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
!>@file sph_exp_div_scl_flux_center.f90
2+
!!@brief module sph_exp_div_scl_flux_center
3+
!!
4+
!!@author H. Matsui
5+
!!@date Programmed in Jan., 2010
6+
!
7+
!>@brief Set scalar advection around center
8+
!!
9+
!!@verbatim
10+
!! subroutine sph_div_flux_4_fix_center &
11+
!! & (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
12+
!! & g_sph_rj, r_CTR1, fix_ICB, fdm2_fix_fld_ctr1, &
13+
!! & flx_rj, adv_rj)
14+
!! subroutine sph_div_flux_4_fill_center &
15+
!! & (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
16+
!! & g_sph_rj, r_CTR1, fdm2_fix_fld_ctr1, flx_rj, adv_rj)
17+
!! integer(kind = kint), intent(in) :: inod_rj_center
18+
!! integer(kind = kint), intent(in) :: idx_rj_degree_zero
19+
!! integer(kind = kint), intent(in) :: nnod_rj, jmax
20+
!! real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
21+
!! real(kind = kreal), intent(in) :: r_CTR1(0:2)
22+
!! real(kind = kreal), intent(in) :: fix_ICB(jmax)
23+
!! real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
24+
!! real(kind = kreal), intent(in) :: flx_rj(nnod_rj,3)
25+
!! real(kind = kreal), intent(inout) :: adv_rj(nnod_rj)
26+
!!@endverbatim
27+
!!
28+
!!@n @param inod_rj_center Local address for center
29+
!!@n @param idx_rj_degree_zero Local address for degree 0
30+
!!@n @param jmax Number of local spherical harmonics mode
31+
!!@n @param r_CTR1(0:2) Radius at innermost point
32+
!!@n @param fdm2_fix_fld_ctr1(-1:1,3)
33+
!! Matrix to evaluate radial derivative
34+
!! for center with fixed field
35+
!!@n @param fdm2_fixed_center(0:2,3)
36+
!! Matrix to evaluate radial derivative
37+
!! for center with fixed field
38+
!!@n @param fix_CTR(jmax) Spectr data for fixed scalar at center
39+
!!@n @param coef_d Coefficient for diffusion term
40+
!!
41+
module sph_exp_div_scl_flux_center
42+
!
43+
use m_precision
44+
use m_constants
45+
!
46+
implicit none
47+
!
48+
! -----------------------------------------------------------------------
49+
!
50+
contains
51+
!
52+
! -----------------------------------------------------------------------
53+
!
54+
subroutine sph_div_flux_4_fix_center &
55+
& (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
56+
& g_sph_rj, r_CTR1, fix_ICB, fdm2_fix_fld_ctr1, &
57+
& flx_rj, adv_rj)
58+
!
59+
integer(kind = kint), intent(in) :: inod_rj_center
60+
integer(kind = kint), intent(in) :: idx_rj_degree_zero
61+
integer(kind = kint), intent(in) :: nnod_rj, jmax
62+
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
63+
real(kind = kreal), intent(in) :: r_CTR1(0:2)
64+
real(kind = kreal), intent(in) :: fix_ICB(jmax)
65+
real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
66+
real(kind = kreal), intent(in) :: flx_rj(nnod_rj,3)
67+
!
68+
real(kind = kreal), intent(inout) :: adv_rj(nnod_rj)
69+
!
70+
real(kind = kreal) :: d1s_dr1
71+
integer(kind = kint) :: i_p1, j
72+
!
73+
!
74+
!$omp parallel do private(i_p1,j,d1s_dr1)
75+
do j = 1, jmax
76+
i_p1 = j + jmax
77+
!
78+
d1s_dr1 = fdm2_fix_fld_ctr1(-1,2) * fix_ICB(j) &
79+
& + fdm2_fix_fld_ctr1( 0,2) * flx_rj(j, 1) &
80+
& + fdm2_fix_fld_ctr1( 1,2) * flx_rj(i_p1,1)
81+
!
82+
adv_rj(j) = (d1s_dr1 - flx_rj(j,2) ) &
83+
& * max(g_sph_rj(j,3),half) * r_CTR1(2)
84+
end do
85+
!$omp end parallel do
86+
!
87+
if(inod_rj_center .gt. 0) then
88+
adv_rj(inod_rj_center) = adv_rj(idx_rj_degree_zero)
89+
end if
90+
!
91+
end subroutine sph_div_flux_4_fix_center
92+
!
93+
! -----------------------------------------------------------------------
94+
!
95+
subroutine sph_div_flux_4_fill_center &
96+
& (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
97+
& g_sph_rj, r_CTR1, fdm2_fix_fld_ctr1, flx_rj, adv_rj)
98+
!
99+
integer(kind = kint), intent(in) :: inod_rj_center
100+
integer(kind = kint), intent(in) :: idx_rj_degree_zero
101+
integer(kind = kint), intent(in) :: nnod_rj, jmax
102+
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
103+
real(kind = kreal), intent(in) :: r_CTR1(0:2)
104+
real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
105+
real(kind = kreal), intent(in) :: flx_rj(nnod_rj,3)
106+
!
107+
real(kind = kreal), intent(inout) :: adv_rj(nnod_rj)
108+
!
109+
real(kind = kreal) :: d1s_dr1
110+
integer(kind = kint) :: i_p1, j
111+
!
112+
!
113+
!$omp parallel do private(i_p1,j,d1s_dr1)
114+
do j = 1, jmax
115+
i_p1 = j + jmax
116+
d1s_dr1 = fdm2_fix_fld_ctr1( 0,2) * flx_rj(j, 1) &
117+
& + fdm2_fix_fld_ctr1( 1,2) * flx_rj(i_p1,1)
118+
!
119+
adv_rj(j) = (d1s_dr1 - flx_rj(j,2)) &
120+
& * max(g_sph_rj(j,3),half) * r_CTR1(2)
121+
end do
122+
!$omp end parallel do
123+
!
124+
if(inod_rj_center .le. 0) return
125+
!
126+
i_p1 = idx_rj_degree_zero + jmax
127+
d1s_dr1 = fdm2_fix_fld_ctr1(-1,2) * flx_rj(inod_rj_center,1) &
128+
& + fdm2_fix_fld_ctr1( 0,2) * flx_rj(idx_rj_degree_zero,1) &
129+
& + fdm2_fix_fld_ctr1( 1,2) * flx_rj(i_p1,1)
130+
adv_rj(idx_rj_degree_zero) &
131+
& = half * r_CTR1(2) * (d1s_dr1 - flx_rj(idx_rj_degree_zero,2))
132+
!
133+
!
134+
adv_rj(inod_rj_center) = zero
135+
!
136+
end subroutine sph_div_flux_4_fill_center
137+
!
138+
! -----------------------------------------------------------------------
139+
!
140+
end module sph_exp_div_scl_flux_center

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

Lines changed: 0 additions & 106 deletions
Original file line numberDiff line numberDiff line change
@@ -21,23 +21,6 @@
2121
!! & g_sph_rj, d_center, fdm2_fix_fld_ctr1, &
2222
!! & fdm2_fixed_center, is_fld, is_grd, &
2323
!! & n_point, ntot_phys_rj, d_rj)
24-
!!
25-
!! subroutine cal_sph_div_flux_4_fix_ctr &
26-
!! & (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
27-
!! & g_sph_rj, r_CTR1, fix_ICB, fdm2_fix_fld_ctr1, &
28-
!! & flx_rj, adv_rj)
29-
!! subroutine cal_sph_div_flux_4_fill_ctr &
30-
!! & (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
31-
!! & g_sph_rj, r_CTR1, fdm2_fix_fld_ctr1, flx_rj, adv_rj)
32-
!! integer(kind = kint), intent(in) :: inod_rj_center
33-
!! integer(kind = kint), intent(in) :: idx_rj_degree_zero
34-
!! integer(kind = kint), intent(in) :: nnod_rj, jmax
35-
!! real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
36-
!! real(kind = kreal), intent(in) :: r_CTR1(0:2)
37-
!! real(kind = kreal), intent(in) :: fix_ICB(jmax)
38-
!! real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
39-
!! real(kind = kreal), intent(in) :: flx_rj(nnod_rj,3)
40-
!! real(kind = kreal), intent(inout) :: adv_rj(nnod_rj)
4124
!!@endverbatim
4225
!!
4326
!!@n @param inod_rj_center Local address for center
@@ -209,94 +192,5 @@ subroutine dsdr_sph_lm0_fixed_ctr_2 &
209192
end subroutine dsdr_sph_lm0_fixed_ctr_2
210193
!
211194
! -----------------------------------------------------------------------
212-
! -----------------------------------------------------------------------
213-
!
214-
subroutine cal_sph_div_flux_4_fix_ctr &
215-
& (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
216-
& g_sph_rj, r_CTR1, fix_ICB, fdm2_fix_fld_ctr1, &
217-
& flx_rj, adv_rj)
218-
!
219-
integer(kind = kint), intent(in) :: inod_rj_center
220-
integer(kind = kint), intent(in) :: idx_rj_degree_zero
221-
integer(kind = kint), intent(in) :: nnod_rj, jmax
222-
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
223-
real(kind = kreal), intent(in) :: r_CTR1(0:2)
224-
real(kind = kreal), intent(in) :: fix_ICB(jmax)
225-
real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
226-
real(kind = kreal), intent(in) :: flx_rj(nnod_rj,3)
227-
!
228-
real(kind = kreal), intent(inout) :: adv_rj(nnod_rj)
229-
!
230-
real(kind = kreal) :: d1s_dr1
231-
integer(kind = kint) :: i_p1, j
232-
!
233-
!
234-
!$omp parallel do private(i_p1,j,d1s_dr1)
235-
do j = 1, jmax
236-
i_p1 = j + jmax
237-
!
238-
d1s_dr1 = fdm2_fix_fld_ctr1(-1,2) * fix_ICB(j) &
239-
& + fdm2_fix_fld_ctr1( 0,2) * flx_rj(j, 1) &
240-
& + fdm2_fix_fld_ctr1( 1,2) * flx_rj(i_p1,1)
241-
!
242-
adv_rj(j) = (d1s_dr1 - flx_rj(j,2) ) &
243-
& * max(g_sph_rj(j,3),half) * r_CTR1(2)
244-
end do
245-
!$omp end parallel do
246-
!
247-
if(inod_rj_center .gt. 0) then
248-
adv_rj(inod_rj_center) = adv_rj(idx_rj_degree_zero)
249-
end if
250-
!
251-
end subroutine cal_sph_div_flux_4_fix_ctr
252-
!
253-
! -----------------------------------------------------------------------
254-
!
255-
subroutine cal_sph_div_flux_4_fill_ctr &
256-
& (inod_rj_center, idx_rj_degree_zero, nnod_rj, jmax, &
257-
& g_sph_rj, r_CTR1, fdm2_fix_fld_ctr1, flx_rj, adv_rj)
258-
!
259-
integer(kind = kint), intent(in) :: inod_rj_center
260-
integer(kind = kint), intent(in) :: idx_rj_degree_zero
261-
integer(kind = kint), intent(in) :: nnod_rj, jmax
262-
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
263-
real(kind = kreal), intent(in) :: r_CTR1(0:2)
264-
real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
265-
real(kind = kreal), intent(in) :: flx_rj(nnod_rj,3)
266-
!
267-
real(kind = kreal), intent(inout) :: adv_rj(nnod_rj)
268-
!
269-
real(kind = kreal) :: d1s_dr1
270-
integer(kind = kint) :: inod, i_n1, i_p1, j
271-
!
272-
!
273-
!$omp parallel do private(i_p1,j,d1s_dr1)
274-
do j = 1, jmax
275-
i_p1 = j + jmax
276-
d1s_dr1 = fdm2_fix_fld_ctr1( 0,2) * flx_rj(j, 1) &
277-
& + fdm2_fix_fld_ctr1( 1,2) * flx_rj(i_p1,1)
278-
!
279-
adv_rj(j) = (d1s_dr1 - flx_rj(j,2)) &
280-
& * max(g_sph_rj(j,3),half) * r_CTR1(2)
281-
end do
282-
!$omp end parallel do
283-
!
284-
if(inod_rj_center .le. 0) return
285-
!
286-
i_n1 = inod_rj_center
287-
inod = idx_rj_degree_zero
288-
i_p1 = inod + jmax
289-
!
290-
d1s_dr1 = fdm2_fix_fld_ctr1(-1,2) * flx_rj(i_n1,1) &
291-
& + fdm2_fix_fld_ctr1( 0,2) * flx_rj(inod,1) &
292-
& + fdm2_fix_fld_ctr1( 1,2) * flx_rj(i_p1,1)
293-
adv_rj(inod) = half * r_CTR1(2) * (d1s_dr1 - flx_rj(inod,2))
294-
!
295-
!
296-
adv_rj(i_n1) = zero
297-
!
298-
end subroutine cal_sph_div_flux_4_fill_ctr
299-
!
300-
! -----------------------------------------------------------------------
301195
!
302196
end module cal_sph_exp_center

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

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -349,7 +349,7 @@ subroutine sel_ICB_sph_scalar_advect &
349349
!
350350
use sph_exp_fix_scalar_ICB
351351
use sph_exp_fixed_flux_ICB
352-
use cal_sph_exp_center
352+
use sph_exp_div_scl_flux_center
353353
!
354354
type(sph_rj_grid), intent(in) :: sph_rj
355355
type(fdm2_center_mat), intent(in) :: fdm2_center
@@ -364,13 +364,13 @@ subroutine sel_ICB_sph_scalar_advect &
364364
!
365365
!
366366
if (sph_bc%iflag_icb .eq. iflag_sph_fill_center) then
367-
call cal_sph_div_flux_4_fill_ctr &
367+
call sph_div_flux_4_fill_center &
368368
& (sph_rj%inod_rj_center, sph_rj%idx_rj_degree_zero, &
369369
& sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
370370
& sph_bc%r_ICB, fdm2_center%dmat_fix_fld, &
371371
& d_rj(1,is_flux), d_rj(1,is_advect))
372372
else if(sph_bc%iflag_icb .eq. iflag_sph_fix_center) then
373-
call cal_sph_div_flux_4_fix_ctr &
373+
call sph_div_flux_4_fix_center &
374374
& (sph_rj%inod_rj_center, sph_rj%idx_rj_degree_zero, &
375375
& sph_rj%nnod_rj, sph_rj%nidx_rj(2), g_sph_rj, &
376376
& sph_bc%r_ICB, ICB_Sspec%S_BC, fdm2_center%dmat_fix_fld, &

0 commit comments

Comments
 (0)