Add radial diffusivity variation for Crank-Nicolson matrix construction

Author: hiromatsui

src/Fortran_libraries/MHD_src/sph_MHD/const_r_mat_4_scalar_sph.f90

@@ -121,7 +121,8 @@ end subroutine const_radial_mat_4_press_sph
 !
       subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect,   &
      &          dt, sph_params, sph_rj, r_2nd, property,                &
-     &          sph_bc, bcs_S, fdm2_center, g_sph_rj, band_s_evo)
+     &          sph_bc, bcs_S, fdm2_center, g_sph_rj, k_ratio, dk_dr,   &
+     &          band_s_evo)
 !
       use m_ludcmp_3band
       use center_sph_matrices
@@ -139,6 +140,8 @@ subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect,   &
       real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
       real(kind = kreal), intent(in) :: dt
       real(kind = kreal), intent(in) :: coef_advect
+      real(kind = kreal), intent(in) :: k_ratio(0:sph_rj%nidx_rj(1))
+      real(kind = kreal), intent(in) :: dk_dr(0:sph_rj%nidx_rj(1))
       character(len=kchara), intent(in) :: mat_name
 !
       type(band_matrices_type), intent(inout) :: band_s_evo
@@ -177,10 +180,18 @@ subroutine const_radial_mat_4_scalar_sph(mat_name, coef_advect,   &
       end if
 !
 !
-      call add_scalar_poisson_mat_sph                                   &
-     &   (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj,        &
-     &    g_sph_rj, sph_bc%kr_in, sph_bc%kr_out, r_coef(1),             &
-     &    r_2nd%fdm(1)%dmat, r_2nd%fdm(2)%dmat, band_s_evo%mat)
+      if(property%flag_val_diffuse) then
+        call add_scalar_r_diffuse_mat_sph                               &
+     &     (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj,      &
+     &      g_sph_rj, sph_bc%kr_in, sph_bc%kr_out,                      &
+     &      r_coef(1), k_ratio(1), dk_dr(1),                            &
+     &      r_2nd%fdm(1)%dmat, r_2nd%fdm(2)%dmat, band_s_evo%mat)
+      else
+        call add_scalar_poisson_mat_sph                                 &
+     &     (sph_rj%nidx_rj(1), sph_rj%nidx_rj(2), sph_rj%ar_1d_rj,      &
+     &      g_sph_rj, sph_bc%kr_in, sph_bc%kr_out, r_coef(1),           &
+     &      r_2nd%fdm(1)%dmat, r_2nd%fdm(2)%dmat, band_s_evo%mat)
+      end if
 !
       call sel_radial_mat_scalar_bc_sph(sph_rj, sph_bc, bcs_S,          &
      &    fdm2_center, g_sph_rj, r_coef, band_s_evo)

src/Fortran_libraries/MHD_src/sph_MHD/const_radial_mat_4_sph.f90

@@ -8,14 +8,15 @@
 !!
 !!@verbatim
 !!      subroutine const_radial_mat_sph_mhd(dt, MHD_prop, sph_MHD_bc,   &
-!!     &          sph, r_2nd, leg, sph_MHD_mat)
+!!     &          sph, r_2nd, leg, refs, sph_MHD_mat)
 !!      subroutine const_radial_mat_sph_snap(MHD_prop, sph_MHD_bc,      &
 !!     &          sph_rj, r_2nd, leg, sph_MHD_mat)
 !!        type(MHD_evolution_param), intent(in) :: MHD_prop
 !!        type(sph_MHD_boundary_data), intent(in) :: sph_MHD_bc
 !!        type(sph_grids), intent(in) :: sph
 !!        type(fdm_matrices), intent(in) :: r_2nd
 !!        type(legendre_4_sph_trans), intent(in) :: leg
+!!        type(radial_reference_field), intent(in) :: refs
 !!        type(MHD_radial_matrices), intent(inout) :: sph_MHD_mat
 !!@endverbatim
 !
@@ -33,7 +34,10 @@ module const_radial_mat_4_sph
       use t_spheric_rj_data
       use t_fdm_coefs
       use t_schmidt_poly_on_rtm
+      use t_diffusion_term_labels
       use t_physical_property
+      use t_phys_data
+      use t_radial_reference_field
       use t_boundary_data_sph_MHD
       use t_boundary_params_sph_MHD
       use t_radial_matrices_sph_MHD
