diff --git a/src/Fortran_libraries/MHD_src/radial_FDM/Makefile.depends b/src/Fortran_libraries/MHD_src/radial_FDM/Makefile.depends index 56e1238f..b6142c4b 100644 --- a/src/Fortran_libraries/MHD_src/radial_FDM/Makefile.depends +++ b/src/Fortran_libraries/MHD_src/radial_FDM/Makefile.depends @@ -1,3 +1,9 @@ +center_sph_matrices.o: $(MHD_FDM_DIR)/center_sph_matrices.f90 m_precision.o m_constants.o + $(F90) -c $(F90OPTFLAGS) $< +select_sph_r_mat_magne_BC.o: $(MHD_FDM_DIR)/select_sph_r_mat_magne_BC.f90 m_precision.o m_constants.o m_machine_parameter.o m_ludcmp_3band.o t_physical_property.o t_spheric_rj_data.o t_sph_matrices.o t_boundary_params_sph_MHD.o t_coef_fdm2_centre.o set_radial_mat_sph.o set_sph_scalar_matrix_ICB.o center_sph_matrices.o set_sph_scalar_matrix_CMB.o + $(F90) -c $(F90OPTFLAGS) $< +select_sph_r_mat_vort_BC.o: $(MHD_FDM_DIR)/select_sph_r_mat_vort_BC.f90 m_precision.o calypso_mpi.o m_constants.o m_machine_parameter.o m_ludcmp_3band.o t_physical_property.o t_spheric_rj_data.o t_sph_matrices.o t_boundary_params_sph_MHD.o t_coef_sph_velocity_BCs.o t_coef_fdm2_centre.o set_radial_mat_sph.o set_sph_scalar_matrix_ICB.o center_sph_matrices.o set_sph_scalar_matrix_CMB.o + $(F90) -c $(F90OPTFLAGS) $< set_sph_scalar_matrix_CMB.o: $(MHD_FDM_DIR)/set_sph_scalar_matrix_CMB.f90 m_precision.o m_constants.o $(F90) -c $(F90OPTFLAGS) $< set_sph_scalar_matrix_ICB.o: $(MHD_FDM_DIR)/set_sph_scalar_matrix_ICB.f90 m_precision.o m_constants.o diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/center_sph_matrices.f90 b/src/Fortran_libraries/MHD_src/radial_FDM/center_sph_matrices.f90 similarity index 100% rename from src/Fortran_libraries/MHD_src/sph_MHD/center_sph_matrices.f90 rename to src/Fortran_libraries/MHD_src/radial_FDM/center_sph_matrices.f90 diff --git a/src/Fortran_libraries/MHD_src/radial_FDM/select_sph_r_mat_magne_BC.f90 b/src/Fortran_libraries/MHD_src/radial_FDM/select_sph_r_mat_magne_BC.f90 new file mode 100644 index 00000000..47715680 --- /dev/null +++ b/src/Fortran_libraries/MHD_src/radial_FDM/select_sph_r_mat_magne_BC.f90 @@ -0,0 +1,165 @@ +!>@file select_sph_r_mat_magne_BC.f90 +!!@brief module select_sph_r_mat_magne_BC +!! +!!@date Programmed by H.Matsui on Apr., 2009 +! +!>@brief Construct matrix for magnetic field at boundaries +!! +!!@verbatim +!! subroutine sel_sph_r_mat_pol_magnetic_ICB(sph_rj, sph_bc_B, & +!! & fdm2_center, g_sph_rj, coef_dbt, band_bp_evo) +!! subroutine sel_sph_r_mat_tor_magnetic_ICB(sph_rj, sph_bc_B, & +!! & fdm2_center, g_sph_rj, coef_dbt, band_bt_evo) +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(sph_boundary_type), intent(in) :: sph_bc_B +!! type(fdm2_center_mat), intent(in) :: fdm2_center +!! real(kind=kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) +!! real(kind = kreal), intent(in) :: coef_dbt +!! type(band_matrices_type), intent(inout) :: band_bt_evo +!! +!! subroutine sel_sph_r_mat_pol_magnetic_CMB(sph_rj, sph_bc_B, & +!! & g_sph_rj, coef_dbt, band_bp_evo) +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(sph_boundary_type), intent(in) :: sph_bc_B +!! real(kind = kreal), intent(in) :: coef_dbt +!! real(kind=kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) +!! type(band_matrices_type), intent(inout) :: band_bp_evo +!!@endverbatim +! + module select_sph_r_mat_magne_BC +! + use m_precision +! + use m_constants + use m_machine_parameter + use m_ludcmp_3band +! + use t_physical_property + use t_spheric_rj_data + use t_sph_matrices + use t_boundary_params_sph_MHD + use t_coef_fdm2_centre +! + use set_radial_mat_sph +! + implicit none +! +! ----------------------------------------------------------------------- +! + contains +! +! ----------------------------------------------------------------------- +! + subroutine sel_sph_r_mat_pol_magnetic_ICB(sph_rj, sph_bc_B, & + & fdm2_center, g_sph_rj, coef_dbt, band_bp_evo) +! + use set_sph_scalar_matrix_ICB + use center_sph_matrices +! + type(sph_rj_grid), intent(in) :: sph_rj + type(sph_boundary_type), intent(in) :: sph_bc_B + type(fdm2_center_mat), intent(in) :: fdm2_center +! + real(kind = kreal), intent(in) :: coef_dbt + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) +! + type(band_matrices_type), intent(inout) :: band_bp_evo +! +! + if(sph_bc_B%iflag_icb .eq. iflag_sph_fill_center) then + call add_vector_poisson_mat_center & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_B%r_ICB, fdm2_center%dmat_fix_fld, & + & coef_dbt, band_bp_evo%mat) + else if(sph_bc_B%iflag_icb .eq. iflag_radial_magne) then + call add_fix_flux_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_B%kr_in, sph_bc_B%r_ICB, sph_bc_B%fdm2_fix_dr_ICB, & + & coef_dbt, band_bp_evo%mat) + else if(sph_bc_B%iflag_icb .eq. iflag_evolve_field) then + call set_fix_fld_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_B%kr_in, band_bp_evo%mat) +! else if(sph_bc_B%iflag_icb .eq. iflag_sph_insulator) then + else + call set_ins_magne_icb_rmat_sph & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_B%kr_in, sph_bc_B%r_ICB, sph_bc_B%fdm2_fix_dr_ICB, & + & coef_dbt, band_bp_evo%mat) + end if +! + end subroutine sel_sph_r_mat_pol_magnetic_ICB +! +! ----------------------------------------------------------------------- +! + subroutine sel_sph_r_mat_tor_magnetic_ICB(sph_rj, sph_bc_B, & + & fdm2_center, g_sph_rj, coef_dbt, band_bt_evo) +! + use set_sph_scalar_matrix_ICB + use center_sph_matrices +! + real(kind = kreal), intent(in) :: coef_dbt + type(sph_rj_grid), intent(in) :: sph_rj + type(sph_boundary_type), intent(in) :: sph_bc_B + type(fdm2_center_mat), intent(in) :: fdm2_center +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) +! + type(band_matrices_type), intent(inout) :: band_bt_evo +! +! + if(sph_bc_B%iflag_icb .eq. iflag_sph_fill_center) then + call add_vector_poisson_mat_center & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_B%r_ICB, fdm2_center%dmat_fix_fld, & + & coef_dbt, band_bt_evo%mat) +! else if(sph_bc_B%iflag_icb .eq. iflag_radial_magne) then +! else if(sph_bc_B%iflag_icb .eq. iflag_evolve_field) then +! else if(sph_bc_B%iflag_icb .eq. iflag_sph_insulator) then + else + call set_fix_fld_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_B%kr_in, band_bt_evo%mat) + end if +! + end subroutine sel_sph_r_mat_tor_magnetic_ICB +! +! ----------------------------------------------------------------------- +! ----------------------------------------------------------------------- +! + subroutine sel_sph_r_mat_pol_magnetic_CMB(sph_rj, sph_bc_B, & + & g_sph_rj, coef_dbt, band_bp_evo) +! + use set_sph_scalar_matrix_CMB +! + type(sph_rj_grid), intent(in) :: sph_rj + type(sph_boundary_type), intent(in) :: sph_bc_B + real(kind = kreal), intent(in) :: coef_dbt + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) +! + type(band_matrices_type), intent(inout) :: band_bp_evo +! +! +! + if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then + call add_fix_flux_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_B%kr_out, sph_bc_B%r_CMB, sph_bc_B%fdm2_fix_dr_CMB, & + & coef_dbt, band_bp_evo%mat) + else if(sph_bc_B%iflag_cmb .eq. iflag_evolve_field) then + call set_fix_fld_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_B%kr_out, band_bp_evo%mat) +! else if(sph_bc_B%iflag_cmb .eq. iflag_sph_insulator) then + else + call set_ins_magne_cmb_rmat_sph & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_B%kr_out, sph_bc_B%r_CMB, sph_bc_B%fdm2_fix_dr_CMB, & + & coef_dbt, band_bp_evo%mat) + end if +! + end subroutine sel_sph_r_mat_pol_magnetic_CMB +! +! ----------------------------------------------------------------------- +! + end module select_sph_r_mat_magne_BC diff --git a/src/Fortran_libraries/MHD_src/radial_FDM/select_sph_r_mat_vort_BC.f90 b/src/Fortran_libraries/MHD_src/radial_FDM/select_sph_r_mat_vort_BC.f90 new file mode 100644 index 00000000..5928b0c7 --- /dev/null +++ b/src/Fortran_libraries/MHD_src/radial_FDM/select_sph_r_mat_vort_BC.f90 @@ -0,0 +1,244 @@ +!>@file select_sph_r_mat_vort_BC.f90 +!!@brief module select_sph_r_mat_vort_BC +!! +!!@date Programmed by H.Matsui on Apr., 2009 +! +!>@brief Construct matrix for time evolution of vector fields +!! +!!@verbatim +!! subroutine sel_sph_r_mat_vort_2step_ICB(sph_rj, sph_bc_U, & +!! & bc_fdms_U, fdm2_center, g_sph_rj, coef_dvt, & +!! & band_vs_poisson, band_wt_evo) +!! subroutine sel_sph_r_mat_tor_flow_ICB(sph_rj, sph_bc_U, & +!! & bc_fdms_U, fdm2_center, g_sph_rj, coef_dvt, & +!! & band_vt_evo) +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(sph_boundary_type), intent(in) :: sph_bc_U +!! type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U +!! type(fdm2_center_mat), intent(in) :: fdm2_center +!! real(kind=kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) +!! real(kind = kreal), intent(in) :: coef_dvt +!! type(band_matrices_type), intent(inout) :: band_wt_evo +!! type(band_matrices_type), intent(inout) :: band_vs_poisson +!! type(band_matrices_type), intent(inout) :: band_vt_evo +!! +!! subroutine sel_sph_r_mat_vort_2step_CMB & +!! & (sph_rj, sph_bc_U, bc_fdms_U, g_sph_rj, coef_dvt, & +!! & band_vs_poisson, band_wt_evo) +!! subroutine sel_sph_r_mat_tor_flow_CMB(sph_rj, sph_bc_U, & +!! & bc_fdms_U, g_sph_rj, coef_dvt, band_vt_evo) +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(sph_boundary_type), intent(in) :: sph_bc_U +!! type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U +!! real(kind=kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) +!! real(kind = kreal), intent(in) :: coef_dvt +!! type(band_matrices_type), intent(inout) :: band_vs_poisson +!! type(band_matrices_type), intent(inout) :: band_wt_evo +!! type(band_matrices_type), intent(inout) :: band_vt_evo +!!@endverbatim +! + module select_sph_r_mat_vort_BC +! + use m_precision + use calypso_mpi +! + use m_constants + use m_machine_parameter + use m_ludcmp_3band +! + use t_physical_property + use t_spheric_rj_data + use t_sph_matrices + use t_boundary_params_sph_MHD + use t_coef_sph_velocity_BCs + use t_coef_fdm2_centre +! + use set_radial_mat_sph +! + implicit none +! +! ----------------------------------------------------------------------- +! + contains +! +! ----------------------------------------------------------------------- +! + subroutine sel_sph_r_mat_vort_2step_ICB(sph_rj, sph_bc_U, & + & bc_fdms_U, fdm2_center, g_sph_rj, coef_dvt, & + & band_vs_poisson, band_wt_evo) +! + use set_sph_scalar_matrix_ICB + use center_sph_matrices +! + type(sph_rj_grid), intent(in) :: sph_rj + type(sph_boundary_type), intent(in) :: sph_bc_U + type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U + type(fdm2_center_mat), intent(in) :: fdm2_center +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) + real(kind = kreal), intent(in) :: coef_dvt +! + type(band_matrices_type), intent(inout) :: band_wt_evo + type(band_matrices_type), intent(inout) :: band_vs_poisson +! +! + if(sph_bc_U%iflag_icb .eq. iflag_sph_fill_center) then + call add_vector_poisson_mat_center & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%r_ICB, fdm2_center%dmat_fix_fld, & + & coef_dvt, band_wt_evo%mat) + call add_vector_poisson_mat_center & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%r_ICB, fdm2_center%dmat_fix_fld, & + & one, band_vs_poisson%mat) + else + call add_fix_flux_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_in, sph_bc_U%r_ICB, sph_bc_U%fdm2_fix_dr_ICB, & + & coef_dvt, band_wt_evo%mat) +! + if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then + call add_fix_flux_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_in, sph_bc_U%r_ICB, & + & bc_fdms_U%fdm2_free_ICB%dmat_vp, & + & one, band_vs_poisson%mat) +! else if(sph_bc_U%iflag_icb .eq. iflag_evolve_field) then +! else if(sph_bc_U%iflag_icb .eq. iflag_fixed_field) then +! else if(sph_bc_U%iflag_icb .eq. iflag_fixed_velo) then + else + call add_fix_flux_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_in, sph_bc_U%r_ICB, sph_bc_U%fdm2_fix_dr_ICB, & + & one, band_vs_poisson%mat) + end if + end if +! + end subroutine sel_sph_r_mat_vort_2step_ICB +! +! ----------------------------------------------------------------------- +! + subroutine sel_sph_r_mat_tor_flow_ICB(sph_rj, sph_bc_U, & + & bc_fdms_U, fdm2_center, g_sph_rj, coef_dvt, & + & band_vt_evo) +! + use set_sph_scalar_matrix_ICB + use center_sph_matrices +! + type(sph_rj_grid), intent(in) :: sph_rj + type(sph_boundary_type), intent(in) :: sph_bc_U + type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U + type(fdm2_center_mat), intent(in) :: fdm2_center +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) + real(kind = kreal), intent(in) :: coef_dvt +! + type(band_matrices_type), intent(inout) :: band_vt_evo +! +! + if(sph_bc_U%iflag_icb .eq. iflag_sph_fill_center) then + call add_vector_poisson_mat_center & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%r_ICB, fdm2_center%dmat_fix_fld, & + & coef_dvt, band_vt_evo%mat) + else + if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then + call add_fix_flux_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_in, sph_bc_U%r_ICB, & + & bc_fdms_U%fdm2_free_ICB%dmat_vt, & + & coef_dvt, band_vt_evo%mat) +! else if(sph_bc_U%iflag_icb .eq. iflag_evolve_field) then +! else if(sph_bc_U%iflag_icb .eq. iflag_fixed_field) then +! else if(sph_bc_U%iflag_icb .eq. iflag_fixed_velo) then +! else if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then + else + call set_fix_fld_icb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_U%kr_in, band_vt_evo%mat) + end if + end if +! + end subroutine sel_sph_r_mat_tor_flow_ICB +! +! ----------------------------------------------------------------------- +! ----------------------------------------------------------------------- +! + subroutine sel_sph_r_mat_vort_2step_CMB & + & (sph_rj, sph_bc_U, bc_fdms_U, g_sph_rj, coef_dvt, & + & band_vs_poisson, band_wt_evo) +! + use set_sph_scalar_matrix_CMB + use center_sph_matrices +! + type(sph_rj_grid), intent(in) :: sph_rj + type(sph_boundary_type), intent(in) :: sph_bc_U + type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) + real(kind = kreal), intent(in) :: coef_dvt +! + type(band_matrices_type), intent(inout) :: band_vs_poisson + type(band_matrices_type), intent(inout) :: band_wt_evo +! +! + if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then + call add_fix_flux_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_out, sph_bc_U%r_CMB, & + & bc_fdms_U%fdm2_free_CMB%dmat_vp, & + & one, band_vs_poisson%mat) +! else if(sph_bc_U%iflag_cmb .eq. iflag_evolve_field) then +! else if(sph_bc_U%iflag_cmb .eq. iflag_fixed_field) then +! else if(sph_bc_U%iflag_cmb .eq. iflag_fixed_velo) then + else + call add_fix_flux_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_out, sph_bc_U%r_CMB, sph_bc_U%fdm2_fix_dr_CMB, & + & one, band_vs_poisson%mat) + end if +! + call add_fix_flux_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_out, sph_bc_U%r_CMB, sph_bc_U%fdm2_fix_dr_CMB, & + & coef_dvt, band_wt_evo%mat) +! + end subroutine sel_sph_r_mat_vort_2step_CMB +! +! ----------------------------------------------------------------------- +! + subroutine sel_sph_r_mat_tor_flow_CMB(sph_rj, sph_bc_U, & + & bc_fdms_U, g_sph_rj, coef_dvt, band_vt_evo) +! + use set_sph_scalar_matrix_CMB +! + type(sph_rj_grid), intent(in) :: sph_rj + type(sph_boundary_type), intent(in) :: sph_bc_U + type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) + real(kind = kreal), intent(in) :: coef_dvt +! + type(band_matrices_type), intent(inout) :: band_vt_evo +! +! + if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then + call add_fix_flux_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & + & sph_bc_U%kr_out, sph_bc_U%r_CMB, & + & bc_fdms_U%fdm2_free_CMB%dmat_vt, & + & coef_dvt, band_vt_evo%mat) +! else if(sph_bc_U%iflag_cmb .eq. iflag_evolve_field) then +! else if(sph_bc_U%iflag_cmb .eq. iflag_fixed_field) then +! else if(sph_bc_U%iflag_cmb .eq. iflag_fixed_velo) then + else + call set_fix_fld_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_U%kr_out, band_vt_evo%mat) + end if +! + end subroutine sel_sph_r_mat_tor_flow_CMB +! +! ----------------------------------------------------------------------- +! + end module select_sph_r_mat_vort_BC diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/Makefile.depends b/src/Fortran_libraries/MHD_src/sph_MHD/Makefile.depends index 351e7867..429e2eec 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/Makefile.depends +++ b/src/Fortran_libraries/MHD_src/sph_MHD/Makefile.depends @@ -82,23 +82,23 @@ cal_vorticity_terms_adams.o: $(MHD_SPH_DIR)/cal_vorticity_terms_adams.f90 m_prec $(F90) -c $(F90OPTFLAGS) $< cal_write_sph_monitor_data.o: $(MHD_SPH_DIR)/cal_write_sph_monitor_data.f90 m_precision.o calypso_mpi.o m_machine_parameter.o t_SPH_mesh_field_data.o t_sph_mhd_monitor_data_IO.o t_spheric_parameter.o t_phys_data.o t_boundary_data_sph_MHD.o t_boundary_sph_spectr.o t_work_4_sph_trans.o t_time_data.o t_rms_4_sph_spectr.o t_sum_sph_rms_data.o t_pickup_sph_spectr_data.o t_no_heat_Nusselt.o t_CMB_dipolarity.o t_sph_typical_scales.o t_energy_label_parameters.o t_fdm_coefs.o t_physical_property.o t_radial_matrices_sph_MHD.o t_sph_matrix.o t_solver_SR.o cal_heat_source_Nu.o field_at_mid_equator.o const_data_4_dynamobench.o cal_rms_fields_by_sph.o pickup_sph_spectr_data.o pickup_gauss_coefficients.o cal_CMB_dipolarity.o cal_typical_scale.o sph_fwd_trans_on_circles.o output_sph_pwr_volume_file.o write_picked_sph_spectr.o write_sph_gauss_coefs.o write_typical_scale.o write_dynamo_benchmark_file.o write_monitors_circle_file.o output_sph_volume_ave_file.o output_sph_pwr_layer_file.o $(F90) -c $(F90OPTFLAGS) $< -center_sph_matrices.o: $(MHD_SPH_DIR)/center_sph_matrices.f90 m_precision.o m_constants.o - $(F90) -c $(F90OPTFLAGS) $< check_dependency_for_MHD.o: $(MHD_SPH_DIR)/check_dependency_for_MHD.f90 m_precision.o m_error_IDs.o m_machine_parameter.o calypso_mpi.o t_control_parameter.o t_spheric_parameter.o t_phys_data.o t_phys_address.o t_base_field_labels.o t_physical_property.o t_SPH_mesh_field_data.o set_control_field_data.o $(F90) -c $(F90OPTFLAGS) $< const_data_4_dynamobench.o: $(MHD_SPH_DIR)/const_data_4_dynamobench.f90 m_precision.o m_constants.o m_machine_parameter.o calypso_mpi.o field_at_mid_equator.o t_spheric_parameter.o t_spheric_rj_data.o t_phys_address.o t_phys_data.o t_time_data.o t_schmidt_poly_on_rtm.o t_rms_4_sph_spectr.o t_boundary_data_sph_MHD.o t_field_on_circle.o t_field_4_dynamobench.o t_work_4_sph_trans.o calypso_mpi_real.o sph_fwd_trans_mid_eq.o cal_rms_fields_by_sph.o global_field_4_dynamobench.o transfer_to_long_integers.o $(F90) -c $(F90OPTFLAGS) $< const_diffusive_profile.o: $(MHD_SPH_DIR)/const_diffusive_profile.f90 m_precision.o m_constants.o m_machine_parameter.o t_spheric_rj_data.o t_phys_data.o t_phys_address.o t_work_4_sph_trans.o t_schmidt_poly_on_rtm.o t_control_parameter.o t_physical_property.o t_boundary_data_sph_MHD.o t_boundary_params_sph_MHD.o t_boundary_sph_spectr.o t_coef_fdm2_centre.o t_fdm_coefs.o t_sph_matrix.o t_sph_center_matrix.o calypso_mpi.o calypso_mpi_real.o const_sph_radial_grad.o fill_scalar_field.o select_exp_scalar_ICB.o select_exp_scalar_CMB.o $(F90) -c $(F90OPTFLAGS) $< +const_r_mat_4_magnetic_sph.o: $(MHD_SPH_DIR)/const_r_mat_4_magnetic_sph.f90 m_precision.o calypso_mpi.o m_constants.o m_machine_parameter.o m_ludcmp_3band.o t_physical_property.o t_spheric_rj_data.o t_sph_matrices.o t_fdm_coefs.o t_boundary_params_sph_MHD.o t_coef_fdm2_centre.o set_radial_mat_sph.o select_sph_r_mat_magne_BC.o set_sph_scalar_matrix_CMB.o center_sph_matrices.o + $(F90) -c $(F90OPTFLAGS) $< const_r_mat_4_scalar_sph.o: $(MHD_SPH_DIR)/const_r_mat_4_scalar_sph.f90 m_precision.o calypso_mpi.o m_constants.o m_machine_parameter.o t_physical_property.o t_spheric_parameter.o t_spheric_rj_data.o t_sph_matrices.o t_sph_center_matrix.o t_fdm_coefs.o t_boundary_params_sph_MHD.o t_coef_fdm2_centre.o m_ludcmp_3band.o cal_inner_core_rotation.o center_sph_matrices.o mat_product_3band_mul.o set_radial_mat_sph.o select_r_mat_scalar_bc_sph.o check_sph_radial_mat.o $(F90) -c $(F90OPTFLAGS) $< -const_r_mat_4_vector_sph.o: $(MHD_SPH_DIR)/const_r_mat_4_vector_sph.f90 m_precision.o calypso_mpi.o m_constants.o m_machine_parameter.o m_ludcmp_3band.o t_physical_property.o t_spheric_rj_data.o t_sph_matrices.o t_fdm_coefs.o t_boundary_params_sph_MHD.o set_radial_mat_sph.o t_coef_sph_velocity_BCs.o t_coef_fdm2_centre.o m_ludcmp_band.o set_sph_scalar_matrix_ICB.o set_sph_scalar_matrix_CMB.o cal_inner_core_rotation.o center_sph_matrices.o mat_product_3band_mul.o check_sph_radial_mat.o +const_r_mat_4_vector_sph.o: $(MHD_SPH_DIR)/const_r_mat_4_vector_sph.f90 m_precision.o m_constants.o m_machine_parameter.o m_ludcmp_3band.o t_physical_property.o t_spheric_rj_data.o t_sph_matrices.o t_fdm_coefs.o t_boundary_params_sph_MHD.o t_coef_sph_velocity_BCs.o t_coef_fdm2_centre.o set_radial_mat_sph.o m_ludcmp_band.o select_sph_r_mat_vort_BC.o center_sph_matrices.o mat_product_3band_mul.o cal_inner_core_rotation.o $(F90) -c $(F90OPTFLAGS) $< const_r_mat_w_center_sph.o: $(MHD_SPH_DIR)/const_r_mat_w_center_sph.f90 m_precision.o calypso_mpi.o m_constants.o m_machine_parameter.o t_physical_property.o t_spheric_rj_data.o t_sph_center_matrix.o t_boundary_params_sph_MHD.o t_coef_fdm2_centre.o t_sph_matrices.o t_sph_matrix.o m_ludcmp_3band.o center_sph_matrices.o $(F90) -c $(F90OPTFLAGS) $< const_radial_forces_on_bc.o: $(MHD_SPH_DIR)/const_radial_forces_on_bc.f90 m_precision.o m_machine_parameter.o m_constants.o t_physical_property.o t_reference_scalar_param.o t_boundary_params_sph_MHD.o t_spheric_rj_data.o t_phys_address.o t_phys_data.o self_buoyancy_on_sphere.o $(F90) -c $(F90OPTFLAGS) $< -const_radial_mat_4_sph.o: $(MHD_SPH_DIR)/const_radial_mat_4_sph.f90 m_precision.o m_constants.o m_machine_parameter.o t_control_parameter.o t_physical_property.o t_spheric_parameter.o t_spheric_rj_data.o t_fdm_coefs.o t_schmidt_poly_on_rtm.o t_boundary_data_sph_MHD.o t_boundary_params_sph_MHD.o t_radial_matrices_sph_MHD.o t_sph_matrices.o t_sph_matrix.o calypso_mpi.o const_r_mat_4_scalar_sph.o const_r_mat_w_center_sph.o const_r_mat_4_vector_sph.o +const_radial_mat_4_sph.o: $(MHD_SPH_DIR)/const_radial_mat_4_sph.f90 m_precision.o m_constants.o m_machine_parameter.o t_control_parameter.o t_physical_property.o t_spheric_parameter.o t_spheric_rj_data.o t_fdm_coefs.o t_schmidt_poly_on_rtm.o t_boundary_data_sph_MHD.o t_boundary_params_sph_MHD.o t_radial_matrices_sph_MHD.o t_sph_matrices.o t_sph_matrix.o calypso_mpi.o const_r_mat_4_scalar_sph.o const_r_mat_w_center_sph.o const_r_mat_4_vector_sph.o const_r_mat_4_magnetic_sph.o $(F90) -c $(F90OPTFLAGS) $< const_radial_references.o: $(MHD_SPH_DIR)/const_radial_references.f90 m_precision.o m_constants.o m_machine_parameter.o t_spheric_rj_data.o t_phys_data.o t_phys_address.o t_work_4_sph_trans.o t_schmidt_poly_on_rtm.o t_control_parameter.o t_physical_property.o t_boundary_data_sph_MHD.o t_boundary_params_sph_MHD.o t_boundary_sph_spectr.o t_coef_fdm2_centre.o t_fdm_coefs.o t_sph_matrix.o t_sph_center_matrix.o const_diffusive_profile.o $(F90) -c $(F90OPTFLAGS) $< @@ -202,7 +202,7 @@ set_address_sph_trans_snap.o: $(MHD_SPH_DIR)/set_address_sph_trans_snap.f90 m_pr $(F90) -c $(F90OPTFLAGS) $< set_bc_flag_sph_velo.o: $(MHD_SPH_DIR)/set_bc_flag_sph_velo.f90 m_precision.o m_constants.o m_machine_parameter.o m_boundary_condition_IDs.o t_boundary_data_sph_MHD.o t_boundary_params_sph_MHD.o t_boundary_sph_spectr.o t_sph_boundary_input_data.o t_spheric_rj_data.o t_group_data.o t_bc_data_list.o set_bc_sph_scalars.o m_base_field_labels.o set_sph_bc_data_by_file.o $(F90) -c $(F90OPTFLAGS) $< -set_bc_sph_mhd.o: $(MHD_SPH_DIR)/set_bc_sph_mhd.f90 m_precision.o m_machine_parameter.o m_boundary_condition_IDs.o t_control_parameter.o t_physical_property.o t_spheric_parameter.o t_group_data.o t_boundary_data_sph_MHD.o t_boundary_params_sph_MHD.o t_spheric_rj_data.o t_bc_data_list.o t_sph_boundary_input_data.o t_coef_sph_velocity_BCs.o t_coef_fdm2_centre.o t_phys_data.o m_base_field_labels.o set_bc_flag_sph_velo.o set_bc_sph_scalars.o set_sph_bc_magne_sph.o t_coef_fdm1_free_rotate_ICB.o t_coef_fdm1_free_rotate_CMB.o +set_bc_sph_mhd.o: $(MHD_SPH_DIR)/set_bc_sph_mhd.f90 m_precision.o m_machine_parameter.o m_boundary_condition_IDs.o t_control_parameter.o t_physical_property.o t_spheric_parameter.o t_group_data.o t_phys_data.o t_boundary_data_sph_MHD.o t_boundary_params_sph_MHD.o t_spheric_rj_data.o t_bc_data_list.o t_sph_boundary_input_data.o t_coef_sph_velocity_BCs.o t_coef_fdm2_centre.o m_base_field_labels.o set_bc_flag_sph_velo.o set_bc_sph_scalars.o set_sph_bc_magne_sph.o t_coef_fdm1_free_rotate_ICB.o t_coef_fdm1_free_rotate_CMB.o $(F90) -c $(F90OPTFLAGS) $< set_bc_sph_scalars.o: $(MHD_SPH_DIR)/set_bc_sph_scalars.f90 m_precision.o calypso_mpi.o m_constants.o m_error_IDs.o m_machine_parameter.o m_boundary_condition_IDs.o t_spheric_rj_data.o t_group_data.o t_boundary_params_sph_MHD.o t_boundary_sph_spectr.o t_sph_boundary_input_data.o t_bc_data_list.o t_field_labels.o m_base_field_labels.o m_base_force_labels.o set_sph_bc_data_by_file.o $(F90) -c $(F90OPTFLAGS) $< diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_heat_source_Nu.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_heat_source_Nu.f90 index 2c956584..8794adbf 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_heat_source_Nu.f90 +++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_heat_source_Nu.f90 @@ -83,7 +83,7 @@ subroutine init_poisson_matrix_for_Nu & if(Nu_type%iflag_Nusselt .eq. iflag_no_source_Nu) return call alloc_Nu_radial_reference(sph%sph_rj, Nu_type) call const_r_mat00_scalar_sph & - & (mat_name, sc_prop%diffusie_reduction_ICB, & + & (my_rank+50, mat_name, sc_prop%diffusie_reduction_ICB, & & sph%sph_params, sph%sph_rj, r_2nd, sph_bc_S, fdm2_center, & & band_s00_poisson_fixS) ! diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/cal_inner_core_rotation.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/cal_inner_core_rotation.f90 index 162da46b..5615c656 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/cal_inner_core_rotation.f90 +++ b/src/Fortran_libraries/MHD_src/sph_MHD/cal_inner_core_rotation.f90 @@ -74,7 +74,7 @@ module cal_inner_core_rotation implicit none ! private :: int_icore_tor_lorentz_l1, cal_icore_viscous_drag_l1 - private :: set_rotate_icb_vt_sph_mat, set_inner_core_rot_l1 + private :: set_inner_core_rot_l1 ! ! ---------------------------------------------------------------------- ! diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_magnetic_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_magnetic_sph.f90 new file mode 100644 index 00000000..c7e1f02a --- /dev/null +++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_magnetic_sph.f90 @@ -0,0 +1,134 @@ +!>@file const_r_mat_4_magnetic_sph.f90 +!!@brief module const_r_mat_4_magnetic_sph +!! +!!@date Programmed by H.Matsui on Apr., 2009 +! +!>@brief Construct matrix for time evolution of vector fields +!! +!!@verbatim +!! subroutine const_radial_mat_4_magne_sph & +!! & (dt, sph_rj, r_2nd, cd_prop, sph_bc_B, fdm2_center, & +!! & g_sph_rj, band_bp_evo, band_bt_evo) +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(fdm_matrices), intent(in) :: r_2nd +!! type(conductive_property), intent(in) :: cd_prop +!! type(sph_boundary_type), intent(in) :: sph_bc_U +!! type(sph_boundary_type), intent(in) :: sph_bc_B +!! type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U +!! type(fdm2_center_mat), intent(in) :: fdm2_center +!! type(band_matrices_type), intent(inout) :: band_bp_evo +!! type(band_matrices_type), intent(inout) :: band_bt_evo +!!@endverbatim +! + module const_r_mat_4_magnetic_sph +! + use m_precision + use calypso_mpi +! + use m_constants + use m_machine_parameter + use m_ludcmp_3band +! + use t_physical_property + use t_spheric_rj_data + use t_sph_matrices + use t_fdm_coefs + use t_boundary_params_sph_MHD + use t_coef_fdm2_centre +! + use set_radial_mat_sph +! + implicit none +! + character(len=kchara), parameter, private & + & :: bp_evo_name = 'poloidal_magne_evolution' + character(len=kchara), parameter, private & + & :: bt_evo_name = 'toroidal_magne_evolution' +! +! ----------------------------------------------------------------------- +! + contains +! +! ----------------------------------------------------------------------- +! + subroutine const_radial_mat_4_magne_sph & + & (dt, sph_rj, r_2nd, cd_prop, sph_bc_B, fdm2_center, & + & g_sph_rj, band_bp_evo, band_bt_evo) +! + use select_sph_r_mat_magne_BC + use set_sph_scalar_matrix_CMB + use center_sph_matrices +! + type(sph_rj_grid), intent(in) :: sph_rj + type(fdm_matrices), intent(in) :: r_2nd + type(conductive_property), intent(in) :: cd_prop + type(sph_boundary_type), intent(in) :: sph_bc_B + type(fdm2_center_mat), intent(in) :: fdm2_center +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) + real(kind = kreal), intent(in) :: dt +! + type(band_matrices_type), intent(inout) :: band_bp_evo + type(band_matrices_type), intent(inout) :: band_bt_evo +! + real(kind = kreal) :: coef_dbt +! +! + band_bp_evo%mat_name = bp_evo_name + band_bt_evo%mat_name = bt_evo_name +! + call alloc_band_mat_sph(ithree, sph_rj, band_bp_evo) + call alloc_band_mat_sph(ithree, sph_rj, band_bt_evo) +! + call set_unit_on_diag(band_bp_evo) + call set_unit_on_diag(band_bt_evo) +! + if(cd_prop%coef_diffuse .eq. zero) then + coef_dbt = one + call set_unit_mat_4_poisson & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_B%kr_in, sph_bc_B%kr_out, band_bp_evo%mat) + call set_unit_mat_4_poisson & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_B%kr_in, sph_bc_B%kr_out, band_bt_evo%mat) + else + coef_dbt = cd_prop%coef_imp * cd_prop%coef_diffuse * dt + call set_unit_mat_4_time_evo & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_bp_evo%mat) + call set_unit_mat_4_time_evo & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_bt_evo%mat) + end if +! + call add_vector_poisson_mat_sph & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, & + & g_sph_rj, sph_bc_B%kr_in, sph_bc_B%kr_out, & + & coef_dbt, r_2nd%fdm(2)%dmat, band_bp_evo%mat) + call add_vector_poisson_mat_sph & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, & + & g_sph_rj, sph_bc_B%kr_in, sph_bc_B%kr_out, & + & coef_dbt, r_2nd%fdm(2)%dmat, band_bt_evo%mat) +! +! Matrices at ICB or center + call sel_sph_r_mat_pol_magnetic_ICB(sph_rj, sph_bc_B, & + & fdm2_center, g_sph_rj, coef_dbt, band_bp_evo) + call sel_sph_r_mat_tor_magnetic_ICB(sph_rj, sph_bc_B, & + & fdm2_center, g_sph_rj, coef_dbt, band_bt_evo) +! +! Matrices at CMB + call sel_sph_r_mat_pol_magnetic_CMB(sph_rj, sph_bc_B, g_sph_rj, & + & coef_dbt, band_bp_evo) + call set_fix_fld_cmb_poisson_mat & + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & + & sph_bc_B%kr_out, band_bt_evo%mat) +! +! + call ludcmp_3band_mul_t & + & (np_smp, sph_rj%istack_rj_j_smp, band_bp_evo) + call ludcmp_3band_mul_t & + & (np_smp, sph_rj%istack_rj_j_smp, band_bt_evo) +! + end subroutine const_radial_mat_4_magne_sph +! +! ----------------------------------------------------------------------- +! + end module const_r_mat_4_magnetic_sph 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 c460e727..e6ab6744 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 @@ -19,12 +19,13 @@ !! type(sph_boundary_type), intent(in) :: sph_bc_U !! type(fdm2_center_mat), intent(in) :: fdm2_center !! -!! subroutine const_r_mat00_scalar_sph & -!! & (mat_name, diffusie_reduction_ICB, sph_params, sph_rj, & -!! & r_2nd, sph_bc, fdm2_center, band_s00_poisson) -!! subroutine const_r_mat00_poisson_fixS & -!! & (mat_name, diffusie_reduction_ICB, sph_params, sph_rj, & -!! & r_2nd, sph_bc, fdm2_center, band_s00_poisson) +!! subroutine const_r_mat00_scalar_sph(id_file, mat_name, & +!! & diffusie_reduction_ICB, sph_params, sph_rj, r_2nd, & +!! & sph_bc, fdm2_center, band_s00_poisson) +!! subroutine const_r_mat00_poisson_fixS(id_file, mat_name, & +!! & diffusie_reduction_ICB, sph_params, sph_rj, r_2nd, & +!! & sph_bc, fdm2_center, band_s00_poisson) +!! integer(kind = kint), intent(in) :: id_file !! type(sph_shell_parameters), intent(in) :: sph_params !! type(sph_rj_grid), intent(in) :: sph_rj !! type(fdm_matrices), intent(in) :: r_2nd @@ -70,7 +71,6 @@ subroutine const_radial_mat_4_press_sph & use mat_product_3band_mul use set_radial_mat_sph use select_r_mat_scalar_bc_sph - use check_sph_radial_mat ! type(fluid_property), intent(in) :: fl_prop type(sph_boundary_type), intent(in) :: sph_bc_U @@ -86,6 +86,7 @@ subroutine const_radial_mat_4_press_sph & real(kind = kreal), allocatable :: r_coef(:) ! ! + write(band_p_poisson%mat_name,'(a)') 'pressure_poisson' coef_p = - fl_prop%coef_press ! allocate(r_coef(sph_rj%nidx_rj(1))) @@ -109,11 +110,6 @@ subroutine const_radial_mat_4_press_sph & ! call ludcmp_3band_mul_t & & (np_smp, sph_rj%istack_rj_j_smp, band_p_poisson) -! - if(i_debug .eq. iflag_full_msg) then - write(band_p_poisson%mat_name,'(a)') 'pressure_poisson' - call check_radial_band_mat(my_rank, sph_rj, band_p_poisson) - end if ! end subroutine const_radial_mat_4_press_sph ! @@ -128,7 +124,6 @@ subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect, & use center_sph_matrices use set_radial_mat_sph use select_r_mat_scalar_bc_sph - use check_sph_radial_mat ! type(sph_shell_parameters), intent(in) :: sph_params type(sph_rj_grid), intent(in) :: sph_rj @@ -188,24 +183,21 @@ subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect, & ! call ludcmp_3band_mul_t & & (np_smp, sph_rj%istack_rj_j_smp, band_s_evo) -! - if(i_debug .eq. iflag_full_msg) then - call check_radial_band_mat(my_rank, sph_rj, band_s_evo) - end if ! end subroutine const_radial_mat_4_scalar_sph ! ! ----------------------------------------------------------------------- ! - subroutine const_r_mat00_scalar_sph & - & (mat_name, diffusie_reduction_ICB, sph_params, sph_rj, & - & r_2nd, sph_bc, fdm2_center, band_s00_poisson) + subroutine const_r_mat00_scalar_sph(id_file, mat_name, & + & diffusie_reduction_ICB, sph_params, sph_rj, r_2nd, & + & sph_bc, fdm2_center, band_s00_poisson) ! use m_ludcmp_3band use set_radial_mat_sph use select_r_mat_scalar_bc_sph use check_sph_radial_mat ! + integer(kind = kint), intent(in) :: id_file type(sph_shell_parameters), intent(in) :: sph_params type(sph_rj_grid), intent(in) :: sph_rj type(fdm_matrices), intent(in) :: r_2nd @@ -264,22 +256,24 @@ subroutine const_r_mat00_scalar_sph & call ludcmp_3band_ctr(band_s00_poisson) ! if(i_debug .eq. iflag_full_msg) then - call check_center_band_matrix(my_rank, sph_rj, band_s00_poisson) + call check_center_band_matrix(id_file, sph_rj, & + & band_s00_poisson) end if ! end subroutine const_r_mat00_scalar_sph ! ! ----------------------------------------------------------------------- ! - subroutine const_r_mat00_poisson_fixS & - & (mat_name, diffusie_reduction_ICB, sph_params, sph_rj, & - & r_2nd, sph_bc, fdm2_center, band_s00_poisson) + subroutine const_r_mat00_poisson_fixS(id_file, mat_name, & + & diffusie_reduction_ICB, sph_params, sph_rj, r_2nd, & + & sph_bc, fdm2_center, band_s00_poisson) ! use m_ludcmp_3band use set_radial_mat_sph use select_r_mat_scalar_bc_sph use check_sph_radial_mat ! + integer(kind = kint), intent(in) :: id_file type(sph_shell_parameters), intent(in) :: sph_params type(sph_rj_grid), intent(in) :: sph_rj type(fdm_matrices), intent(in) :: r_2nd @@ -329,7 +323,8 @@ subroutine const_r_mat00_poisson_fixS & call ludcmp_3band_ctr(band_s00_poisson) ! if(i_debug .eq. iflag_full_msg) then - call check_center_band_matrix(my_rank, sph_rj, band_s00_poisson) + call check_center_band_matrix(id_file, sph_rj, & + & band_s00_poisson) end if ! end subroutine const_r_mat00_poisson_fixS 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 b3757ae8..152611bd 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 @@ -8,10 +8,10 @@ !!@verbatim !! subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, & !! & fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, & -!! & band_vs_poisson, band_vp_evo, band_vt_evo, band_wt_evo) -!! subroutine const_radial_mat_4_magne_sph & -!! & (dt, sph_rj, r_2nd, cd_prop, sph_bc_B, fdm2_center, & -!! & g_sph_rj, band_bp_evo, band_bt_evo) +!! & band_vs_poisson, band_vp_evo, band_wt_evo) +!! subroutine const_radial_mat_toroidal_flow(dt, sph_rj, r_2nd, & +!! & fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, & +!! & band_vt_evo) !! type(sph_rj_grid), intent(in) :: sph_rj !! type(fdm_matrices), intent(in) :: r_2nd !! type(fluid_property), intent(in) :: fl_prop @@ -24,14 +24,11 @@ !! type(band_matrices_type), intent(inout) :: band_vt_evo !! type(band_matrices_type), intent(inout) :: band_wt_evo !! type(band_matrices_type), intent(inout) :: band_vs_poisson -!! type(band_matrices_type), intent(inout) :: band_bp_evo -!! type(band_matrices_type), intent(inout) :: band_bt_evo !!@endverbatim ! module const_r_mat_4_vector_sph ! use m_precision - use calypso_mpi ! use m_constants use m_machine_parameter @@ -42,6 +39,8 @@ module const_r_mat_4_vector_sph use t_sph_matrices use t_fdm_coefs use t_boundary_params_sph_MHD + use t_coef_sph_velocity_BCs + use t_coef_fdm2_centre ! use set_radial_mat_sph ! @@ -55,11 +54,6 @@ module const_r_mat_4_vector_sph & :: vp_evo_name = 'poloidal_velocity_evolution' character(len=kchara), parameter, private & & :: vsp_evo_name = 'velocity_pressure_evolution' -! - character(len=kchara), parameter, private & - & :: bp_evo_name = 'poloidal_magne_evolution' - character(len=kchara), parameter, private & - & :: bt_evo_name = 'toroidal_magne_evolution' ! ! ----------------------------------------------------------------------- ! @@ -68,19 +62,13 @@ module const_r_mat_4_vector_sph ! ----------------------------------------------------------------------- ! subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, & - & fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, & - & band_vs_poisson, band_vp_evo, band_vt_evo, band_wt_evo) -! - use t_coef_sph_velocity_BCs - use t_coef_fdm2_centre + & fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, & + & band_vs_poisson, band_vp_evo, band_wt_evo) ! use m_ludcmp_band - use set_sph_scalar_matrix_ICB - use set_sph_scalar_matrix_CMB - use cal_inner_core_rotation + use select_sph_r_mat_vort_BC use center_sph_matrices use mat_product_3band_mul - use check_sph_radial_mat ! type(sph_rj_grid), intent(in) :: sph_rj type(fdm_matrices), intent(in) :: r_2nd @@ -93,7 +81,6 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, & real(kind = kreal), intent(in) :: dt ! type(band_matrices_type), intent(inout) :: band_vp_evo - type(band_matrices_type), intent(inout) :: band_vt_evo type(band_matrices_type), intent(inout) :: band_wt_evo type(band_matrices_type), intent(inout) :: band_vs_poisson ! @@ -101,32 +88,24 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, & real(kind = kreal) :: coef_dvt ! ! - band_vt_evo%mat_name = vt_evo_name band_wt_evo%mat_name = wt_evo_name band_vp_evo%mat_name = vp_evo_name ! call alloc_band_mat_sph(ifive, sph_rj, band_vp_evo) - call alloc_band_mat_sph(ithree, sph_rj, band_vt_evo) call alloc_band_mat_sph(ithree, sph_rj, band_wt_evo) call alloc_band_mat_sph(ithree, sph_rj, band_vs_poisson) ! call set_unit_on_diag(band_vp_evo) - call set_unit_on_diag(band_vt_evo) call set_unit_on_diag(band_wt_evo) ! if(fl_prop%coef_diffuse .eq. zero) then coef_dvt = one call set_unit_mat_4_poisson & & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_U%kr_in, sph_bc_U%kr_out, band_vt_evo%mat) - call set_unit_mat_4_poisson & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & & sph_bc_U%kr_in, sph_bc_U%kr_out, band_wt_evo%mat) else coef_dvt = fl_prop%coef_imp * fl_prop%coef_diffuse * dt call set_unit_mat_4_time_evo & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_vt_evo%mat) - call set_unit_mat_4_time_evo & & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_wt_evo%mat) end if ! @@ -137,10 +116,6 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, & call add_vector_poisson_mat_sph & & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, & & g_sph_rj, sph_bc_U%kr_in, sph_bc_U%kr_out, & - & coef_dvt, r_2nd%fdm(2)%dmat, band_vt_evo%mat) - call add_vector_poisson_mat_sph & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, & - & g_sph_rj, sph_bc_U%kr_in, sph_bc_U%kr_out, & & coef_dvt, r_2nd%fdm(2)%dmat, band_wt_evo%mat) call add_vector_poisson_mat_sph & & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, & @@ -148,89 +123,12 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, & & one, r_2nd%fdm(2)%dmat, band_vs_poisson%mat) ! ! Boundary condition for ICB -! - if(sph_bc_U%iflag_icb .eq. iflag_sph_fill_center) then - call add_vector_poisson_mat_center & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%r_ICB, fdm2_center%dmat_fix_fld, & - & coef_dvt, band_vt_evo%mat) - call add_vector_poisson_mat_center & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%r_ICB, fdm2_center%dmat_fix_fld, & - & coef_dvt, band_wt_evo%mat) - call add_vector_poisson_mat_center & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%r_ICB, fdm2_center%dmat_fix_fld, & - & one, band_vs_poisson%mat) - else - call add_fix_flux_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_in, sph_bc_U%r_ICB, sph_bc_U%fdm2_fix_dr_ICB, & - & coef_dvt, band_wt_evo%mat) -! - if(sph_bc_U%iflag_icb .eq. iflag_free_slip) then - call add_fix_flux_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_in, sph_bc_U%r_ICB, & - & bc_fdms_U%fdm2_free_ICB%dmat_vt, & - & coef_dvt, band_vt_evo%mat) - call add_fix_flux_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_in, sph_bc_U%r_ICB, & - & bc_fdms_U%fdm2_free_ICB%dmat_vp, & - & one, band_vs_poisson%mat) -! else if(sph_bc_U%iflag_icb .eq. iflag_evolve_field) then -! else if(sph_bc_U%iflag_icb .eq. iflag_fixed_field) then -! else if(sph_bc_U%iflag_icb .eq. iflag_fixed_velo) then - else - call set_fix_fld_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_U%kr_in, band_vt_evo%mat) - call add_fix_flux_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_in, sph_bc_U%r_ICB, sph_bc_U%fdm2_fix_dr_ICB, & - & one, band_vs_poisson%mat) - end if - end if -! -! Overwrite rotation for inner core -! - if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then - call set_icore_viscous_matrix & - & (sph_bc_U%kr_in, bc_fdms_U%fdm1_fix_fld_ICB, & - & dt, sph_rj, fl_prop, band_vt_evo) - end if + call sel_sph_r_mat_vort_2step_ICB(sph_rj, sph_bc_U, bc_fdms_U, & + & fdm2_center, g_sph_rj, coef_dvt, band_vs_poisson, band_wt_evo) ! ! Boundary condition for CMB -! - call add_fix_flux_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_out, sph_bc_U%r_CMB, sph_bc_U%fdm2_fix_dr_CMB, & - & coef_dvt, band_wt_evo%mat) -! - if(sph_bc_U%iflag_cmb .eq. iflag_free_slip) then - call add_fix_flux_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_out, sph_bc_U%r_CMB, & - & bc_fdms_U%fdm2_free_CMB%dmat_vt, & - & coef_dvt, band_vt_evo%mat) - call add_fix_flux_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_out, sph_bc_U%r_CMB, & - & bc_fdms_U%fdm2_free_CMB%dmat_vp, & - & one, band_vs_poisson%mat) -! else if(sph_bc_U%iflag_cmb .eq. iflag_evolve_field) then -! else if(sph_bc_U%iflag_cmb .eq. iflag_fixed_field) then -! else if(sph_bc_U%iflag_cmb .eq. iflag_fixed_velo) then - else - call set_fix_fld_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_U%kr_out, band_vt_evo%mat) - call add_fix_flux_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_U%kr_out, sph_bc_U%r_CMB, sph_bc_U%fdm2_fix_dr_CMB, & - & one, band_vs_poisson%mat) - end if + call sel_sph_r_mat_vort_2step_CMB(sph_rj, sph_bc_U, bc_fdms_U, & + & g_sph_rj, coef_dvt, band_vs_poisson, band_wt_evo) ! ! call cal_mat_product_3band_mul & @@ -241,160 +139,146 @@ subroutine const_radial_mat_vort_2step(dt, sph_rj, r_2nd, & call ludcmp_5band_mul_t & & (np_smp, sph_rj%istack_rj_j_smp, band_vp_evo) call ludcmp_3band_mul_t & - & (np_smp, sph_rj%istack_rj_j_smp, band_vt_evo) - call ludcmp_3band_mul_t & & (np_smp, sph_rj%istack_rj_j_smp, band_wt_evo) call ludcmp_3band_mul_t & & (np_smp, sph_rj%istack_rj_j_smp, band_vs_poisson) -! - if(i_debug .eq. iflag_full_msg) then - call check_radial_band_mat(my_rank, sph_rj, band_vt_evo) - call check_radial_band_mat(my_rank, sph_rj, band_wt_evo) - call check_radial_band_mat(my_rank, sph_rj, band_vp_evo) - end if -! -! do j = 1, sph_rj%nidx_rj(2) -! do k = 1, sph_rj%nidx_rj(1) -! band_vp_evo%det(j) & -! & = band_vp_evo%det(j) * band_vp_evo%lu(5,k,j) -! band_vt_evo%det(j) & -! & = band_vt_evo%det(j) * band_vt_evo%lu(3,k,j) -! end do -! write(my_rank+60,*) 'det vp', j, & -! & band_vp_evo%det(j), band_vt_evo%det(j) -! end do ! end subroutine const_radial_mat_vort_2step ! ! ----------------------------------------------------------------------- ! - subroutine const_radial_mat_4_magne_sph & - & (dt, sph_rj, r_2nd, cd_prop, sph_bc_B, fdm2_center, & - & g_sph_rj, band_bp_evo, band_bt_evo) + subroutine const_radial_mat_toroidal_flow(dt, sph_rj, r_2nd, & + & fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, & + & band_vt_evo) ! - use t_coef_fdm2_centre - use set_sph_scalar_matrix_ICB - use set_sph_scalar_matrix_CMB + use m_ludcmp_band + use select_sph_r_mat_vort_BC + use cal_inner_core_rotation use center_sph_matrices - use check_sph_radial_mat + use mat_product_3band_mul ! type(sph_rj_grid), intent(in) :: sph_rj type(fdm_matrices), intent(in) :: r_2nd - type(conductive_property), intent(in) :: cd_prop - type(sph_boundary_type), intent(in) :: sph_bc_B + type(fluid_property), intent(in) :: fl_prop + type(sph_boundary_type), intent(in) :: sph_bc_U + type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U type(fdm2_center_mat), intent(in) :: fdm2_center ! real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) real(kind = kreal), intent(in) :: dt ! - type(band_matrices_type), intent(inout) :: band_bp_evo - type(band_matrices_type), intent(inout) :: band_bt_evo -! - real(kind = kreal) :: coef_dbt -! + type(band_matrices_type), intent(inout) :: band_vt_evo ! - band_bp_evo%mat_name = bp_evo_name - band_bt_evo%mat_name = bt_evo_name +! integer(kind = kint) :: j + real(kind = kreal) :: coef_dvt ! - call alloc_band_mat_sph(ithree, sph_rj, band_bp_evo) - call alloc_band_mat_sph(ithree, sph_rj, band_bt_evo) ! - call set_unit_on_diag(band_bp_evo) - call set_unit_on_diag(band_bt_evo) + band_vt_evo%mat_name = vt_evo_name + call alloc_band_mat_sph(ithree, sph_rj, band_vt_evo) + call set_unit_on_diag(band_vt_evo) ! - if(cd_prop%coef_diffuse .eq. zero) then - coef_dbt = one - call set_unit_mat_4_poisson & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_in, sph_bc_B%kr_out, band_bp_evo%mat) + if(fl_prop%coef_diffuse .eq. zero) then + coef_dvt = one call set_unit_mat_4_poisson & & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_in, sph_bc_B%kr_out, band_bt_evo%mat) + & sph_bc_U%kr_in, sph_bc_U%kr_out, band_vt_evo%mat) else - coef_dbt = cd_prop%coef_imp * cd_prop%coef_diffuse * dt - call set_unit_mat_4_time_evo & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_bp_evo%mat) + coef_dvt = fl_prop%coef_imp * fl_prop%coef_diffuse * dt call set_unit_mat_4_time_evo & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_bt_evo%mat) + & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), band_vt_evo%mat) end if +! ! call add_vector_poisson_mat_sph & & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, & - & g_sph_rj, sph_bc_B%kr_in, sph_bc_B%kr_out, & - & coef_dbt, r_2nd%fdm(2)%dmat, band_bp_evo%mat) - call add_vector_poisson_mat_sph & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj, & - & g_sph_rj, sph_bc_B%kr_in, sph_bc_B%kr_out, & - & coef_dbt, r_2nd%fdm(2)%dmat, band_bt_evo%mat) -! -! - if(sph_bc_B%iflag_icb .eq. iflag_sph_fill_center) then - call add_vector_poisson_mat_center & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_B%r_ICB, fdm2_center%dmat_fix_fld, & - & coef_dbt, band_bp_evo%mat) - call add_vector_poisson_mat_center & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_B%r_ICB, fdm2_center%dmat_fix_fld, & - & coef_dbt, band_bt_evo%mat) - else if(sph_bc_B%iflag_icb .eq. iflag_radial_magne) then - call add_fix_flux_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_B%kr_in, sph_bc_B%r_ICB, sph_bc_B%fdm2_fix_dr_ICB, & - & coef_dbt, band_bp_evo%mat) - call set_fix_fld_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_in, band_bt_evo%mat) - else if(sph_bc_B%iflag_icb .eq. iflag_evolve_field) then - call set_fix_fld_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_in, band_bp_evo%mat) - call set_fix_fld_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_in, band_bt_evo%mat) -! else if(sph_bc_B%iflag_icb .eq. iflag_sph_insulator) then - else - call set_ins_magne_icb_rmat_sph & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_B%kr_in, sph_bc_B%r_ICB, sph_bc_B%fdm2_fix_dr_ICB, & - & coef_dbt, band_bp_evo%mat) - call set_fix_fld_icb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_in, band_bt_evo%mat) - end if + & g_sph_rj, sph_bc_U%kr_in, sph_bc_U%kr_out, & + & coef_dvt, r_2nd%fdm(2)%dmat, band_vt_evo%mat) ! - if(sph_bc_B%iflag_cmb .eq. iflag_radial_magne) then - call add_fix_flux_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_B%kr_out, sph_bc_B%r_CMB, sph_bc_B%fdm2_fix_dr_CMB, & - & coef_dbt, band_bp_evo%mat) - else if(sph_bc_B%iflag_cmb .eq. iflag_evolve_field) then - call set_fix_fld_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_out, band_bp_evo%mat) -! else if(sph_bc_B%iflag_cmb .eq. iflag_sph_insulator) then - else - call set_ins_magne_cmb_rmat_sph & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), g_sph_rj, & - & sph_bc_B%kr_out, sph_bc_B%r_CMB, sph_bc_B%fdm2_fix_dr_CMB, & - & coef_dbt, band_bp_evo%mat) +! Boundary condition for ICB + call sel_sph_r_mat_tor_flow_ICB(sph_rj, sph_bc_U, & + & bc_fdms_U, fdm2_center, g_sph_rj, coef_dvt, & + & band_vt_evo) +! +! Overwrite rotation of inner core for degree 1 + if(sph_bc_U%iflag_icb .eq. iflag_rotatable_ic) then + call set_icore_viscous_matrix & + & (sph_bc_U%kr_in, bc_fdms_U%fdm1_fix_fld_ICB, & + & dt, sph_rj, fl_prop, band_vt_evo) end if - call set_fix_fld_cmb_poisson_mat & - & (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), & - & sph_bc_B%kr_out, band_bt_evo%mat) ! +! Boundary condition for CMB + call sel_sph_r_mat_tor_flow_CMB(sph_rj, sph_bc_U, & + & bc_fdms_U, g_sph_rj, coef_dvt, band_vt_evo) ! +! LU decomposition call ludcmp_3band_mul_t & - & (np_smp, sph_rj%istack_rj_j_smp, band_bp_evo) - call ludcmp_3band_mul_t & - & (np_smp, sph_rj%istack_rj_j_smp, band_bt_evo) + & (np_smp, sph_rj%istack_rj_j_smp, band_vt_evo) +! + end subroutine const_radial_mat_toroidal_flow +! +! ----------------------------------------------------------------------- +! + subroutine const_radial_mat7_vpol_press(dt, sph_rj, r_2nd, & + & fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, & + & band7_vsp_evo) ! - if(i_debug .eq. iflag_full_msg) then - call check_radial_band_mat(my_rank, sph_rj, band_bp_evo) - call check_radial_band_mat(my_rank, sph_rj, band_bt_evo) + type(sph_rj_grid), intent(in) :: sph_rj + type(fdm_matrices), intent(in) :: r_2nd + type(fluid_property), intent(in) :: fl_prop + type(sph_boundary_type), intent(in) :: sph_bc_U + type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U + type(fdm2_center_mat), intent(in) :: fdm2_center +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) + real(kind = kreal), intent(in) :: dt +! + type(band_matrices_type), intent(inout) :: band7_vsp_evo +! +! integer(kind = kint) :: j + real(kind = kreal) :: coef_dvt +! +! + band7_vsp_evo%mat_name = vsp_evo_name + call alloc_band_matrices_type(iseven, (2*sph_rj%nidx_rj(1)), & + & sph_rj%nidx_rj(2), band7_vsp_evo) + call set_unit_on_diag(band7_vsp_evo) +! + end subroutine const_radial_mat7_vpol_press +! +! ----------------------------------------------------------------------- +! + subroutine const_radial_mat9_vpol_press(dt, sph_rj, r_2nd, & + & fl_prop, sph_bc_U, bc_fdms_U, fdm2_center, g_sph_rj, & + & band9_vsp_evo) +! + type(sph_rj_grid), intent(in) :: sph_rj + type(fdm_matrices), intent(in) :: r_2nd + type(fluid_property), intent(in) :: fl_prop + type(sph_boundary_type), intent(in) :: sph_bc_U + type(velocity_boundary_FDMs), intent(in) :: bc_fdms_U + type(fdm2_center_mat), intent(in) :: fdm2_center +! + real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13) + real(kind = kreal), intent(in) :: dt +! + type(band_matrices_type), intent(inout) :: band9_vsp_evo +! +! integer(kind = kint) :: j + real(kind = kreal) :: coef_dvt +! +! + if(fl_prop%coef_diffuse .eq. zero) then + coef_dvt = one + else + coef_dvt = fl_prop%coef_imp * fl_prop%coef_diffuse * dt end if ! - end subroutine const_radial_mat_4_magne_sph + band9_vsp_evo%mat_name = vsp_evo_name + call alloc_band_matrices_type(inine, (2*sph_rj%nidx_rj(1)), & + & sph_rj%nidx_rj(2), band9_vsp_evo) + call set_unit_on_diag(band9_vsp_evo) +! + end subroutine const_radial_mat9_vpol_press ! ! ----------------------------------------------------------------------- ! diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_w_center_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_w_center_sph.f90 index 25e26565..a2ad25af 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_w_center_sph.f90 +++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_w_center_sph.f90 @@ -149,7 +149,7 @@ subroutine const_rmat_poisson00_sph(mat_name, sph_rj, & & poisson_mat%n_vect, poisson_mat%n_comp, poisson_mat%mat, & & band_s00_poisson) end if - call check_center_band_matrix(6, sph_rj, band_s00_poisson) + call check_center_band_matrix(50+my_rank, sph_rj, band_s00_poisson) ! end subroutine const_rmat_poisson00_sph ! @@ -201,7 +201,7 @@ subroutine const_rmat_press00_sph & call ludcmp_3band_ctr(band_p00_poisson) ! if(i_debug .ne. iflag_full_msg) return - call check_center_band_matrix(my_rank, sph_rj, band_p00_poisson) + call check_center_band_matrix(50+my_rank, sph_rj, band_p00_poisson) ! end subroutine const_rmat_press00_sph ! @@ -250,7 +250,7 @@ subroutine const_rmat_scalar00_sph(sph_rj, fdm2_center, & call ludcmp_3band_ctr(band_s00_evo) ! if(i_debug .ne. iflag_full_msg) return - call check_center_band_matrix(my_rank, sph_rj, band_s00_evo) + call check_center_band_matrix(50+my_rank, sph_rj, band_s00_evo) ! end subroutine const_rmat_scalar00_sph ! @@ -277,7 +277,7 @@ subroutine copy_radial_mat_scalar00_sph & call ludcmp_3band_ctr(band_s00_poisson) ! if(i_debug .ne. iflag_full_msg) return - call check_center_band_matrix(my_rank, sph_rj, band_s00_poisson) + call check_center_band_matrix(50+my_rank, sph_rj, band_s00_poisson) ! end subroutine copy_radial_mat_scalar00_sph ! diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/const_radial_mat_4_sph.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/const_radial_mat_4_sph.f90 index db30edff..5d0aa5e7 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/const_radial_mat_4_sph.f90 +++ b/src/Fortran_libraries/MHD_src/sph_MHD/const_radial_mat_4_sph.f90 @@ -77,9 +77,11 @@ subroutine const_radial_mat_sph_mhd(dt, MHD_prop, sph_MHD_bc, & ! type(MHD_radial_matrices), intent(inout) :: sph_MHD_mat ! + integer(kind = kint) :: id_file ! + id_file = 50 + my_rank call const_radial_matrices_sph & - & (dt, sph%sph_params, sph%sph_rj, r_2nd, & + & (id_file, dt, sph%sph_params, sph%sph_rj, r_2nd, & & MHD_prop, sph_MHD_bc, leg%g_sph_rj, sph_MHD_mat) ! call const_radial_mat_sph_w_center & @@ -104,6 +106,9 @@ subroutine const_radial_mat_sph_snap(MHD_prop, sph_MHD_bc, & ! type(MHD_radial_matrices), intent(inout) :: sph_MHD_mat ! + integer(kind = kint) :: id_file +! + id_file = 50 + my_rank ! if(MHD_prop%fl_prop%iflag_scheme .lt. id_Crank_nicolson) return if(iflag_debug .gt. 0) & @@ -112,6 +117,10 @@ subroutine const_radial_mat_sph_snap(MHD_prop, sph_MHD_bc, & & (sph_rj, r_2nd, MHD_prop%fl_prop, & & sph_MHD_bc%sph_bc_U, sph_MHD_bc%fdm2_center, & & leg%g_sph_rj, sph_MHD_mat%band_p_poisson) + if(i_debug .eq. iflag_full_msg) then + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_p_poisson) + end if ! if(sph_rj%inod_rj_center .eq. 0) return ! @@ -125,12 +134,15 @@ end subroutine const_radial_mat_sph_snap ! ----------------------------------------------------------------------- ! ----------------------------------------------------------------------- ! - subroutine const_radial_matrices_sph(dt, sph_params, sph_rj, & - & r_2nd, MHD_prop, sph_MHD_bc, g_sph_rj, sph_MHD_mat) + subroutine const_radial_matrices_sph & + & (id_file, dt, sph_params, sph_rj, r_2nd, MHD_prop, & + & sph_MHD_bc, g_sph_rj, sph_MHD_mat) ! use const_r_mat_4_scalar_sph use const_r_mat_4_vector_sph + use const_r_mat_4_magnetic_sph ! + integer(kind = kint), intent(in) :: id_file type(sph_shell_parameters), intent(in) :: sph_params type(sph_rj_grid), intent(in) :: sph_rj type(fdm_matrices), intent(in) :: r_2nd @@ -145,17 +157,47 @@ subroutine const_radial_matrices_sph(dt, sph_params, sph_rj, & ! if(MHD_prop%fl_prop%iflag_scheme .ge. id_Crank_nicolson) then if(iflag_debug .gt. 0) & + & write(*,*) 'const_radial_mat_toroidal_flow' + call const_radial_mat_toroidal_flow & + & (dt, sph_rj, r_2nd, MHD_prop%fl_prop, & + & sph_MHD_bc%sph_bc_U, sph_MHD_bc%bc_fdms_U, & + & sph_MHD_bc%fdm2_center, g_sph_rj, sph_MHD_mat%band_vt_evo) + if(iflag_debug .gt. 0) & & write(*,*) 'const_radial_mat_vort_2step' call const_radial_mat_vort_2step & & (dt, sph_rj, r_2nd, MHD_prop%fl_prop, & & sph_MHD_bc%sph_bc_U, sph_MHD_bc%bc_fdms_U, & & sph_MHD_bc%fdm2_center, g_sph_rj, & & sph_MHD_mat%band_vs_poisson, sph_MHD_mat%band_vp_evo, & - & sph_MHD_mat%band_vt_evo, sph_MHD_mat%band_wt_evo) + & sph_MHD_mat%band_wt_evo) +! + if(iflag_debug .gt. 0) & + & write(*,*) 'const_radial_mat_4_press_sph' call const_radial_mat_4_press_sph & & (sph_rj, r_2nd, MHD_prop%fl_prop, & & sph_MHD_bc%sph_bc_U, sph_MHD_bc%fdm2_center, & & g_sph_rj, sph_MHD_mat%band_p_poisson) +! + if(i_debug .eq. iflag_full_msg) then + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_wt_evo) + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_vp_evo) + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_vt_evo) + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_p_poisson) + end if +! +! do j = 1, sph_rj%nidx_rj(2) +! do k = 1, sph_rj%nidx_rj(1) +! sph_MHD_mat%band_vp_evo%det(j) & +! & = sph_MHD_mat%band_vp_evo%det(j) & +! & * sph_MHD_mat%band_vp_evo%lu(5,k,j) +! end do +! write(my_rank+60,*) 'det vp', j, & +! & sph_MHD_mat%band_vp_evo%det(j) +! end do end if ! if(MHD_prop%cd_prop%iflag_Bevo_scheme & @@ -164,6 +206,13 @@ subroutine const_radial_matrices_sph(dt, sph_params, sph_rj, & & (dt, sph_rj, r_2nd, MHD_prop%cd_prop, & & sph_MHD_bc%sph_bc_B, sph_MHD_bc%fdm2_center, & & g_sph_rj, sph_MHD_mat%band_bp_evo, sph_MHD_mat%band_bt_evo) +! + if(i_debug .eq. iflag_full_msg) then + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_bp_evo) + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_bt_evo) + end if end if ! if(MHD_prop%ht_prop%iflag_scheme .ge. id_Crank_nicolson) then @@ -172,6 +221,10 @@ subroutine const_radial_matrices_sph(dt, sph_params, sph_rj, & & sph_params, sph_rj, r_2nd, MHD_prop%ht_prop, & & sph_MHD_bc%sph_bc_T, sph_MHD_bc%fdm2_center, & & g_sph_rj, sph_MHD_mat%band_temp_evo) + if(i_debug .eq. iflag_full_msg) then + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_temp_evo) + end if end if ! if(MHD_prop%cp_prop%iflag_scheme .ge. id_Crank_nicolson) then @@ -180,6 +233,10 @@ subroutine const_radial_matrices_sph(dt, sph_params, sph_rj, & & sph_params, sph_rj, r_2nd, MHD_prop%cp_prop, & & sph_MHD_bc%sph_bc_C, sph_MHD_bc%fdm2_center, & & g_sph_rj, sph_MHD_mat%band_comp_evo) + if(i_debug .eq. iflag_full_msg) then + call check_radial_band_mat(id_file, sph_rj, & + & sph_MHD_mat%band_comp_evo) + end if end if ! end subroutine const_radial_matrices_sph diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/solve_sph_fluid_crank.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/solve_sph_fluid_crank.f90 index f7ed370c..cfd6c19f 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/solve_sph_fluid_crank.f90 +++ b/src/Fortran_libraries/MHD_src/sph_MHD/solve_sph_fluid_crank.f90 @@ -111,12 +111,12 @@ subroutine solve_pressure_by_div_v(sph_rj, band_p_poisson, & open(id_file,file='ave_press_test.txt') write(id_file,*) 'Matrix for average pressure' call check_single_radial_3band_mat & - & (id_offset, sph_rj%nidx_rj(1), sph_rj%radius_1d_rj_r, & + & (id_file, sph_rj%nidx_rj(1), sph_rj%radius_1d_rj_r, & & band_p_poisson%mat(1:band_p_poisson%n_band, & & 1:band_p_poisson%n_vect,j)) write(id_file,*) 'LU decomposition for average pressure' call check_single_radial_5band_mat & - & (id_offset, sph_rj%nidx_rj(1), sph_rj%radius_1d_rj_r, & + & (id_file, sph_rj%nidx_rj(1), sph_rj%radius_1d_rj_r, & & band_p_poisson%lu(1:band_p_poisson%n_band_lu, & & 1:band_p_poisson%n_vect,j)) write(id_file,*) 'RHS for average pressure' @@ -192,7 +192,7 @@ subroutine solve_scalar_sph_crank & ! j = find_local_sph_address(sph_rj, 0,0) ! if(j.gt.0) then ! write(*,*) 'matrix' -! call check_single_radial_3band_mat(my_rank, sph_rj%nidx_rj(1), & +! call check_single_radial_3band_mat(6, sph_rj%nidx_rj(1), & ! & sph_rj%radius_1d_rj_r, band_s_evo%mat(1,1,j)) ! end if ! @@ -210,7 +210,7 @@ subroutine solve_scalar_sph_crank & if(sph_rj%inod_rj_center .eq. 0) return ! ! write(50+my_rank,*) 'matrix for l=m=0' -! call check_single_radial_3band_mat(my_rank, sph_rj%nidx_rj(1), & +! call check_single_radial_3band_mat(6, sph_rj%nidx_rj(1), & ! & sph_rj%radius_1d_rj_r, band_s00_evo%mat) ! call lubksb_3band_ctr(band_s00_evo, sol_00) diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/t_radial_matrices_sph_MHD.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/t_radial_matrices_sph_MHD.f90 index 1fbd6095..0fc104f7 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/t_radial_matrices_sph_MHD.f90 +++ b/src/Fortran_libraries/MHD_src/sph_MHD/t_radial_matrices_sph_MHD.f90 @@ -12,8 +12,10 @@ !! !! subroutine allocate_press_vpol_mat_sph(sph_rj, sph_MHD_mat) !! subroutine check_velocity_matrices_sph & -!! & (id_rank, sph_rj, sph_MHD_mat) -!! type(MHD_radial_matrices), intent(inout) :: sph_MHD_mat +!! & (id_file, sph_rj, sph_MHD_mat) +!! integer(kind = kint), intent(in) :: id_file +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(MHD_radial_matrices), intent(in) :: sph_MHD_mat !!@endverbatim !! module t_radial_matrices_sph_MHD @@ -38,7 +40,7 @@ module t_radial_matrices_sph_MHD !> Structure of band matrices for pressure poisson type(band_matrices_type) :: band_p_poisson ! -!> Structure of band matrices for poloidal velocity +!> Structure of band matrices for poloidal velocity and pressure type(band_matrices_type) :: band_vsp_evo ! !> Structure of band matrices for poloidal magnetic field @@ -143,11 +145,11 @@ end subroutine allocate_press_vpol_mat_sph ! ----------------------------------------------------------------------- ! subroutine check_velocity_matrices_sph & - & (id_rank, sph_rj, sph_MHD_mat) + & (id_file, sph_rj, sph_MHD_mat) ! use check_sph_radial_mat ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file type(sph_rj_grid), intent(in) :: sph_rj type(MHD_radial_matrices), intent(in) :: sph_MHD_mat ! @@ -161,11 +163,11 @@ subroutine check_velocity_matrices_sph & end do ! call check_radial_7band_mat & - & (id_rank, (2*sph_rj%nidx_rj(1)), sph_rj%nidx_rj(2), & + & (id_file, (2*sph_rj%nidx_rj(1)), sph_rj%nidx_rj(2), & & sph_rj%idx_gl_1d_rj_j, rr, sph_MHD_mat%band_vsp_evo%mat) ! call check_radial_band_mat & - & (id_rank, sph_rj, sph_MHD_mat%band_vt_evo) + & (id_file, sph_rj, sph_MHD_mat%band_vt_evo) ! end subroutine check_velocity_matrices_sph ! diff --git a/src/Fortran_libraries/MHD_src/sph_MHD/t_sph_center_matrix.f90 b/src/Fortran_libraries/MHD_src/sph_MHD/t_sph_center_matrix.f90 index b6fc2454..fbf938c0 100644 --- a/src/Fortran_libraries/MHD_src/sph_MHD/t_sph_center_matrix.f90 +++ b/src/Fortran_libraries/MHD_src/sph_MHD/t_sph_center_matrix.f90 @@ -21,6 +21,9 @@ !! type(band_matrix_type), intent(inout) :: smat !! !! subroutine check_center_band_matrix(id_rank, sph_rj, smat) +!! integer(kind = kint), intent(in) :: id_file +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(band_matrix_type), intent(in) :: smat !!@endverbatim ! module t_sph_center_matrix @@ -183,18 +186,18 @@ end subroutine lubksb_7band_ctr ! ---------------------------------------------------------------------- ! ---------------------------------------------------------------------- ! - subroutine check_center_band_matrix(id_rank, sph_rj, smat) + subroutine check_center_band_matrix(id_file, sph_rj, smat) ! use t_spheric_rj_data use check_sph_radial_mat ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file type(sph_rj_grid), intent(in) :: sph_rj type(band_matrix_type), intent(in) :: smat ! ! - write(50+id_rank,'(a)') 'Matrix for ', trim(smat%mat_name) - call check_radial_3band_mat_w_ctr(id_rank, smat%n_vect, & + write(id_file,'(a)') 'Matrix for ', trim(smat%mat_name) + call check_radial_3band_mat_w_ctr(id_file, smat%n_vect, & & sph_rj%radius_1d_rj_r, smat%mat) ! end subroutine check_center_band_matrix diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/check_sph_radial_mat.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/check_sph_radial_mat.f90 index dee2d80c..e4e5ee4a 100644 --- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/check_sph_radial_mat.f90 +++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/check_sph_radial_mat.f90 @@ -7,26 +7,50 @@ !>@brief Output band matrix data to check !! !!@verbatim -!! subroutine check_radial_3band_mat_w_ctr(id_rank, nri, rr, mat00) +!! subroutine check_radial_3band_mat_w_ctr(id_file, nri, rr, mat00) +!! integer(kind = kint), intent(in) :: id_file +!! integer(kind = kint), intent(in) :: nri +!! real(kind = kreal), intent(in) :: rr(nri) +!! real(kind = kreal), intent(in) :: mat00(3,0:nri) !! -!! subroutine check_radial_3band_mat(id_rank, nri, jmax, j_sph, & -!! & rr, mat) -!! subroutine check_radial_5band_mat(id_rank, nri, jmax, j_sph, & -!! & rr, mat) -!! subroutine check_radial_7band_mat(id_rank, nri, jmax, j_sph, & -!! & rr, mat) +!! subroutine check_radial_3band_mat(id_file, nri, jmax, & +!! & j_sph, rr, mat) +!! integer, intent(in) :: id_rank +!! integer(kind = kint), intent(in) :: id_file +!! integer(kind = kint), intent(in) :: nri, jmax +!! integer(kind = kint), intent(in) :: j_sph(jmax,3) +!! real(kind = kreal), intent(in) :: rr(nri) +!! real(kind = kreal), intent(in) :: mat(3,nri,jmax) +!! subroutine check_radial_5band_mat(id_file, nri, jmax, & +!! & j_sph, rr, mat) +!! real(kind = kreal), intent(in) :: mat(5,nri,jmax) +!! subroutine check_radial_7band_mat(id_file, nri, jmax, & +!! & j_sph, rr, mat) +!! real(kind = kreal), intent(in) :: mat(7,nri,jmax) +!! subroutine check_radial_9band_mat(id_file, nri, jmax, & +!! & j_sph, rr, mat) +!! real(kind = kreal), intent(in) :: mat(9,nri,jmax) !! -!! subroutine check_single_radial_3band_mat(id_rank, nri, rr, mat) -!! subroutine check_single_radial_5band_mat(id_rank, nri, rr, mat) -!! subroutine check_single_radial_7band_mat(id_rank, nri, rr, mat) +!! subroutine check_single_radial_3band_mat(id_file, nri, rr, mat) +!! integer(kind = kint), intent(in) :: id_file +!! integer(kind = kint), intent(in) :: nri, jmax +!! integer(kind = kint), intent(in) :: j_sph(jmax,3) +!! real(kind = kreal), intent(in) :: rr(nri) +!! real(kind = kreal), intent(in) :: mat(3,nri) +!! subroutine check_single_radial_5band_mat(id_file, nri, rr, mat) +!! real(kind = kreal), intent(in) :: mat(5,nri) +!! subroutine check_single_radial_7band_mat(id_file, nri, rr, mat) +!! real(kind = kreal), intent(in) :: mat(7,nri) +!! subroutine check_single_radial_9band_mat(id_file, nri, rr, mat) +!! real(kind = kreal), intent(in) :: mat(9,nri) !!@endverbatim !! -!!@n @param id_rank radius to three next points of ICB -!!@n @param nri radius to three next points of ICB -!!@n @param jmax radius to three next points of ICB -!!@n @param j_sph radius to three next points of ICB -!!@n @param rr radius to three next points of ICB -!!@n @param mat radius to three next points of ICB +!!@n @param id_file File ID +!!@n @param nri Number of radial points +!!@n @param jmax Number of modes +!!@n @param j_sph Modes +!!@n @param rr radius +!!@n @param mat Band matrix ! module check_sph_radial_mat ! @@ -41,9 +65,9 @@ module check_sph_radial_mat ! ! ----------------------------------------------------------------------- ! - subroutine check_radial_3band_mat_w_ctr(id_rank, nri, rr, mat00) + subroutine check_radial_3band_mat_w_ctr(id_file, nri, rr, mat00) ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file integer(kind = kint), intent(in) :: nri real(kind = kreal), intent(in) :: rr(nri) real(kind = kreal), intent(in) :: mat00(3,0:nri) @@ -51,14 +75,14 @@ subroutine check_radial_3band_mat_w_ctr(id_rank, nri, rr, mat00) integer(kind = kint) :: k ! ! - write(50+id_rank,'(a)') 'k, r, a(k,k-1), a(k,k), a(k,k+1)' - write(50+id_rank,'(i6,1p4E25.15e3)') izero, zero, & + write(id_file,'(a)') 'k, r, a(k,k-1), a(k,k), a(k,k+1)' + write(id_file,'(i6,1p4E25.15e3)') izero, zero, & & -1.0d30, mat00(2,0), mat00(1,1) do k = 1, nri-1 - write(50+id_rank,'(i6,1p4E25.15e3)') k, rr(k), & + write(id_file,'(i6,1p4E25.15e3)') k, rr(k), & & mat00(3,k-1), mat00(2,k), mat00(1,k+1) end do - write(50+id_rank,'(i6,1p4E25.15e3)') nri, rr(nri), & + write(id_file,'(i6,1p4E25.15e3)') nri, rr(nri), & & mat00(3,nri-1), mat00(2,nri), 1.0d30 ! end subroutine check_radial_3band_mat_w_ctr @@ -66,10 +90,10 @@ end subroutine check_radial_3band_mat_w_ctr ! ----------------------------------------------------------------------- ! ----------------------------------------------------------------------- ! - subroutine check_radial_3band_mat(id_rank, nri, jmax, j_sph, & - & rr, mat) + subroutine check_radial_3band_mat(id_file, nri, jmax, & + & j_sph, rr, mat) ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file integer(kind = kint), intent(in) :: nri, jmax integer(kind = kint), intent(in) :: j_sph(jmax,3) real(kind = kreal), intent(in) :: rr(nri) @@ -79,20 +103,20 @@ subroutine check_radial_3band_mat(id_rank, nri, jmax, j_sph, & ! ! do j = 1, jmax - write(50+id_rank,'(a,4i6)') '(j, global_j, l, m): ', & + write(id_file,'(a,4i6)') '(j, global_j, l, m): ', & & j, j_sph(jmax,1:3) call check_single_radial_3band_mat & - & (id_rank, nri, rr, mat(1,1,j)) + & (id_file, nri, rr, mat(1,1,j)) end do ! end subroutine check_radial_3band_mat ! ! ----------------------------------------------------------------------- ! - subroutine check_radial_5band_mat(id_rank, nri, jmax, j_sph, & - & rr, mat) + subroutine check_radial_5band_mat(id_file, nri, jmax, & + & j_sph, rr, mat) ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file integer(kind = kint), intent(in) :: nri, jmax integer(kind = kint), intent(in) :: j_sph(jmax,3) real(kind = kreal), intent(in) :: rr(nri) @@ -102,20 +126,20 @@ subroutine check_radial_5band_mat(id_rank, nri, jmax, j_sph, & ! ! do j = 1, jmax - write(50+id_rank,'(a,4i6)') '(j, global_j, l, m): ', & + write(id_file,'(a,4i6)') '(j, global_j, l, m): ', & & j, j_sph(jmax,1:3) - call check_single_radial_5band_mat & - & (id_rank, nri, rr, mat(1,1,j)) + call check_single_radial_5band_mat(id_file, nri, rr, & + & mat(1,1,j)) end do ! end subroutine check_radial_5band_mat ! ! ----------------------------------------------------------------------- ! - subroutine check_radial_7band_mat(id_rank, nri, jmax, j_sph, & - & rr, mat) + subroutine check_radial_7band_mat(id_file, nri, jmax, & + & j_sph, rr, mat) ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file integer(kind = kint), intent(in) :: nri, jmax integer(kind = kint), intent(in) :: j_sph(jmax,3) real(kind = kreal), intent(in) :: rr(nri) @@ -125,20 +149,43 @@ subroutine check_radial_7band_mat(id_rank, nri, jmax, j_sph, & ! ! do j = 1, jmax - write(50+id_rank,'(a,4i6)') '(j, global_j, l, m): ', & + write(id_file,'(a,4i6)') '(j, global_j, l, m): ', & & j, j_sph(jmax,1:3) - call check_single_radial_7band_mat & - & (id_rank, nri, rr, mat(1,1,j)) + call check_single_radial_7band_mat(id_file, nri, rr, & + & mat(1,1,j)) end do ! end subroutine check_radial_7band_mat ! ! ----------------------------------------------------------------------- +! + subroutine check_radial_9band_mat(id_file, nri, jmax, & + & j_sph, rr, mat) +! + integer(kind = kint), intent(in) :: id_file + integer(kind = kint), intent(in) :: nri, jmax + integer(kind = kint), intent(in) :: j_sph(jmax,3) + real(kind = kreal), intent(in) :: rr(nri) + real(kind = kreal), intent(in) :: mat(9,nri,jmax) +! + integer(kind = kint) :: j +! +! + do j = 1, jmax + write(id_file,'(a,4i6)') '(j, global_j, l, m): ', & + & j, j_sph(jmax,1:3) + call check_single_radial_9band_mat(id_file, nri, rr, & + & mat(1,1,j)) + end do +! + end subroutine check_radial_9band_mat +! +! ----------------------------------------------------------------------- ! ----------------------------------------------------------------------- ! - subroutine check_single_radial_3band_mat(id_rank, nri, rr, mat) + subroutine check_single_radial_3band_mat(id_file, nri, rr, mat) ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file integer(kind = kint), intent(in) :: nri real(kind = kreal), intent(in) :: rr(nri) real(kind = kreal), intent(in) :: mat(3,nri) @@ -146,24 +193,24 @@ subroutine check_single_radial_3band_mat(id_rank, nri, rr, mat) integer(kind = kint) :: k ! ! - write(50+id_rank,'(a)') & + write(id_file,'(a)') & & 'k, r, a(k,k-1), a(k,k), a(k,k+1)' - write(50+id_rank,'(i6,1p4E25.15e3)') ione, rr(1), & + write(id_file,'(i6,1p4E25.15e3)') ione, rr(1), & & -1.0d30, mat(2,1), mat(1,2) do k = 2, nri-1 - write(50+id_rank,'(i6,1p4E25.15e3)') k, rr(k), & + write(id_file,'(i6,1p4E25.15e3)') k, rr(k), & & mat(3,k-1), mat(2,k), mat(1,k+1) end do - write(50+id_rank,'(i6,1p4E25.15e3)') nri, rr(nri), & + write(id_file,'(i6,1p4E25.15e3)') nri, rr(nri), & & mat(3,nri-1), mat(2,nri), 1.0d30 ! end subroutine check_single_radial_3band_mat ! ! ----------------------------------------------------------------------- ! - subroutine check_single_radial_5band_mat(id_rank, nri, rr, mat) + subroutine check_single_radial_5band_mat(id_file, nri, rr, mat) ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file integer(kind = kint), intent(in) :: nri real(kind = kreal), intent(in) :: rr(nri) real(kind = kreal), intent(in) :: mat(5,nri) @@ -171,28 +218,28 @@ subroutine check_single_radial_5band_mat(id_rank, nri, rr, mat) integer(kind = kint) :: k, l ! ! - write(50+id_rank,'(a)') & + write(id_file,'(a)') & & 'k, r, a(k,k-2), a(k,k-1), a(k,k), a(k,k+1), a(k,k+2)' - write(50+id_rank,'(i6,1p6E25.15e3)') ione, rr(1), & - & -1.0d30, -1.0d30, (mat(3-l,1+l),l=0,2) - write(50+id_rank,'(i6,1p6E25.15e3)') itwo, rr(2), & - & -1.0d30, (mat(3-l,2+l),l=-1,2) + write(id_file,'(i6,1p6E25.15e3)') ione, rr(1), & + & -1.0d30, -1.0d30, (mat(3-l,1+l),l=0,2) + write(id_file,'(i6,1p6E25.15e3)') itwo, rr(2), & + & -1.0d30, (mat(3-l,2+l),l=-1,2) do k = 3, nri-2 - write(50+id_rank,'(i6,1p6E25.15e3)') k, rr(k), & - & (mat(3-l,k+l),l=-2,2) + write(id_file,'(i6,1p6E25.15e3)') k, rr(k), & + & (mat(3-l,k+l),l=-2,2) end do - write(50+id_rank,'(i6,1p6E25.15e3)') (nri-1), rr(nri-1), & - & (mat(3-l,nri-1+l),l=-2,1), 1.0d30 - write(50+id_rank,'(i6,1p6E25.15e3)') nri, rr(nri), & - & (mat(3-l,nri+l),l=-2,0), 1.0d30, 1.0d30 + write(id_file,'(i6,1p6E25.15e3)') (nri-1), rr(nri-1), & + & (mat(3-l,nri-1+l),l=-2,1), 1.0d30 + write(id_file,'(i6,1p6E25.15e3)') nri, rr(nri), & + & (mat(3-l,nri+l),l=-2,0), 1.0d30, 1.0d30 ! end subroutine check_single_radial_5band_mat ! ! ----------------------------------------------------------------------- ! - subroutine check_single_radial_7band_mat(id_rank, nri, rr, mat) + subroutine check_single_radial_7band_mat(id_file, nri, rr, mat) ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file integer(kind = kint), intent(in) :: nri real(kind = kreal), intent(in) :: rr(nri) real(kind = kreal), intent(in) :: mat(7,nri) @@ -200,29 +247,68 @@ subroutine check_single_radial_7band_mat(id_rank, nri, rr, mat) integer(kind = kint) :: k, l ! ! - write(50+id_rank,'(a,a)') & + write(id_file,'(a,a)') & & 'k, r, a(k,k-3), a(k,k-2), a(k,k-1), a(k,k), ', & & 'a(k,k+1), a(k,k+2), a(k,k+3)' ! - write(50+id_rank,'(i6,1p8E25.15e3)') ione, rr(1), & - & -1.0d30, -1.0d30, -1.0d30, (mat(4-l,1+l),l=0,3) - write(50+id_rank,'(i6,1p8E25.15e3)') itwo, rr(2), & - & -1.0d30, -1.0d30, (mat(4-l,2+l),l=-1,3) - write(50+id_rank,'(i6,1p8E25.15e3)') ithree, rr(3), & - & -1.0d30, (mat(4-l,3+l),l=-2,3) + write(id_file,'(i6,1p8E25.15e3)') ione, rr(1), & + & -1.0d30, -1.0d30, -1.0d30, (mat(4-l,1+l),l=0,3) + write(id_file,'(i6,1p8E25.15e3)') itwo, rr(2), & + & -1.0d30, -1.0d30, (mat(4-l,2+l),l=-1,3) + write(id_file,'(i6,1p8E25.15e3)') ithree, rr(3), & + & -1.0d30, (mat(4-l,3+l),l=-2,3) do k = 4, nri-3 - write(50+id_rank,'(i6,1p8E25.15e3)') k, rr(k), & - & (mat(4-l,k+l),l=-3,3) + write(id_file,'(i6,1p8E25.15e3)') k, rr(k), & + & (mat(4-l,k+l),l=-3,3) end do - write(50+id_rank,'(i6,1p8E25.15e3)') (nri-2), rr(nri-2), & - & (mat(4-l,nri-2+l),l=-3,2), 1.0d30 - write(50+id_rank,'(i6,1p8E25.15e3)') (nri-1), rr(nri-1), & - & (mat(4-l,nri-1+l),l=-3,1), 1.0d30, 1.0d30 - write(50+id_rank,'(i6,1p8E25.15e3)') nri, rr(nri ), & - & (mat(4-l,nri +l),l=-3,0), 1.0d30, 1.0d30, 1.0d30 + write(id_file,'(i6,1p8E25.15e3)') (nri-2), rr(nri-2), & + & (mat(4-l,nri-2+l),l=-3,2), 1.0d30 + write(id_file,'(i6,1p8E25.15e3)') (nri-1), rr(nri-1), & + & (mat(4-l,nri-1+l),l=-3,1), 1.0d30, 1.0d30 + write(id_file,'(i6,1p8E25.15e3)') nri, rr(nri ), & + & (mat(4-l,nri +l),l=-3,0), 1.0d30, 1.0d30, 1.0d30 ! end subroutine check_single_radial_7band_mat ! ! ----------------------------------------------------------------------- +! + subroutine check_single_radial_9band_mat(id_file, nri, rr, mat) +! + integer(kind = kint), intent(in) :: id_file + integer(kind = kint), intent(in) :: nri + real(kind = kreal), intent(in) :: rr(nri) + real(kind = kreal), intent(in) :: mat(9,nri) +! + integer(kind = kint) :: k, l +! +! + write(id_file,'(a,a)') & + & 'k, r, a(k,k-4), a(k,k-3), a(k,k-2), a(k,k-1), a(k,k), ', & + & 'a(k,k+1), a(k,k+2), a(k,k+3), a(k,k+4)' +! + write(id_file,'(i6,1p10E25.15e3)') ione, rr(1), & + & -1.0d30, -1.0d30, -1.0d30, -1.0d30, (mat(5-l,1+l),l= 0,4) + write(id_file,'(i6,1p10E25.15e3)') itwo, rr(2), & + & -1.0d30, -1.0d30, -1.0d30, (mat(5-l,2+l),l=-1,4) + write(id_file,'(i6,1p10E25.15e3)') ithree, rr(3), & + & -1.0d30, -1.0d30, (mat(5-l,3+l),l=-2,4) + write(id_file,'(i6,1p10E25.15e3)') ifour, rr(4), & + & -1.0d30, (mat(5-l,4+l),l=-3,4) + do k = 5, nri-4 + write(id_file,'(i6,1p10E25.15e3)') k, rr(k), & + & (mat(5-l,k+l),l=-4,4) + end do + write(id_file,'(i6,1p10E25.15e3)') (nri-3), rr(nri-3), & + & (mat(5-l,nri-3+l),l=-4,3), 1.0d30 + write(id_file,'(i6,1p10E25.15e3)') (nri-2), rr(nri-2), & + & (mat(5-l,nri-2+l),l=-4,2), 1.0d30, 1.0d30 + write(id_file,'(i6,1p10E25.15e3)') (nri-1), rr(nri-1), & + & (mat(5-l,nri-1+l),l=-4,1), 1.0d30, 1.0d30, 1.0d30 + write(id_file,'(i6,1p10E25.15e3)') nri, rr(nri ), & + & (mat(5-l,nri +l),l=-4,0), 1.0d30, 1.0d30, 1.0d30, 1.0d30 +! + end subroutine check_single_radial_9band_mat +! +! ----------------------------------------------------------------------- ! end module check_sph_radial_mat diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrices.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrices.f90 index 67a27c4d..b6d9f497 100644 --- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrices.f90 +++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrices.f90 @@ -28,7 +28,10 @@ !! type(sph_rj_grid), intent(in) :: sph_rj !! type(band_matrices_type), intent(inout) :: smat !! -!! subroutine check_radial_band_mat(id_rank, sph_rj, smat) +!! subroutine check_radial_band_mat(id_file, sph_rj, smat) +!! integer(kind = kint), intent(in) :: id_file +!! type(sph_rj_grid), intent(in) :: sph_rj +!! type(band_matrices_type), intent(in) :: smat !!@endverbatim ! module t_sph_matrices @@ -250,28 +253,32 @@ end subroutine lubksb_7band_mul_t ! ---------------------------------------------------------------------- ! ---------------------------------------------------------------------- ! - subroutine check_radial_band_mat(id_rank, sph_rj, smat) + subroutine check_radial_band_mat(id_file, sph_rj, smat) ! use t_spheric_rj_data use check_sph_radial_mat ! - integer, intent(in) :: id_rank + integer(kind = kint), intent(in) :: id_file type(sph_rj_grid), intent(in) :: sph_rj type(band_matrices_type), intent(in) :: smat ! ! - write(50+id_rank,'(a,a)') 'Matrix ', trim(smat%mat_name) + write(id_file,'(a,a)') 'Matrix ', trim(smat%mat_name) ! if(smat%n_band .eq. ithree) then - call check_radial_3band_mat(id_rank, smat%n_vect, smat%n_comp, & + call check_radial_3band_mat(id_file, smat%n_vect, smat%n_comp, & & sph_rj%idx_gl_1d_rj_j, sph_rj%radius_1d_rj_r, smat%mat) ! else if(smat%n_band .eq. ifive) then - call check_radial_5band_mat(id_rank, smat%n_vect, smat%n_comp, & + call check_radial_5band_mat(id_file, smat%n_vect, smat%n_comp, & & sph_rj%idx_gl_1d_rj_j, sph_rj%radius_1d_rj_r, smat%mat) ! else if(smat%n_band .eq. iseven) then - call check_radial_7band_mat(id_rank, smat%n_vect, smat%n_comp, & + call check_radial_7band_mat(id_file, smat%n_vect, smat%n_comp, & + & sph_rj%idx_gl_1d_rj_j, sph_rj%radius_1d_rj_r, smat%mat) +! + else if(smat%n_band .eq. inine) then + call check_radial_9band_mat(id_file, smat%n_vect, smat%n_comp, & & sph_rj%idx_gl_1d_rj_j, sph_rj%radius_1d_rj_r, smat%mat) end if ! diff --git a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrix.f90 b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrix.f90 index 4d591098..dbc66562 100644 --- a/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrix.f90 +++ b/src/Fortran_libraries/PARALLEL_src/SPH_SHELL_src/t_sph_matrix.f90 @@ -220,15 +220,19 @@ subroutine check_radial_band_matrix(id_rank, sph_rj, smat) ! ! if(smat%n_band .eq. ithree) then - call check_single_radial_3band_mat(id_rank, smat%n_vect, & + call check_single_radial_3band_mat((50+id_rank), smat%n_vect, & & sph_rj%radius_1d_rj_r, smat%mat) ! else if(smat%n_band .eq. ifive) then - call check_single_radial_5band_mat(id_rank, smat%n_vect, & + call check_single_radial_5band_mat((50+id_rank), smat%n_vect, & & sph_rj%radius_1d_rj_r, smat%mat) ! else if(smat%n_band .eq. iseven) then - call check_single_radial_7band_mat(id_rank, smat%n_vect, & + call check_single_radial_7band_mat((50+id_rank), smat%n_vect, & + & sph_rj%radius_1d_rj_r, smat%mat) +! + else if(smat%n_band .eq. inine) then + call check_single_radial_9band_mat((50+id_rank), smat%n_vect, & & sph_rj%radius_1d_rj_r, smat%mat) end if ! diff --git a/src/Fortran_libraries/SERIAL_src/IO_ZLIB/init_reference_scalar.f90 b/src/Fortran_libraries/SERIAL_src/IO_ZLIB/init_reference_scalar.f90 index bdcf4e55..09d194de 100644 --- a/src/Fortran_libraries/SERIAL_src/IO_ZLIB/init_reference_scalar.f90 +++ b/src/Fortran_libraries/SERIAL_src/IO_ZLIB/init_reference_scalar.f90 @@ -118,7 +118,7 @@ subroutine s_init_reference_scalar & else if(ref_param%iflag_reference & & .eq. id_numerical_solution) then call const_r_mat00_scalar_sph & - & (mat_name, sc_prop%diffusie_reduction_ICB, & + & ((my_rank+50), mat_name, sc_prop%diffusie_reduction_ICB, & & sph_params, sph_rj, r_2nd, sph_bc_S, fdm2_center, & & band_s00_poisson) file_name = add_dat_extension(mat_name)