Skip to content

Commit 7846b1e

Browse files
hiromatsuihiromatsui
committed
Split routine for matrix products at bounaries
1 parent 5f2d838 commit 7846b1e

File tree

9 files changed

+588
-302
lines changed

9 files changed

+588
-302
lines changed

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

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -47,13 +47,15 @@ module const_r_mat_4_vector_sph
4747
implicit none
4848
!
4949
character(len=kchara), parameter, private &
50-
& :: vt_evo_name = 'toroidal_velocity_evolution'
50+
& :: vt_evo_name = 'toroidal_velocity_evolution'
5151
character(len=kchara), parameter, private &
52-
& :: wt_evo_name = 'toroidal_vorticity_evolution'
52+
& :: wt_evo_name = 'toroidal_vorticity_evolution'
5353
character(len=kchara), parameter, private &
54-
& :: vp_evo_name = 'poloidal_velocity_evolution'
54+
& :: wt_poison_name = 'toroidal_vorticity_Poisson'
5555
character(len=kchara), parameter, private &
56-
& :: vsp_evo_name = 'velocity_pressure_evolution'
56+
& :: vp_evo_name = 'poloidal_velocity_evolution'
57+
character(len=kchara), parameter, private &
58+
& :: vsp_evo_name = 'velocity_pressure_evolution'
5759
!
5860
! -----------------------------------------------------------------------
5961
!
@@ -65,10 +67,12 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, &
6567
& fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, &
6668
& band_vs_poisson, band_vp_evo, band_wt_evo)
6769
!
70+
use calypso_mpi
6871
use m_ludcmp_band
6972
use select_sph_r_mat_vort_BC
7073
use center_sph_matrices
7174
use mat_product_3band_mul
75+
use mat_product_3band_mul_bc
7276
!
7377
type(sph_rj_grid), intent(in) :: sph_rj
7478
type(fdm_matrices), intent(in) :: r_2nd
@@ -88,6 +92,7 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, &
8892
real(kind = kreal) :: coef_dvt
8993
!
9094
!
95+
band_vs_poisson%mat_name = wt_poison_name
9196
band_wt_evo%mat_name = wt_evo_name
9297
band_vp_evo%mat_name = vp_evo_name
9398
!
@@ -129,12 +134,21 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, &
129134
! Boundary condition for CMB
130135
call sel_sph_r_mat_vort_2step_CMB(sph_rj, sph_bc_U, bc_fdms_U, &
131136
& g_sph_rj, coef_dvt, band_vs_poisson, band_wt_evo)
132-
!
133137
!
134138
call cal_mat_product_3band_mul &
135139
& (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), &
136140
& sph_bc_U%kr_in, sph_bc_U%kr_out, band_wt_evo%mat, &
137141
& band_vs_poisson%mat, band_vp_evo%mat)
142+
call cal_vp_evo_mat_product_bc &
143+
& (sph_bc_U, sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), &
144+
& band_wt_evo%mat, band_vs_poisson%mat, band_vp_evo%mat)
145+
!
146+
! call check_specific_radial_band_mat(my_rank, (100+my_rank), 2, 0,&
147+
! & sph_rj, band_vs_poisson)
148+
! call check_specific_radial_band_mat(my_rank, (200+my_rank), 2, 0,&
149+
! & sph_rj, band_wt_evo)
150+
! call check_specific_radial_band_mat(my_rank, (300+my_rank), 2, 0,&
151+
! & sph_rj, band_vp_evo)
138152
!
139153
call ludcmp_5band_mul_t &
140154
& (np_smp, sph_rj%istack_rj_j_smp, band_vp_evo)

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