@@ -66,13 +70,14 @@ module const_radial_mat_4_sph
 ! -----------------------------------------------------------------------
 !
       subroutine const_radial_mat_sph_mhd(dt, MHD_prop, sph_MHD_bc,     &
-     &          sph, r_2nd, leg, sph_MHD_mat)
+     &          sph, r_2nd, leg, refs, sph_MHD_mat)
 !
       type(MHD_evolution_param), intent(in) :: MHD_prop
       type(sph_MHD_boundary_data), intent(in) :: sph_MHD_bc
       type(sph_grids), intent(in) :: sph
       type(fdm_matrices), intent(in) :: r_2nd
       type(legendre_4_sph_trans), intent(in) :: leg
+      type(radial_reference_field), intent(in) :: refs
 !
       real(kind = kreal), intent(in) :: dt
 !
@@ -83,7 +88,9 @@ subroutine const_radial_mat_sph_mhd(dt, MHD_prop, sph_MHD_bc,     &
       id_file = 50 + my_rank
       call const_radial_matrices_sph                                    &
      &   (id_file, dt, sph%sph_params, sph%sph_rj, r_2nd,               &
-     &    MHD_prop, sph_MHD_bc, leg%g_sph_rj, sph_MHD_mat)
+     &    MHD_prop, sph_MHD_bc, leg%g_sph_rj,                           &
+     &    refs%iref_diffusivity, refs%iref_grad_diffusivity,            &
+     &    refs%ref_field, sph_MHD_mat)
 !
       call const_radial_mat_sph_w_center                                &
      &   (dt, sph%sph_rj, MHD_prop, sph_MHD_bc, sph_MHD_mat)
@@ -136,8 +143,10 @@ end subroutine const_radial_mat_sph_snap
 ! -----------------------------------------------------------------------
 !
       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)
+     &         (id_file, dt, sph_params, sph_rj,                        &
+     &          r_2nd, MHD_prop, sph_MHD_bc, g_sph_rj,                  &
+     &          iref_diffusivity, iref_grad_diffusivity, ref_field,     &
+     &          sph_MHD_mat)
 !
       use const_r_mat_4_scalar_sph
       use const_r_mat_4_vector_sph
@@ -149,6 +158,10 @@ subroutine const_radial_matrices_sph                              &
       type(fdm_matrices), intent(in) :: r_2nd
       type(MHD_evolution_param), intent(in) :: MHD_prop
       type(sph_MHD_boundary_data), intent(in) :: sph_MHD_bc
+!
+      type(diffusivity_adress), intent(in) :: iref_diffusivity
+      type(diffusivity_adress), intent(in) :: iref_grad_diffusivity
+      type(phys_data), intent(in) :: ref_field
 !
       real(kind = kreal), intent(in) :: dt
       real(kind = kreal), intent(in) :: g_sph_rj(sph_rj%nidx_rj(2),13)
@@ -223,6 +236,8 @@ subroutine const_radial_matrices_sph                              &
      &      sph_params, sph_rj, r_2nd, MHD_prop%ht_prop,                &
      &      sph_MHD_bc%sph_bc_T, sph_MHD_bc%bcs_T,                      &
      &      sph_MHD_bc%fdm2_center, g_sph_rj,                           &
+     &      ref_field%d_fld(1,iref_diffusivity%i_T_diffusivity),        &
+     &      ref_field%d_fld(1,iref_grad_diffusivity%i_T_diffusivity),   &
      &      sph_MHD_mat%band_temp_evo)
         if(i_debug .eq. iflag_full_msg) then
           call check_radial_band_mat(id_file, sph_rj,                   &
@@ -236,6 +251,8 @@ subroutine const_radial_matrices_sph                              &
      &      sph_params, sph_rj, r_2nd, MHD_prop%cp_prop,                &
      &      sph_MHD_bc%sph_bc_C, sph_MHD_bc%bcs_C,                      &
      &      sph_MHD_bc%fdm2_center, g_sph_rj,                           &
+     &      ref_field%d_fld(1,iref_diffusivity%i_T_diffusivity),        &
+     &      ref_field%d_fld(1,iref_grad_diffusivity%i_T_diffusivity),   &
      &      sph_MHD_mat%band_comp_evo)
         if(i_debug .eq. iflag_full_msg) then
           call check_radial_band_mat(id_file, sph_rj,                   &

src/programs/SPH_MHD/SPH_analyzer_MHD.f90

@@ -144,7 +144,7 @@ subroutine SPH_initialize_MHD(MHD_files, SPH_model, FEM_dat,      &
       call const_radial_mat_sph_mhd                                     &
      &   (MHD_step%time_d%dt, SPH_model%MHD_prop, SPH_model%sph_MHD_bc, &
      &    SPH_MHD%sph, SPH_WK%r_2nd, SPH_WK%trans_p%leg,                &
-     &    SPH_WK%MHD_mats)
+     &    SPH_model%refs, SPH_WK%MHD_mats)
 !*
 !* obtain linear terms for starting
 !*