[cig-commits] [commit] Hiro_latest: Use structures for boundary condition in matrix construction (3bca14d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Nov 18 16:22:07 PST 2013
Repository : ssh://geoshell/calypso
On branch : Hiro_latest
Link : https://github.com/geodynamics/calypso/compare/93e9f8f974c7a247c8f02e54ec18de063f86c8fb...3c548304673360ddedd7d68c8095b3fb74a2b9ce
>---------------------------------------------------------------
commit 3bca14d61bcd9159efebc8650dbb22d26ffe81f0
Author: Hiroaki Matsui <h_kemono at mac.com>
Date: Sun Nov 17 09:33:26 2013 -0800
Use structures for boundary condition in matrix construction
>---------------------------------------------------------------
3bca14d61bcd9159efebc8650dbb22d26ffe81f0
.../MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90 | 2 +-
.../MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90 | 18 +-
.../MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90 | 43 +++--
.../MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90 | 9 +
.../MHD_src/sph_MHD/set_sph_magne_mat_bc.f90 | 164 +++++++++++-------
.../MHD_src/sph_MHD/set_sph_mom_mat_bc.f90 | 186 ++++++++++++---------
.../MHD_src/sph_MHD/set_sph_scalar_mat_bc.f90 | 82 +++++----
7 files changed, 301 insertions(+), 203 deletions(-)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90
index 982b96e..d2d7e0b 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_sph_exp_nod_cmb_ins.f90
@@ -25,7 +25,7 @@
!!@n @param kr_out Radial ID for outer boundary
!!@n @param r_CMB(0:2) Radius at CMB
!!@n @param fdm2_fix_fld_CMB(0:2,3)
-!!! Matrix to evaluate radial derivative at CMB with fiexed field
+!! Matrix to evaluate radial derivative at CMB with fiexed field
!!@n @param fdm2_fix_dr_CMB(-1:1,3)
!! Matrix to evaluate field at CMB with fiexed radial derivative
!!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90
index 23ea119..b1031e1 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90
@@ -49,18 +49,20 @@
!
if (sph_bc_T%iflag_icb .eq. iflag_fixed_flux) then
call set_fix_flux_icb_rmat_sph(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_T%kr_in, sph_bc_T%r_ICB, sph_bc_T%fdm2_fix_dr_ICB, &
& coef_imp_t, coef_d_temp, temp_evo_mat)
else
call set_fix_scalar_icb_rmat_sph(nidx_rj(1), nidx_rj(2), &
- & temp_evo_mat)
+ & sph_bc_T%kr_in, temp_evo_mat)
end if
!
if (sph_bc_T%iflag_cmb .eq. iflag_fixed_flux) then
call set_fix_flux_cmb_rmat_sph(nidx_rj(1), nidx_rj(2), &
- & coef_imp_t, coef_d_temp, temp_evo_mat)
+ & sph_bc_T%kr_out, sph_bc_T%r_CMB, sph_bc_T%fdm2_fix_dr_CMB, &
+ & coef_imp_t, coef_d_temp, temp_evo_mat)
else
call set_fix_scalar_cmb_rmat_sph(nidx_rj(1), nidx_rj(2), &
- & temp_evo_mat)
+ & sph_bc_C%kr_out, temp_evo_mat)
end if
!
!$omp parallel do private(jst,jed,j)
@@ -95,18 +97,20 @@
!
if (sph_bc_C%iflag_icb .eq. iflag_fixed_flux) then
call set_fix_flux_icb_rmat_sph(nidx_rj(1), nidx_rj(2), &
- & coef_imp_c, coef_d_light, composit_evo_mat)
+ & sph_bc_C%kr_in, sph_bc_C%r_ICB, sph_bc_C%fdm2_fix_dr_ICB, &
+ & coef_imp_c, coef_d_light, composit_evo_mat)
else
call set_fix_scalar_icb_rmat_sph(nidx_rj(1), nidx_rj(2), &
- & composit_evo_mat)
+ & sph_bc_C%kr_in, composit_evo_mat)
end if
!
if (sph_bc_C%iflag_cmb .eq. iflag_fixed_flux) then
call set_fix_flux_cmb_rmat_sph(nidx_rj(1), nidx_rj(2), &
- & coef_imp_c, coef_d_light, composit_evo_mat)
+ & sph_bc_C%kr_out, sph_bc_C%r_CMB, sph_bc_C%fdm2_fix_dr_CMB, &
+ & coef_imp_c, coef_d_light, composit_evo_mat)
else
call set_fix_scalar_cmb_rmat_sph(nidx_rj(1), nidx_rj(2), &
- & composit_evo_mat)
+ & sph_bc_C%kr_out, composit_evo_mat)
end if
!
!$omp parallel do private(jst,jed,j)
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90
index 8cc4a37..2a67a6b 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_vector_sph.f90
@@ -58,8 +58,12 @@
!
! Boundary condition for ICB
!
- call set_icb_wt_sph_evo_mat
- call set_icb_p_sph_poisson_mat
+ call set_icb_wt_sph_evo_mat(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & wt_evo_mat)
+ call set_icb_p_sph_poisson_mat(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_U%kr_in, sph_bc_U%r_ICB, sph_bc_U%fdm2_fix_dr_ICB, &
+ & p_poisson_mat)
!
if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then
call set_free_slip_icb_vt_sph_mat
@@ -77,8 +81,12 @@
!
! Boundary condition for CMB
!
- call set_cmb_wt_sph_evo_mat
- call set_cmb_p_sph_poisson_mat
+ call set_cmb_wt_sph_evo_mat(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_U%kr_out, sph_bc_U%r_CMB, sph_bc_U%fdm2_fix_dr_CMB, &
+ & wt_evo_mat)
+ call set_cmb_p_sph_poisson_mat(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_U%kr_out, sph_bc_U%r_CMB, sph_bc_U%fdm2_fix_dr_CMB, &
+ & p_poisson_mat)
!
if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then
call set_free_slip_cmb_vt_sph_mat
@@ -135,32 +143,39 @@
use m_boundary_params_sph_MHD
use set_sph_magne_mat_bc
!
- integer(kind = kint) :: kr_in, ip, jst, jed, j
+ integer(kind = kint) :: ip, jst, jed, j
integer(kind = kint) :: ierr
!
!
if(sph_bc_B%iflag_icb .eq. iflag_sph_fill_center) then
- kr_in = ione
call set_magne_center_rmat_sph
else if(sph_bc_B%iflag_icb .eq. iflag_radial_magne) then
- kr_in = nlayer_ICB
- call set_qvacume_magne_icb_rmat_sph
+ call no_r_poynting_cmb_rmat_sph(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_B%kr_in, sph_bc_B%r_ICB, sph_bc_B%fdm2_fix_dr_ICB, &
+ & bs_evo_mat, bt_evo_mat)
else
- kr_in = nlayer_ICB
- call set_ins_magne_icb_rmat_sph
+ call set_ins_magne_icb_rmat_sph(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_B%kr_in, sph_bc_B%r_ICB, sph_bc_B%fdm2_fix_dr_ICB, &
+ & bs_evo_mat, bt_evo_mat)
end if
!
if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then
- call set_qvacume_magne_cmb_rmat_sph
+ call no_r_poynting_cmb_rmat_sph(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_B%kr_out, sph_bc_B%r_CMB, sph_bc_B%fdm2_fix_dr_CMB, &
+ & bs_evo_mat, bt_evo_mat)
else
- call set_ins_magne_cmb_rmat_sph
+ call set_ins_magne_cmb_rmat_sph(nidx_rj(1), nidx_rj(2), &
+ & sph_bc_B%kr_out, sph_bc_B%r_CMB, sph_bc_B%fdm2_fix_dr_CMB, &
+ & bs_evo_mat, bt_evo_mat)
end if
!
!$omp parallel
call set_radial_vect_evo_mat_sph(nidx_rj(1), nidx_rj(2), &
- & kr_in, nlayer_CMB, coef_imp_b, coef_d_magne, bs_evo_mat)
+ & sph_bc_B%kr_in, sph_bc_B%kr_out, coef_imp_b, coef_d_magne, &
+ & bs_evo_mat)
call set_radial_vect_evo_mat_sph(nidx_rj(1), nidx_rj(2), &
- & kr_in, nlayer_CMB, coef_imp_b, coef_d_magne, bt_evo_mat)
+ & sph_bc_B%kr_in, sph_bc_B%kr_out, coef_imp_b, coef_d_magne, &
+ & bt_evo_mat)
!$omp end parallel
!
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90
index e360dc3..8f63b5d 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_exp_rigid_CMB.f90
@@ -21,6 +21,15 @@
!! & fdm2_fix_fld_CMB, coef_d, is_fld, is_diffuse)
!!@endverbatim
!!
+!!@n @param jmax Number of modes for spherical harmonics @f$L*(L+2)@f$
+!!@n @param kr_out Radial ID for outer boundary
+!!@n @param r_CMB(0:2) Radius at CMB
+!!
+!!@n @param fdm2_fix_fld_CMB(0:2,3)
+!! Matrix to evaluate radial derivative at CMB with fiexed field
+!!@n @param fdm2_fix_dr_CMB(-1:1,3)
+!! Matrix to evaluate field at CMB with fiexed radial derivative
+!!
!!@n @param coef_d Coefficient for diffusion term
!!@n @param is_fld Address of poloidal velocity in d_rj
!!@n @param is_rot Address of poloidal vorticity in d_rj
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_magne_mat_bc.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_magne_mat_bc.f90
index 3ebea9a..21d0ea9 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_magne_mat_bc.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_magne_mat_bc.f90
@@ -11,14 +11,33 @@
!! subroutine set_magne_center_rmat_sph
!!
!! Boundary condition to connect potential field
-!! subroutine set_ins_magne_icb_rmat_sph
-!! subroutine set_ins_magne_cmb_rmat_sph
+!! subroutine set_ins_magne_icb_rmat_sph(nri, jmax, kr_in, r_ICB, &
+!! & fdm2_fix_dr_ICB, bs_evo_mat, bt_evo_mat)
+!! subroutine set_ins_magne_cmb_rmat_sph(nri, jmax, kr_out, r_CMB, &
+!! & fdm2_fix_dr_CMB, bs_evo_mat, bt_evo_mat)
!!
!! Boundary condition for radial magnetic field
-!! subroutine set_qvacume_magne_icb_rmat_sph
-!! subroutine set_qvacume_magne_cmb_rmat_sph
+!! subroutine no_r_poynting_icb_rmat_sph(nri, jmax, kr_in, r_ICB, &
+!! & fdm2_fix_dr_ICB, bs_evo_mat, bt_evo_mat)
+!! subroutine no_r_poynting_cmb_rmat_sph(nri, jmax, kr_out, r_CMB, &
+!! & fdm2_fix_dr_CMB, bs_evo_mat, bt_evo_mat)
!!@endverbatim
!
+!!@n @param jmax Number of local spherical harmonics mode
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param kr_out Radial ID for outer boundary
+!!@n @param r_ICB(0:2) Radius at ICB
+!!@n @param r_CMB(0:2) Radius at CMB
+!!@n @param fdm2_fix_dr_ICB(-1:1,3)
+!! Matrix to evaluate field at ICB with fiexed radial derivative
+!!@n @param fdm2_fix_dr_CMB(-1:1,3)
+!! Matrix to evaluate field at CMB with fiexed radial derivative
+!!
+!!@n @param bs_evo_mat(3,nri,jmax) 3-band matrix for evolution of
+!! poloidal magnetic field
+!!@n @param bt_evo_mat(3,nri,jmax) 3-band matrix for Poisson equation
+!! toroidal magnetic field
+!
module set_sph_magne_mat_bc
!
use m_precision
@@ -26,9 +45,7 @@
use m_constants
use m_t_int_parameter
use m_physical_property
- use m_spheric_parameter
use m_schmidt_poly_on_rtm
- use m_radial_matrices_sph
!
implicit none
!
@@ -40,6 +57,8 @@
!
subroutine set_magne_center_rmat_sph
!
+ use m_spheric_parameter
+ use m_radial_matrices_sph
use m_fdm_coefs
!
integer(kind = kint) :: j
@@ -63,53 +82,60 @@
!
! -----------------------------------------------------------------------
!
- subroutine set_ins_magne_icb_rmat_sph
+ subroutine set_ins_magne_icb_rmat_sph(nri, jmax, kr_in, r_ICB, &
+ & fdm2_fix_dr_ICB, bs_evo_mat, bt_evo_mat)
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_in
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
- use m_coef_fdm_fixed_ICB
+ real(kind = kreal), intent(inout) :: bs_evo_mat(3,nri,jmax)
+ real(kind = kreal), intent(inout) :: bt_evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- bs_evo_mat(2,nlayer_ICB, j) &
- & = one + coef_imp_b*dt*coef_d_magne &
- & * ( -coef_fdm_fix_dr_ICB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2) &
- & - (g_sph_rj(j,1)+one) &
- & * ar_1d_rj(nlayer_ICB,1) &
- & * coef_fdm_fix_dr_ICB_2(-1,3) )
- bs_evo_mat(1,nlayer_ICB+1,j) &
- & = - coef_imp_b*dt*coef_d_magne &
- & * coef_fdm_fix_dr_ICB_2( 1,3)
-!
- bt_evo_mat(2,nlayer_ICB, j) = one
- bt_evo_mat(1,nlayer_ICB+1,j) = zero
+ do j = 1, jmax
+ bs_evo_mat(2,kr_in, j) = one + coef_imp_b*dt*coef_d_magne &
+ & * ( -fdm2_fix_dr_ICB( 0,3) &
+ & + g_sph_rj(j,3)*r_ICB(2) &
+ & - (g_sph_rj(j,1)+one) * r_ICB(1) &
+ & * fdm2_fix_dr_ICB(-1,3) )
+ bs_evo_mat(1,kr_in+1,j) = - coef_imp_b*dt*coef_d_magne &
+ & * fdm2_fix_dr_ICB( 1,3)
+!
+ bt_evo_mat(2,kr_in, j) = one
+ bt_evo_mat(1,kr_in+1,j) = zero
end do
!
end subroutine set_ins_magne_icb_rmat_sph
!
! -----------------------------------------------------------------------
!
- subroutine set_ins_magne_cmb_rmat_sph
+ subroutine set_ins_magne_cmb_rmat_sph(nri, jmax, kr_out, r_CMB, &
+ & fdm2_fix_dr_CMB, bs_evo_mat, bt_evo_mat)
!
- use m_coef_fdm_fixed_CMB
+ integer(kind = kint), intent(in) :: nri, jmax, kr_out
+ real(kind = kreal), intent(in) :: r_CMB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
+!
+ real(kind = kreal), intent(inout) :: bs_evo_mat(3,nri,jmax)
+ real(kind = kreal), intent(inout) :: bt_evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- bs_evo_mat(3,nlayer_CMB-1,j) &
- & = - coef_imp_b*dt*coef_d_magne &
- & * coef_fdm_fix_dr_CMB_2(-1,3)
- bs_evo_mat(2,nlayer_CMB, j) &
- & = one + coef_imp_b*dt*coef_d_magne &
- & * ( -coef_fdm_fix_dr_CMB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2) &
- & + g_sph_rj(j,1)*ar_1d_rj(nlayer_CMB,1) &
- & * coef_fdm_fix_dr_CMB_2( 1,3) )
-!
- bt_evo_mat(3,nlayer_CMB-1,j) = zero
- bt_evo_mat(2,nlayer_CMB, j) = one
+ do j = 1, jmax
+ bs_evo_mat(3,kr_out-1,j) = - coef_imp_b*dt*coef_d_magne &
+ & * fdm2_fix_dr_CMB(-1,3)
+ bs_evo_mat(2,kr_out, j) = one + coef_imp_b*dt*coef_d_magne &
+ & * ( -fdm2_fix_dr_CMB( 0,3) &
+ & + g_sph_rj(j,3)*r_CMB(2) &
+ & + g_sph_rj(j,1)*r_CMB(1) &
+ & * fdm2_fix_dr_CMB( 1,3) )
+!
+ bt_evo_mat(3,kr_out-1,j) = zero
+ bt_evo_mat(2,kr_out, j) = one
end do
!
end subroutine set_ins_magne_cmb_rmat_sph
@@ -117,51 +143,59 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine set_qvacume_magne_icb_rmat_sph
+ subroutine no_r_poynting_icb_rmat_sph(nri, jmax, kr_in, r_ICB, &
+ & fdm2_fix_dr_ICB, bs_evo_mat, bt_evo_mat)
!
- use m_coef_fdm_fixed_ICB
+ integer(kind = kint), intent(in) :: nri, jmax, kr_in
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
+!
+ real(kind = kreal), intent(inout) :: bs_evo_mat(3,nri,jmax)
+ real(kind = kreal), intent(inout) :: bt_evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- bs_evo_mat(2,nlayer_ICB, j) &
- & = one + coef_imp_b*dt*coef_d_magne &
- & * ( -coef_fdm_fix_dr_ICB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2) )
- bs_evo_mat(1,nlayer_ICB+1,j) &
- & = - coef_imp_b*dt*coef_d_magne &
- & * coef_fdm_fix_dr_ICB_2( 1,3)
-!
- bt_evo_mat(2,nlayer_ICB, j) = one
- bt_evo_mat(1,nlayer_ICB+1,j) = zero
+ do j = 1, jmax
+ bs_evo_mat(2,kr_in, j) = one + coef_imp_b*dt*coef_d_magne &
+ & * ( -fdm2_fix_dr_ICB( 0,3) &
+ & + g_sph_rj(j,3)*r_ICB(2) )
+ bs_evo_mat(1,kr_in+1,j) = - coef_imp_b*dt*coef_d_magne &
+ & * fdm2_fix_dr_ICB( 1,3)
+!
+ bt_evo_mat(2,kr_in, j) = one
+ bt_evo_mat(1,kr_in+1,j) = zero
end do
!
- end subroutine set_qvacume_magne_icb_rmat_sph
+ end subroutine no_r_poynting_icb_rmat_sph
!
! -----------------------------------------------------------------------
!
- subroutine set_qvacume_magne_cmb_rmat_sph
+ subroutine no_r_poynting_cmb_rmat_sph(nri, jmax, kr_out, r_CMB, &
+ & fdm2_fix_dr_CMB, bs_evo_mat, bt_evo_mat)
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_out
+ real(kind = kreal), intent(in) :: r_CMB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
- use m_coef_fdm_fixed_CMB
+ real(kind = kreal), intent(inout) :: bs_evo_mat(3,nri,jmax)
+ real(kind = kreal), intent(inout) :: bt_evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- bs_evo_mat(3,nlayer_CMB-1,j) &
- & = - coef_imp_b*dt*coef_d_magne &
- & * coef_fdm_fix_dr_CMB_2(-1,3)
- bs_evo_mat(2,nlayer_CMB, j) &
- & = one + coef_imp_b*dt*coef_d_magne &
- & * ( -coef_fdm_fix_dr_CMB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2) )
-!
- bt_evo_mat(3,nlayer_CMB-1,j) = zero
- bt_evo_mat(2,nlayer_CMB, j) = one
+ do j = 1, jmax
+ bs_evo_mat(3,kr_out-1,j) = - coef_imp_b*dt*coef_d_magne &
+ & * fdm2_fix_dr_CMB(-1,3)
+ bs_evo_mat(2,kr_out, j) = one + coef_imp_b*dt*coef_d_magne &
+ & * ( -fdm2_fix_dr_CMB( 0,3) &
+ & + g_sph_rj(j,3)*r_CMB(2) )
+!
+ bt_evo_mat(3,kr_out-1,j) = zero
+ bt_evo_mat(2,kr_out, j) = one
end do
!
- end subroutine set_qvacume_magne_cmb_rmat_sph
+ end subroutine no_r_poynting_cmb_rmat_sph
!
! -----------------------------------------------------------------------
!
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_mom_mat_bc.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_mom_mat_bc.f90
index 59c3791..4ba678d 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_mom_mat_bc.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_mom_mat_bc.f90
@@ -7,14 +7,35 @@
!>@brief Construct matrix for fixed velocity at boundaries
!!
!!@verbatim
-!! subroutine set_icb_wt_sph_evo_mat
-!! subroutine set_icb_vp_sph_poisson_mat
-!! subroutine set_icb_p_sph_poisson_mat
+!! subroutine set_icb_wt_sph_evo_mat(nri, jmax, kr_in, r_ICB, &
+!! & fdm2_fix_dr_ICB, wt_evo_mat)
+!! subroutine set_icb_vp_sph_poisson_mat(nri, jmax, kr_in, &
+!! & poisson_mat)
+!! subroutine set_icb_p_sph_poisson_mat(nri, jmax, kr_in, r_ICB, &
+!! & fdm2_fix_dr_ICB, poisson_mat)
!!
-!! subroutine set_cmb_wt_sph_evo_mat
-!! subroutine set_cmb_vp_sph_poisson_mat
-!! subroutine set_cmb_p_sph_poisson_mat
+!! subroutine set_cmb_wt_sph_evo_mat(nri, jmax, kr_out, r_CMB, &
+!! & fdm2_fix_dr_CMB, wt_evo_mat)
+!! subroutine set_cmb_vp_sph_poisson_mat(nri, jmax, kr_out, &
+!! & poisson_mat)
+!! subroutine set_cmb_p_sph_poisson_mat(nri, jmax, kr_out, r_CMB, &
+!! & fdm2_fix_dr_CMB, poisson_mat)
!!@endverbatim
+!!
+!!@n @param nri Number of radial points
+!!@n @param jmax Number of local spherical harmonics mode
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param kr_out Radial ID for outer boundary
+!!@n @param r_ICB(0:2) Radius at ICB
+!!@n @param r_CMB(0:2) Radius at CMB
+!!@n @param fdm2_fix_dr_ICB(-1:1,3)
+!! Matrix to evaluate field at ICB with fiexed radial derivative
+!!@n @param fdm2_fix_dr_CMB(-1:1,3)
+!! Matrix to evaluate field at CMB with fiexed radial derivative
+!!
+!!@n @param wt_evo_mat(3,nri,jmax) 3-band matrix for evolution of
+!! toroidal vorticity
+!!@n @param poisson_mat(3,nri,jmax) 3-band matrix for Poisson equation
!
module set_sph_mom_mat_bc
!
@@ -23,10 +44,7 @@
use m_constants
use m_t_int_parameter
use m_physical_property
- use m_spheric_parameter
use m_schmidt_poly_on_rtm
- use m_radial_matrices_sph
- use m_fdm_coefs
!
implicit none
!
@@ -36,77 +54,75 @@
!
! -----------------------------------------------------------------------
!
- subroutine set_icb_wt_sph_evo_mat
+ subroutine set_icb_wt_sph_evo_mat(nri, jmax, kr_in, r_ICB, &
+ & fdm2_fix_dr_ICB, wt_evo_mat)
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_in
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
- use m_coef_fdm_fixed_ICB
+ real(kind = kreal), intent(inout) :: wt_evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
-! wt_evo_mat(2,nlayer_ICB, j) &
+ do j = 1, jmax
+! wt_evo_mat(2,kr_in, j) &
! & = one + coef_imp_v*dt*coef_d_velo &
-! & * g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2)
-! wt_evo_mat(1,nlayer_ICB+1,j) = zero
+! & * g_sph_rj(j,3)*r_ICB(2)
+! wt_evo_mat(1,kr_in+1,j) = zero
!
- wt_evo_mat(2,nlayer_ICB, j) &
+ wt_evo_mat(2,kr_in, j) &
& = one + coef_imp_v*dt*coef_d_velo &
- & * ( -coef_fdm_fix_dr_ICB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2))
- wt_evo_mat(1,nlayer_ICB+1,j) &
+ & * ( -fdm2_fix_dr_ICB( 0,3) &
+ & + g_sph_rj(j,3)*r_ICB(2))
+ wt_evo_mat(1,kr_in+1,j) &
& = - coef_imp_v*dt*coef_d_velo &
- & * coef_fdm_fix_dr_ICB_2( 1,3)
+ & * fdm2_fix_dr_ICB( 1,3)
end do
!
end subroutine set_icb_wt_sph_evo_mat
!
! -----------------------------------------------------------------------
!
- subroutine set_icb_vp_sph_poisson_mat
+ subroutine set_icb_vp_sph_poisson_mat(nri, jmax, kr_in, &
+ & poisson_mat)
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_in
+ real(kind = kreal), intent(inout) :: poisson_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- vs_poisson_mat(2,nlayer_ICB, j) = one
- vs_poisson_mat(1,nlayer_ICB+1,j) = zero
+ do j = 1, jmax
+ poisson_mat(2,kr_in, j) = one
+ poisson_mat(1,kr_in+1,j) = zero
end do
!
end subroutine set_icb_vp_sph_poisson_mat
!
! -----------------------------------------------------------------------
!
- subroutine set_icb_p_sph_poisson_mat
+ subroutine set_icb_p_sph_poisson_mat(nri, jmax, kr_in, &
+ & r_ICB, fdm2_fix_dr_ICB, poisson_mat)
+!
+ use m_fdm_coefs
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_in
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
- use m_coef_fdm_fixed_ICB
+ real(kind = kreal), intent(inout) :: poisson_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- p_poisson_mat(2,nlayer_ICB, j) = coef_press &
- & * (coef_fdm_fix_dr_ICB_2( 0,3) &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2) &
- & + two * ar_1d_rj(nlayer_ICB,1) &
- & * coef_fdm_fix_dr_ICB_2( 0,2) )
- p_poisson_mat(1,nlayer_ICB+1,j) = coef_press &
- & * (coef_fdm_fix_dr_ICB_2( 1,3) &
- & + two * ar_1d_rj(nlayer_ICB,1) &
- & * coef_fdm_fix_dr_ICB_2( 1,2) )
-!
- p_poisson_mat(3,nlayer_ICB,j ) = coef_press &
- & * (d2nod_mat_fdm_2(nlayer_ICB+1,-1) &
- & + two*ar_1d_rj(nlayer_ICB+1,1) &
- & * d1nod_mat_fdm_2(nlayer_ICB+1,-1) )
- p_poisson_mat(2,nlayer_ICB+1,j) = coef_press &
- & * ( d2nod_mat_fdm_2(nlayer_ICB+1, 0) &
- & + two*ar_1d_rj(nlayer_ICB+1,1) &
- & * d1nod_mat_fdm_2(nlayer_ICB+1, 0) &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB+1,2) )
- p_poisson_mat(1,nlayer_ICB+2,j) = coef_press &
- & * ( d2nod_mat_fdm_2(nlayer_ICB+1, 1) &
- & + two*ar_1d_rj(nlayer_ICB+1,1) &
- & * d1nod_mat_fdm_2(nlayer_ICB+1, 1) )
+ do j = 1, jmax
+ poisson_mat(2,kr_in, j) = coef_press * (fdm2_fix_dr_ICB( 0,3) &
+ & - g_sph_rj(j,3)*r_ICB(2) &
+ & + two * r_ICB(1) * fdm2_fix_dr_ICB( 0,2) )
+ poisson_mat(1,kr_in+1,j) = coef_press * (fdm2_fix_dr_ICB( 1,3) &
+ & + two * r_ICB(1) * fdm2_fix_dr_ICB( 1,2) )
end do
!
end subroutine set_icb_p_sph_poisson_mat
@@ -114,63 +130,73 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine set_cmb_wt_sph_evo_mat
+ subroutine set_cmb_wt_sph_evo_mat(nri, jmax, kr_out, r_CMB, &
+ & fdm2_fix_dr_CMB, wt_evo_mat)
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_out
+ real(kind = kreal), intent(in) :: r_CMB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
- use m_coef_fdm_fixed_CMB
+ real(kind = kreal), intent(inout) :: wt_evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
-! wt_evo_mat(3,nlayer_CMB-1,j) = zero
-! wt_evo_mat(2,nlayer_CMB, j) &
+ do j = 1, jmax
+! wt_evo_mat(3,kr_out-1,j) = zero
+! wt_evo_mat(2,kr_out, j) &
! & = one + coef_imp_v*dt*coef_d_velo &
-! & * g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2)
-!
- wt_evo_mat(3,nlayer_CMB-1,j) &
- & = - coef_imp_v*dt*coef_d_velo &
- & * coef_fdm_fix_dr_CMB_2(-1,3)
- wt_evo_mat(2,nlayer_CMB, j) &
- & = one + coef_imp_v*dt*coef_d_velo &
- & * ( -coef_fdm_fix_dr_CMB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2))
+! & * g_sph_rj(j,3)*r_CMB(2)
+!
+ wt_evo_mat(3,kr_out-1,j) &
+ & = - coef_imp_v*dt*coef_d_velo &
+ & * fdm2_fix_dr_CMB(-1,3)
+ wt_evo_mat(2,kr_out, j) &
+ & = one + coef_imp_v*dt*coef_d_velo &
+ & * ( -fdm2_fix_dr_CMB( 0,3) &
+ & + g_sph_rj(j,3)*r_CMB(2))
end do
!
end subroutine set_cmb_wt_sph_evo_mat
!
! -----------------------------------------------------------------------
!
- subroutine set_cmb_vp_sph_poisson_mat
+ subroutine set_cmb_vp_sph_poisson_mat(nri, jmax, kr_out, &
+ & poisson_mat)
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_out
+ real(kind = kreal), intent(inout) :: poisson_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- vs_poisson_mat(3,nlayer_CMB-1,j) = zero
- vs_poisson_mat(2,nlayer_CMB, j) = one
+ do j = 1, jmax
+ poisson_mat(3,kr_out-1,j) = zero
+ poisson_mat(2,kr_out, j) = one
end do
!
end subroutine set_cmb_vp_sph_poisson_mat
!
! -----------------------------------------------------------------------
!
- subroutine set_cmb_p_sph_poisson_mat
+ subroutine set_cmb_p_sph_poisson_mat(nri, jmax, kr_out, r_CMB, &
+ & fdm2_fix_dr_CMB, poisson_mat)
+!
+ integer(kind = kint), intent(in) :: nri, jmax, kr_out
+ real(kind = kreal), intent(in) :: r_CMB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
- use m_coef_fdm_fixed_CMB
+ real(kind = kreal), intent(inout) :: poisson_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
- do j = 1, nidx_rj(2)
- p_poisson_mat(3,nlayer_CMB-1,j) = coef_press &
- & * ( coef_fdm_fix_dr_CMB_2(-1,3) &
- & + two* ar_1d_rj(nlayer_CMB,1) &
- & * coef_fdm_fix_dr_CMB_2(-1,2) )
- p_poisson_mat(2,nlayer_CMB, j) = coef_press &
- & * ( coef_fdm_fix_dr_CMB_2( 0,3) &
- & - g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2) &
- & + two * ar_1d_rj(nlayer_CMB,1) &
- & * coef_fdm_fix_dr_CMB_2( 0,2) )
+ do j = 1, jmax
+ poisson_mat(3,kr_out-1,j) = coef_press * (fdm2_fix_dr_CMB(-1,3) &
+ & + two* r_CMB(1) * fdm2_fix_dr_CMB(-1,2) )
+ poisson_mat(2,kr_out, j) = coef_press * (fdm2_fix_dr_CMB( 0,3) &
+ & - g_sph_rj(j,3)*r_CMB(2) &
+ & + two * r_CMB(1) * fdm2_fix_dr_CMB( 0,2) )
end do
!
end subroutine set_cmb_p_sph_poisson_mat
diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_scalar_mat_bc.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_scalar_mat_bc.f90
index 5d63208..f828feb 100644
--- a/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_scalar_mat_bc.f90
+++ b/src/Fortran_libraries/MHD_src/sph_MHD/set_sph_scalar_mat_bc.f90
@@ -7,21 +7,29 @@
!>@brief Construct matrix for scalar fields at boundaries
!!
!!@verbatim
-!! subroutine set_fix_scalar_icb_rmat_sph(nri, jmax, evo_mat)
-!! subroutine set_fix_flux_icb_rmat_sph(nri, jmax, &
-!! & coef_imp, coef_d, evo_mat)
+!! subroutine set_fix_scalar_icb_rmat_sph(nri, jmax, kr_in, &
+!! & evo_mat)
+!! subroutine set_fix_flux_icb_rmat_sph(nri, jmax, kr_in, r_ICB, &
+!! & fdm2_fix_dr_ICB, coef_imp, coef_d, evo_mat)
!!
-!! subroutine set_fix_scalar_cmb_rmat_sph(nri, jmax, evo_mat)
-!! subroutine set_fix_flux_cmb_rmat_sph(nri, jmax, &
-!! & coef_imp, coef_d, evo_mat)
+!! subroutine set_fix_scalar_cmb_rmat_sph(nri, jmax, kr_out, &
+!! & evo_mat)
+!! subroutine set_fix_flux_cmb_rmat_sph(nri, jmax, kr_out, r_CMB, &
+!! & fdm2_fix_dr_CMB, coef_imp, coef_d, evo_mat)
!!@endverbatim
!!
!!@n @param nri Number of radial points
!!@n @param jmax Number of spherical harmonics modes
-!!@n @param kr_st Start radial address to construct matrix
-!!@n @param kr_ed End radial address to construct matrix
+!!@n @param kr_in Radial ID for inner boundary
+!!@n @param kr_out Radial ID for outer boundary
+!!@n @param r_ICB(0:2) Radius at ICB
+!!@n @param r_CMB(0:2) Radius at CMB
!!@n @param coef_imp Coefficient for contribution of implicit term
!!@n @param coef_d Coefficient of diffusiotn term
+!!@n @param fdm2_fix_dr_ICB(-1:1,3)
+!! Matrix to evaluate field at ICB with fiexed radial derivative
+!!@n @param fdm2_fix_dr_CMB(-1:1,3)
+!! Matrix to evaluate field at CMB with fiexed radial derivative
!!
!!@n @param evo_mat(3,nri,jmax) Band matrix for time evolution
!
@@ -43,30 +51,31 @@
!
! -----------------------------------------------------------------------
!
- subroutine set_fix_scalar_icb_rmat_sph(nri, jmax, evo_mat)
+ subroutine set_fix_scalar_icb_rmat_sph(nri, jmax, kr_in, &
+ & evo_mat)
!
- integer(kind = kint), intent(in) :: jmax, nri
+ integer(kind = kint), intent(in) :: jmax, nri, kr_in
real(kind = kreal), intent(inout) :: evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
do j = 1, jmax
- evo_mat(2,nlayer_ICB, j) = one
- evo_mat(1,nlayer_ICB+1,j) = zero
+ evo_mat(2,kr_in, j) = one
+ evo_mat(1,kr_in+1,j) = zero
end do
!
end subroutine set_fix_scalar_icb_rmat_sph
!
! -----------------------------------------------------------------------
!
- subroutine set_fix_flux_icb_rmat_sph(nri, jmax, &
- & coef_imp, coef_d, evo_mat)
+ subroutine set_fix_flux_icb_rmat_sph(nri, jmax, kr_in, r_ICB, &
+ & fdm2_fix_dr_ICB, coef_imp, coef_d, evo_mat)
!
- use m_coef_fdm_fixed_ICB
-!
- integer(kind = kint), intent(in) :: jmax, nri
+ integer(kind = kint), intent(in) :: jmax, nri, kr_in
real(kind = kreal), intent(in) :: coef_imp, coef_d
+ real(kind = kreal), intent(in) :: r_ICB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_ICB(-1:1,3)
!
real(kind = kreal), intent(inout) :: evo_mat(3,nri,jmax)
!
@@ -74,11 +83,11 @@
!
!
do j = 1, jmax
- evo_mat(2,nlayer_ICB, j) = one + coef_imp*dt*coef_d &
- & * ( -coef_fdm_fix_dr_ICB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_ICB,2))
- evo_mat(1,nlayer_ICB+1,j) = - coef_imp*dt*coef_d &
- & * coef_fdm_fix_dr_ICB_2( 1,3)
+ evo_mat(2,kr_in, j) = one + coef_imp*dt*coef_d &
+ & * ( -fdm2_fix_dr_ICB( 0,3) &
+ & + g_sph_rj(j,3)*r_ICB(2))
+ evo_mat(1,kr_in+1,j) = - coef_imp*dt*coef_d &
+ & * fdm2_fix_dr_ICB( 1,3)
end do
!
end subroutine set_fix_flux_icb_rmat_sph
@@ -86,30 +95,31 @@
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
!
- subroutine set_fix_scalar_cmb_rmat_sph(nri, jmax, evo_mat)
+ subroutine set_fix_scalar_cmb_rmat_sph(nri, jmax, kr_out, &
+ & evo_mat)
!
- integer(kind = kint), intent(in) :: jmax, nri
+ integer(kind = kint), intent(in) :: jmax, nri, kr_out
real(kind = kreal), intent(inout) :: evo_mat(3,nri,jmax)
!
integer(kind = kint) :: j
!
!
do j = 1, jmax
- evo_mat(3,nlayer_CMB-1,j) = zero
- evo_mat(2,nlayer_CMB, j) = one
+ evo_mat(3,kr_out-1,j) = zero
+ evo_mat(2,kr_out, j) = one
end do
!
end subroutine set_fix_scalar_cmb_rmat_sph
!
! -----------------------------------------------------------------------
!
- subroutine set_fix_flux_cmb_rmat_sph(nri, jmax, &
- & coef_imp, coef_d, evo_mat)
-!
- use m_coef_fdm_fixed_CMB
+ subroutine set_fix_flux_cmb_rmat_sph(nri, jmax, kr_out, r_CMB, &
+ & fdm2_fix_dr_CMB, coef_imp, coef_d, evo_mat)
!
- integer(kind = kint), intent(in) :: jmax, nri
+ integer(kind = kint), intent(in) :: jmax, nri, kr_out
real(kind = kreal), intent(in) :: coef_imp, coef_d
+ real(kind = kreal), intent(in) :: r_CMB(0:2)
+ real(kind = kreal), intent(in) :: fdm2_fix_dr_CMB(-1:1,3)
!
real(kind = kreal), intent(inout) :: evo_mat(3,nri,jmax)
!
@@ -117,11 +127,11 @@
!
!
do j = 1, jmax
- evo_mat(3,nlayer_CMB-1,j) = - coef_imp*dt*coef_d &
- & * coef_fdm_fix_dr_CMB_2(-1,3)
- evo_mat(2,nlayer_CMB, j) = one + coef_imp*dt*coef_d &
- & * ( -coef_fdm_fix_dr_CMB_2( 0,3) &
- & + g_sph_rj(j,3)*ar_1d_rj(nlayer_CMB,2))
+ evo_mat(3,kr_out-1,j) = - coef_imp*dt*coef_d &
+ & * fdm2_fix_dr_CMB(-1,3)
+ evo_mat(2,kr_out, j) = one + coef_imp*dt*coef_d &
+ & * ( -fdm2_fix_dr_CMB( 0,3) &
+ & + g_sph_rj(j,3)*r_CMB(2))
end do
!
end subroutine set_fix_flux_cmb_rmat_sph
More information about the CIG-COMMITS
mailing list