Lines changed: 186 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
!>@file mat_product_3band_mul_bc.f90
2+
!!@brief module mat_product_3band_mul_bc
3+
!!
4+
!!@author H. Matsui
5+
!!@date Programmed on 2007
6+
!!@n Modified on May, 2013
7+
!!@n Modified on July, 2025
8+
!
9+
!>@brief Take product of band matrices at boudaries
10+
!!
11+
!!@verbatim
12+
!!---------------------------------------------------------------------
13+
!! subroutine cal_vp_evo_mat_product_bc(sph_bc_U, n, mcomp, &
14+
!! & a_left, a_right, a_prod)
15+
!! type(sph_boundary_type), intent(in) :: sph_bc_U
16+
!! integer(kind = kint), intent(in) :: n, mcomp
17+
!! real(kind = kreal), intent(in) :: a_left(3,n,mcomp)
18+
!! real(kind = kreal), intent(in) :: a_right(3,n,mcomp)
19+
!! real(kind = kreal), intent(inout) :: a_prod(5,n,mcomp)
20+
!! Evaluate product of 3-band matrices at boudaries
21+
!! (a_prod) = (a_left)(a_right)
22+
!!
23+
!! subroutine cal_mat_prod_3b5b_mat_bc(n, mcomp, kr_st, kr_ed, &
24+
!! & a3_left, a5_right, a7_prod)
25+
!! integer(kind = kint), intent(in) :: kr_st, kr_ed
26+
!! integer(kind = kint), intent(in) :: n, mcomp
27+
!! real(kind = kreal), intent(in) :: a3_left(3,n,mcomp)
28+
!! real(kind = kreal), intent(in) :: a5_right(5,n,mcomp)
29+
!! real(kind = kreal), intent(inout) :: a7_prod(7,n,mcomp)
30+
!! Evaluate product of 3-band and 5-band matrix at boudaries
31+
!! (a7_prod) = (a3_left)(a5_right)
32+
!!
33+
!! Format of 3-band matrix
34+
!! | a(2,1) a(1,2) ........ 0 0 |
35+
!! | a(3,1) a(2,2) ........ . . |
36+
!! | 0 a(3,2) ........ . . |
37+
!! a(i,j) = | . 0 ........ 0 . |
38+
!! | ...... a(3,k-1) a(2,k) a(1,k+1) .......... |
39+
!! | . . ........ a(1,N-2) 0 |
40+
!! | . . ........ a(2,N-2) a(1,N-1) |
41+
!! | 0 0 ........ a(3,N-2) a(2,N-1) |
42+
!!
43+
!! Arbitorary band matrix
44+
!! band_a(i-j+iband+1,j) = a(i,j)
45+
!! band_a(k,j) = a(k+j-iband-1,j)
46+
!! 3-band matrix
47+
!! band_a(i-j+2,j) = a(i,j)
48+
!! band_a(k,j) = a(k+j-2,j)
49+
!! 5-band matrix
50+
!! band_lu(i-j+3,j) = a(i,j)
51+
!! band_lu(k,j) = a(k+j-3,j)
52+
!! 7-band matrix
53+
!! band_lu(i-j+4,j) = a(i,j)
54+
!! band_lu(k,j) = a(k+j-4,j)
55+
!! 9-band matrix
56+
!! band_lu(i-j+5,j) = a(i,j)
57+
!! band_lu(k,j) = a(k+j-5,j)
58+
!!---------------------------------------------------------------------
59+
!!@endverbatim
60+
!!
61+
!!@n @param n size of matrix
62+
!!@n @param ncomp number of matrix
63+
!!@n @param kr_st start address
64+
!!@n @param kr_ed end address
65+
!!
66+
!!@n @param a_left(3,n,mcomp) input matrix on left (3-band)
67+
!!@n @param a_right(3,n,mcomp) input matrix on right (3-band)
68+
!!@n @param a3_left(3,n,mcomp) input matrix on left (3-band)
69+
!!@n @param a5_right(5,n,mcomp) input matrix on left (5-band)
70+
!!@n @param a_prod(5,n,mcomp) produced matrix (5-band)
71+
!!@n @param a7_prod(7,n,mcomp) produced matrix (7-band)
72+
!
73+
!
74+
module mat_product_3band_mul_bc
75+
!
76+
use m_precision
77+
use m_constants
78+
!
79+
implicit none
80+
!
81+
! -----------------------------------------------------------------------
82+
!
83+
contains
84+
!
85+
! -----------------------------------------------------------------------
86+
!
87+
subroutine cal_vp_evo_mat_product_bc(sph_bc_U, n, mcomp, &
88+
& a_left, a_right, a_prod)
89+
!
90+
use t_boundary_params_sph_MHD
91+
!
92+
type(sph_boundary_type), intent(in) :: sph_bc_U
93+
integer(kind = kint), intent(in) :: n, mcomp
94+
real(kind = kreal), intent(in) :: a_left(3,n,mcomp)
95+
real(kind = kreal), intent(in) :: a_right(3,n,mcomp)
96+
!
97+
real(kind = kreal), intent(inout) :: a_prod(5,n,mcomp)
98+
!
99+
integer(kind = kint) :: k, j
100+
!
101+
!$omp parallel
102+
if(sph_bc_U%iflag_icb .eq. iflag_sph_fill_center) then
103+
k = 1
104+
!$omp do private(j)
105+
do j = 1, mcomp
106+
! a_prod(5,k-2,j) = zero
107+
! a_prod(4,k-1,j) = zero
108+
a_prod(3,k, j) = a_left(2,k ,j) * a_right(2,k ,j) &
109+
& + a_left(1,k+1,j) * a_right(3,k ,j)
110+
a_prod(2,k+1,j) = a_left(2,k ,j) * a_right(1,k+1,j) &
111+
& + a_left(1,k+1,j) * a_right(2,k+1,j)
112+
a_prod(1,k+2,j) = a_left(1,k+1,j) * a_right(1,k+2,j)
113+
end do
114+
!$omp end do nowait
115+
else
116+
k = sph_bc_U%kr_in
117+
!$omp do private (j)
118+
do j = 1, mcomp
119+
a_prod(3,k, j) = one
120+
a_prod(2,k+1,j) = zero
121+
a_prod(1,k+2,j) = zero
122+
end do
123+
!$omp end do nowait
124+
end if
125+
!
126+
!
127+
k = sph_bc_U%kr_out
128+
!$omp do private(j)
129+
do j = 1, mcomp
130+
a_prod(5,k-2,j) = zero
131+
a_prod(4,k-1,j) = zero
132+
a_prod(3,k, j) = one
133+
end do
134+
!$omp end do
135+
!$omp end parallel
136+
!
137+
end subroutine cal_vp_evo_mat_product_bc
138+
!
139+
! -----------------------------------------------------------------------
140+
!
141+
subroutine cal_mat_prod_3b5b_mat_bc(n, mcomp, kr_st, kr_ed, &
142+
& a3_left, a5_right, a7_prod)
143+
!
144+
integer(kind = kint), intent(in) :: kr_st, kr_ed
145+
integer(kind = kint), intent(in) :: n, mcomp
146+
real(kind = kreal), intent(in) :: a3_left(3,n,mcomp)
147+
real(kind = kreal), intent(in) :: a5_right(5,n,mcomp)
148+
!
149+
real(kind = kreal), intent(inout) :: a7_prod(7,n,mcomp)
150+
!
151+
integer(kind = kint) :: k, j
152+
!
153+
!
154+
!$omp parallel
155+
k = kr_st
156+
!$omp do private(j)
157+
do j = 1, mcomp
158+
a7_prod(4,k, j) = a3_left(2,k ,j) * a5_right(3,k ,j) &
159+
& + a3_left(1,k+1,j) * a5_right(4,k ,j)
160+
a7_prod(3,k+1,j) = a3_left(2,k ,j) * a5_right(2,k+1,j) &
161+
& + a3_left(1,k+1,j) * a5_right(3,k+1,j)
162+
a7_prod(2,k+2,j) = a3_left(2,k ,j) * a5_right(1,k+2,j) &
163+
& + a3_left(1,k+1,j) * a5_right(2,k+2,j)
164+
a7_prod(1,k+3,j) = a3_left(1,k+1,j) * a5_right(1,k+3,j)
165+
end do
166+
!$omp end do nowait
167+
!
168+
k = kr_ed
169+
!$omp do private(j)
170+
do j = 1, mcomp
171+
a7_prod(7,k-3,j) = a3_left(3,k-1,j) * a5_right(5,k-3,j)
172+
a7_prod(6,k-2,j) = a3_left(3,k-1,j) * a5_right(4,k-2,j) &
173+
& + a3_left(2,k ,j) * a5_right(5,k-2,j)
174+
a7_prod(5,k-1,j) = a3_left(3,k-1,j) * a5_right(3,k-1,j) &
175+
& + a3_left(2,k ,j) * a5_right(4,k-1,j)
176+
a7_prod(4,k, j) = a3_left(3,k-1,j) * a5_right(2,k ,j) &
177+
& + a3_left(2,k ,j) * a5_right(3,k ,j)
178+
end do
179+
!$omp end do
180+
!$omp end parallel
181+
!
182+
end subroutine cal_mat_prod_3b5b_mat_bc
183+
!
184+
! -----------------------------------------------------------------------
185+
!
186+
end module mat_product_3band_mul_bc

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -93,7 +93,7 @@ end subroutine solve_velo_by_vort_sph_crank
9393
subroutine solve_pressure_by_div_v(sph_rj, band_p_poisson, &
9494
& is_press, n_point, ntot_phys_rj, d_rj)
9595
!
96-
use check_sph_radial_mat
96+
use check_single_radial_mat
9797
!
9898
type(sph_rj_grid), intent(in) :: sph_rj
9999
type(band_matrices_type), intent(in) :: band_p_poisson

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -189,7 +189,7 @@ end subroutine lubksb_7band_ctr
189189
subroutine check_center_band_matrix(id_file, sph_rj, smat)
190190
!
191191
use t_spheric_rj_data
192-
use check_sph_radial_mat
192+
use check_single_radial_mat
193193
!
194194
integer(kind = kint), intent(in) :: id_file
195195
type(sph_rj_grid), intent(in) :: sph_rj

0 commit comments

Comments
 (0)