Skip to content

Commit 1b4e61f

Browse files
hiromatsuihiromatsui
committed
Add radial diffusivity variation for matrix construction at boundaries
1 parent f87f3e8 commit 1b4e61f

File tree

5 files changed

+140
-72
lines changed

5 files changed

+140
-72
lines changed

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

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,14 @@
1111
!! & r_CTR1, fdm2_fix_fld_ctr1, coef_p, mat3)
1212
!! subroutine add_scalar_poisson_mat_ctr1(nri, jmax, g_sph_rj, &
1313
!! & r_CTR1, fdm2_fix_fld_ctr1, coef_p, mat3)
14+
!! subroutine add_scl_val_diffuse_mat_ctr1 &
15+
!! & (nri, jmax, g_sph_rj, r_CTR1, fdm2_fix_fld_ctr1, &
16+
!! & coef_p, k_ratio, dk_dr, mat3)
1417
!! integer(kind = kint), intent(in) :: jmax, nri
1518
!! real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
1619
!! real(kind = kreal), intent(in) :: r_CTR1(0:2)
1720
!! real(kind = kreal), intent(in) :: coef_p
21+
!! real(kind = kreal), intent(in) :: k_ratio, dk_dr
1822
!! real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
1923
!! real(kind = kreal), intent(inout) :: mat3(3,nri,jmax)
2024
!!
@@ -98,6 +102,38 @@ subroutine add_scalar_poisson_mat_ctr1(nri, jmax, g_sph_rj, &
98102
end subroutine add_scalar_poisson_mat_ctr1
99103
!
100104
! -----------------------------------------------------------------------
105+
!
106+
subroutine add_scl_val_diffuse_mat_ctr1 &
107+
& (nri, jmax, g_sph_rj, r_CTR1, fdm2_fix_fld_ctr1, &
108+
& coef_p, k_ratio, dk_dr, mat3)
109+
!
110+
integer(kind = kint), intent(in) :: jmax, nri
111+
real(kind = kreal), intent(in) :: g_sph_rj(jmax,13)
112+
real(kind = kreal), intent(in) :: r_CTR1(0:2)
113+
real(kind = kreal), intent(in) :: coef_p
114+
real(kind = kreal), intent(in) :: k_ratio, dk_dr
115+
real(kind = kreal), intent(in) :: fdm2_fix_fld_ctr1(-1:1,3)
116+
!
117+
real(kind = kreal), intent(inout) :: mat3(3,nri,jmax)
118+
!
119+
integer(kind = kint) :: j
120+
!
121+
!
122+
do j = 1, jmax
123+
mat3(2,1,j) = mat3(2,1,j) &
124+
& - coef_p * k_ratio * (fdm2_fix_fld_ctr1(0,3) &
125+
& + two*r_CTR1(1) * fdm2_fix_fld_ctr1(0,2) &
126+
& - g_sph_rj(j,3)*r_CTR1(2)) &
127+
& - coef_p * dk_dr * fdm2_fix_fld_ctr1(0,2)
128+
mat3(1,2,j) = mat3(1,2,j) &
129+
& - coef_p * k_ratio * (fdm2_fix_fld_ctr1(1,3) &
130+
& + two*r_CTR1(1) * fdm2_fix_fld_ctr1(1,2)) &
131+
& - coef_p * dk_dr * fdm2_fix_fld_ctr1(1,2)
132+
end do
133+
!
134+
end subroutine add_scl_val_diffuse_mat_ctr1
135+
!
136+
! -----------------------------------------------------------------------
101137
! -----------------------------------------------------------------------
102138
!
103139
subroutine set_unit_mat3_filter_to_center(nri, jmax, ICB_Vspec, &

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

Lines changed: 30 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -9,16 +9,24 @@
99
!! subroutine const_radial_mat_4_press_sph &
1010
!! & (sph_rj, r_2nd, fl_prop, sph_bc_U, fdm2_center, &
1111
!! & g_sph_rj, band_p_poisson)
12-
!! subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect, &
13-
!! & dt, sph_params, sph_rj, r_2nd, property, &
14-
!! & sph_bc, bcs_S, fdm2_center, g_sph_rj, band_s_evo)
15-
!! type(scalar_property), intent(in) :: property
16-
!! type(sph_shell_parameters), intent(in) :: sph_params
17-
!! type(sph_rj_grid), intent(in) :: sph_rj
12+
!! subroutine const_radial_mat_4_scalar_sph(mat_name, dt, sph_rj, &
13+
!! & g_sph_rj, r_2nd, fdm2_center, property, sph_bc, bcs_S,&
14+
!! & k_ratio, dk_dr, band_s_evo)
15+
!! type(sph_rj_grid), intent(in) :: sph_rj
1816
!! type(fdm_matrices), intent(in) :: r_2nd
17+
!! type(scalar_property), intent(in) :: property
18+
!! type(fluid_property), intent(in) :: fl_prop
1919
!! type(sph_boundary_type), intent(in) :: sph_bc_U
20+
!! type(sph_boundary_type), intent(in) :: sph_bc
2021
!! type(sph_scalar_boundary_data) :: bcs_S
2122
!! type(fdm2_center_mat), intent(in) :: fdm2_center
23+
!! real(kind = kreal), intent(in) :: dt
24+
!! real(kind=kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
25+
!! real(kind = kreal), intent(in) :: k_ratio(0:sph_rj%nidx_rj(1))
26+
!! real(kind = kreal), intent(in) :: dk_dr(0:sph_rj%nidx_rj(1))
27+
!! character(len=kchara), intent(in) :: mat_name
28+
!! type(band_matrices_type), intent(inout) :: band_p_poisson
29+
!! type(band_matrices_type), intent(inout) :: band_s_evo
2230
!!
2331
!! subroutine const_r_mat00_scalar_sph(id_file, mat_name, &
2432
!! & diffuse_reduction_ratio_ICB, sph_params, &
@@ -86,30 +94,23 @@ subroutine const_radial_mat_4_press_sph &
8694
type(band_matrices_type), intent(inout) :: band_p_poisson
8795
!
8896
real(kind = kreal) :: coef_p
89-
real(kind = kreal), allocatable :: r_coef(:)
9097
!
9198
!
9299
write(band_p_poisson%mat_name,'(a)') 'pressure_poisson'
93100
coef_p = - fl_prop%coef_press
94-
!
95-
allocate(r_coef(sph_rj%nidx_rj(1)))
96-
!$omp parallel workshare
97-
r_coef(1:sph_rj%nidx_rj(1)) = coef_p
98-
!$omp end parallel workshare
99101
!
100102
call alloc_band_mat_sph(ithree, sph_rj, band_p_poisson)
101-
103+
!
102104
call set_unit_mat_4_poisson &
103105
& (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), &
104106
& sph_bc_U%kr_in, sph_bc_U%kr_out, band_p_poisson%mat)
105107
call add_scalar_poisson_mat_sph &
106108
& (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, &
107-
& g_sph_rj, sph_bc_U%kr_in, sph_bc_U%kr_out, r_coef(1), &
109+
& g_sph_rj, sph_bc_U%kr_in, sph_bc_U%kr_out, coef_p, &
108110
& r_2nd%fdm(1)%dmat, r_2nd%fdm(2)%dmat, band_p_poisson%mat)
109111
!
110112
call sel_radial_mat_press_bc_sph(sph_rj, sph_bc_U, fdm2_center, &
111-
& g_sph_rj, r_coef, band_p_poisson)
112-
deallocate(r_coef)
113+
& g_sph_rj, coef_p, band_p_poisson)
113114
!
114115
call ludcmp_3band_mul_t &
115116
& (np_smp, sph_rj%istack_rj_j_smp, band_p_poisson)
@@ -119,42 +120,38 @@ end subroutine const_radial_mat_4_press_sph
119120
! -----------------------------------------------------------------------
120121
! -----------------------------------------------------------------------
121122
!
122-
subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect, &
123-
& dt, sph_params, sph_rj, r_2nd, property, &
124-
& sph_bc, bcs_S, fdm2_center, g_sph_rj, k_ratio, dk_dr, &
125-
& band_s_evo)
123+
subroutine const_radial_mat_4_scalar_sph(mat_name, dt, sph_rj, &
124+
& g_sph_rj, r_2nd, fdm2_center, property, sph_bc, bcs_S, &
125+
& k_ratio, dk_dr, band_s_evo)
126126
!
127127
use m_ludcmp_3band
128128
use center_sph_matrices
129129
use set_radial_mat_sph
130130
use select_r_mat_scalar_bc_sph
131131
!
132-
type(sph_shell_parameters), intent(in) :: sph_params
133132
type(sph_rj_grid), intent(in) :: sph_rj
134133
type(fdm_matrices), intent(in) :: r_2nd
135134
type(scalar_property), intent(in) :: property
136135
type(sph_boundary_type), intent(in) :: sph_bc
137136
type(sph_scalar_boundary_data) :: bcs_S
138137
type(fdm2_center_mat), intent(in) :: fdm2_center
139138
!
140-
real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
141139
real(kind = kreal), intent(in) :: dt
142-
real(kind = kreal), intent(in) :: coef_advect
140+
real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
143141
real(kind = kreal), intent(in) :: k_ratio(0:sph_rj%nidx_rj(1))
144142
real(kind = kreal), intent(in) :: dk_dr(0:sph_rj%nidx_rj(1))
145143
character(len=kchara), intent(in) :: mat_name
146144
!
147145
type(band_matrices_type), intent(inout) :: band_s_evo
148146
!
149147
real(kind = kreal) :: coef
150-
real(kind = kreal), allocatable :: r_coef(:)
151148
!
152149
!
153150
band_s_evo%mat_name = mat_name
154151
call alloc_band_mat_sph(ithree, sph_rj, band_s_evo)
155152
call set_unit_on_diag(band_s_evo)
156153
!
157-
if(coef_advect .eq. zero) then
154+
if(property%coef_advect .eq. zero) then
158155
coef = one
159156
call set_unit_mat_4_poisson &
160157
& (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), &
@@ -164,38 +161,26 @@ subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect, &
164161
call set_unit_mat_4_time_evo &
165162
& (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_s_evo%mat)
166163
end if
167-
!
168-
allocate(r_coef(sph_rj%nidx_rj(1)))
169-
!$omp parallel workshare
170-
r_coef(1:sph_rj%nidx_rj(1)) = coef
171-
!$omp end parallel workshare
172-
!
173-
if(property%diffuse_reduction_ratio_ICB .lt. one) then
174-
r_coef(sph_params%nlayer_ICB) &
175-
& = property%diffuse_reduction_ratio_ICB &
176-
& * r_coef(sph_params%nlayer_ICB)
177-
if(my_rank .eq. 0) write(*,*) 'reduction of diffusivity at', &
178-
& sph_params%nlayer_ICB, ' to ', r_coef(sph_params%nlayer_ICB), &
179-
& ' from ' , coef
180-
end if
181-
!
182164
!
183165
if(property%flag_val_diffuse) then
184166
call add_scalar_r_diffuse_mat_sph &
185167
& (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, &
186168
& g_sph_rj, sph_bc%kr_in, sph_bc%kr_out, &
187-
& r_coef(1), k_ratio(1), dk_dr(1), &
169+
& coef, k_ratio(1), dk_dr(1), &
188170
& r_2nd%fdm(1)%dmat, r_2nd%fdm(2)%dmat, band_s_evo%mat)
189171
else
190172
call add_scalar_poisson_mat_sph &
191173
& (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, &
192-
& g_sph_rj, sph_bc%kr_in, sph_bc%kr_out, r_coef(1), &
174+
& g_sph_rj, sph_bc%kr_in, sph_bc%kr_out, coef, &
193175
& r_2nd%fdm(1)%dmat, r_2nd%fdm(2)%dmat, band_s_evo%mat)
194176
end if
195177
!
196-
call sel_radial_mat_scalar_bc_sph(sph_rj, sph_bc, bcs_S, &
197-
& fdm2_center, g_sph_rj, r_coef, band_s_evo)
198-
deallocate(r_coef)
178+
call sel_sph_radial_mat_scalar_ICB(property%flag_val_diffuse, &
179+
& sph_rj, sph_bc, bcs_S, fdm2_center, g_sph_rj, &
180+
& coef, k_ratio(sph_bc%kr_in), dk_dr(sph_bc%kr_in), band_s_evo)
181+
call sel_sph_radial_mat_scalar_CMB(property%flag_val_diffuse, &
182+
& sph_rj, sph_bc, bcs_S, fdm2_center, g_sph_rj, &
183+
& coef, k_ratio(sph_bc%kr_out), band_s_evo)
199184
!
200185
call ludcmp_3band_mul_t &
201186
& (np_smp, sph_rj%istack_rj_j_smp, band_s_evo)

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

Lines changed: 4 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -231,11 +231,9 @@ subroutine const_radial_matrices_sph &
231231
end if
232232
!
233233
if(MHD_prop%ht_prop%iflag_scheme .ge. id_Crank_nicolson) then
234-
call const_radial_mat_4_scalar_sph &
235-
& (temp_evo_name, MHD_prop%ht_prop%coef_advect, dt, &
236-
& sph_params, sph_rj, r_2nd, MHD_prop%ht_prop, &
234+
call const_radial_mat_4_scalar_sph(temp_evo_name, dt, sph_rj, &
235+
& g_sph_rj, r_2nd, sph_MHD_bc%fdm2_center, MHD_prop%ht_prop, &
237236
& sph_MHD_bc%sph_bc_T, sph_MHD_bc%bcs_T, &
238-
& sph_MHD_bc%fdm2_center, g_sph_rj, &
239237
& ref_field%d_fld(1,iref_diffusivity%i_T_diffusivity), &
240238
& ref_field%d_fld(1,iref_grad_diffusivity%i_T_diffusivity), &
241239
& sph_MHD_mat%band_temp_evo)
@@ -246,11 +244,9 @@ subroutine const_radial_matrices_sph &
246244
end if
247245
!
248246
if(MHD_prop%cp_prop%iflag_scheme .ge. id_Crank_nicolson) then
249-
call const_radial_mat_4_scalar_sph &
250-
& (comp_evo_name, MHD_prop%cp_prop%coef_advect, dt, &
251-
& sph_params, sph_rj, r_2nd, MHD_prop%cp_prop, &
247+
call const_radial_mat_4_scalar_sph(comp_evo_name, dt, sph_rj, &
248+
& g_sph_rj, r_2nd, sph_MHD_bc%fdm2_center, MHD_prop%cp_prop, &
252249
& sph_MHD_bc%sph_bc_C, sph_MHD_bc%bcs_C, &
253-
& sph_MHD_bc%fdm2_center, g_sph_rj, &
254250
& ref_field%d_fld(1,iref_diffusivity%i_T_diffusivity), &
255251
& ref_field%d_fld(1,iref_grad_diffusivity%i_T_diffusivity), &
256252
& sph_MHD_mat%band_comp_evo)

0 commit comments

Comments
 (0)