Merge branch 'main' of github.com:cpinte/mcfost

src/MRW.f90

@@ -8,6 +8,9 @@ module MRW
   integer, parameter :: n = 10000
   real(dp), dimension(n), save :: zeta, y_MRW
 
+  real, parameter :: gamma_MRW = 2.0
+  real(dp), parameter :: cst_ct = 3 / pi**2
+
 contains
 
  subroutine initialize_cumulative_zeta()
@@ -68,37 +71,44 @@ end function sample_zeta
 
  !----------------------------------------------
 
- subroutine make_MRW_step(id,icell,lambda, x,y,z,R0, E)
+ subroutine make_MRW_step(id,icell, x,y,z,E, d, rec_Planck_opacity)
 
    use naleat, only : random_isotropic_direction
 
    integer, intent(in) :: id, icell
-   integer, intent(inout) :: lambda
    real(kind=dp), intent(inout) :: x,y,z
-   real(kind=dp), intent(in) :: R0, E
+   real(kind=dp), intent(in) :: E, d, rec_Planck_opacity
 
-   real(dp) :: u,v,w, ct, diff_coeff
+   real(dp) :: u,v,w, ct
 
    ! Place photon randomly on sphere of radius R0 around current position
    call random_isotropic_direction(id, u,v,w)
-   x = x + u*R0
-   y = y + v*R0
-   z = z + w*R0
+   x = x + u*d
+   y = y + v*d
+   z = z + w*d
 
    ! Select new direction : isotropic or emitting sphere ?
 
    ! sample P0
    y = sample_zeta(id)
 
    ! Solve Eq. 8 of Min et al 2009
-   ct = -log(y) / diff_coeff * (R0/pi)**2.
+   ! D = 1/(3*rec_Planck_opacity)
+   ! cst_ct = 3/pi**3
+   ct = -log(y) * cst_ct * rec_Planck_opacity * d**2.
+
 
+   ! Steps that needs to be done
    ! Deposit energy using Planck mean opacity
 
    ! Update temperature and diffusion coefficient ?
 
+
+   ! Only at end of MRW, to switch to MC
    ! Select new wavelength
 
+   ! Select new photon direction
+
 
    return
 

src/Voronoi.f90

@@ -993,10 +993,10 @@ end subroutine cross_Voronoi_cell
 
   !----------------------------------------
 
-  real(dp) function distance_to_closest_wall_Voronoi(id,icell,x,y,z) result(s)
+  real(dp) function distance_to_closest_wall_Voronoi(icell,x,y,z) result(s)
 
 
-    integer, intent(in) :: id, icell
+    integer, intent(in) :: icell
     real(kind=dp), intent(in) :: x,y,z
 
     ! n = normale a la face, p = point sur la face, r = position du photon, k = direction de vol

src/cylindrical_grid.f90

@@ -10,7 +10,8 @@ module cylindrical_grid
        move_to_grid_cyl, define_cylindrical_grid, build_cylindrical_cell_mapping,  cell_map, cell_map_i, cell_map_j,&
        cell_map_k, lexit_cell, r_lim, r_lim_2, r_lim_3, delta_z, r_grid, z_grid, phi_grid, tab_region, &
        z_lim, w_lim, theta_lim, tan_theta_lim, tan_phi_lim, cos_phi_lim, sin_phi_lim, volume, l_dark_zone, zmax, &
-       delta_cell_dark_zone, ri_in_dark_zone, ri_out_dark_zone, zj_sup_dark_zone, zj_inf_dark_zone, l_is_dark_zone
+       delta_cell_dark_zone, ri_in_dark_zone, ri_out_dark_zone, zj_sup_dark_zone, zj_inf_dark_zone, l_is_dark_zone, &
+       distance_to_closest_wall_cyl
 
   real(kind=dp), parameter, public :: prec_grille=1.0e-14_dp
 
@@ -1177,9 +1178,9 @@ end subroutine cross_cylindrical_cell
 
   !***********************************************************
 
-  real(dp) function distance_to_closest_wall_cyl(id,icell,x,y,z) result(s)
+  real(dp) function distance_to_closest_wall_cyl(icell,x,y,z) result(s)
 
-    integer, intent(in) :: id, icell
+    integer, intent(in) :: icell
     real(kind=dp), intent(in) :: x,y,z
 
     real(dp) :: r,s1,s2,s3,s4,s5,s6,z0
@@ -1196,8 +1197,8 @@ real(dp) function distance_to_closest_wall_cyl(id,icell,x,y,z) result(s)
     ! z walls
     z0 = abs(z)
     zj0 = abs(zj0)
-    s3 = z_lim(ri0,abs(zj0)+1) - z
-    s4 = z - z_lim(ri0,abs(zj0))
+    s3 = z_lim(ri0,abs(zj0)+1) - z0
+    s4 = z0 - z_lim(ri0,abs(zj0))
 
     if (l3D) then
        ! phi walls

src/diffusion.f90

@@ -628,7 +628,7 @@ end subroutine diffusion_approx_nLTE_nRE
 
 !************************************************************
 
-subroutine diffusion_opacity(icell, Planck_opacity,rec_Planck_opacity)
+subroutine compute_Planck_opacities(icell, Planck_opacity,rec_Planck_opacity)
   ! Compute current Planck reciprocal mean opacity for all cells
   ! (note : diffusion coefficient needs to be defined with Rosseland opacity
   ! in B&W mode)
@@ -648,16 +648,17 @@ subroutine diffusion_opacity(icell, Planck_opacity,rec_Planck_opacity)
   real(dp), intent(out) :: Planck_opacity,rec_Planck_opacity ! cm2/g (ie per gram of gas)
 
   integer :: lambda
-  real(dp) :: somme, somme2, cst, cst_wl, B, dB_dT, coeff_exp, wl, delta_wl, norm, Temp
+  real(dp) :: somme, somme2, cst, cst_wl, B, dB_dT, coeff_exp, wl, delta_wl, norm, T
 
   integer, pointer :: p_icell
   integer, target :: icell0
 
   icell0 = icell
 
-  temp = Tdust(icell)
-  if ((temp > 1) .and. (Voronoi(icell)%original_id > 0)) then
-
+  T = Tdust(icell) ! WARNING : this is 0 until Temp_finale
+  T = 20
+!  if ((T > 1) .and. (Voronoi(icell)%original_id > 0)) then
+  if (T > 1) then
      if (lvariable_dust) then
         p_icell => icell0
      else
@@ -667,7 +668,7 @@ subroutine diffusion_opacity(icell, Planck_opacity,rec_Planck_opacity)
      somme  = 0.0_dp
      somme2 = 0.0_dp
      norm = 0.0_dp
-     cst    = cst_th/temp
+     cst    = cst_th/T
      do lambda = 1,n_lambda
         ! longueur d'onde en metre
         wl       = tab_lambda(lambda)*1.e-6
@@ -681,18 +682,18 @@ subroutine diffusion_opacity(icell, Planck_opacity,rec_Planck_opacity)
            B = 0.0_dp
            !dB_dT = 0.0_dp
         endif
-        somme  = somme  + B/(kappa(p_icell,lambda) * kappa_factor(icell0))*delta_wl
-        somme2  = somme2  + B * (kappa(p_icell,lambda) * kappa_factor(icell0))*delta_wl
+        somme  = somme  + B/kappa(p_icell,lambda)
+        somme2  = somme2  + B * kappa(p_icell,lambda)
         norm = norm + B*delta_wl
      enddo
-     rec_Planck_opacity = norm/somme
-     Planck_opacity = somme2/norm
+     rec_Planck_opacity = norm/somme * kappa_factor(icell0)
+     Planck_opacity = somme2/norm * kappa_factor(icell0)
   else
-     rec_Planck_opacity = 0.
-     Planck_opacity = 0.
+     rec_Planck_opacity = 1e-30
+     Planck_opacity = 1e-30
   endif
 
-end subroutine diffusion_opacity
+end subroutine compute_Planck_opacities
 
 !************************************************************
 

src/dust_transfer.f90

@@ -43,7 +43,7 @@ subroutine transfert_poussiere()
 ! 20/04/2023
 
   use thermal_emission, only : frac_E_stars, frac_E_disk
-  use MRW, only : make_MRW_step, initialize_cumulative_zeta
+  use MRW, only : initialize_cumulative_zeta
   use utils
   implicit none
 
@@ -966,6 +966,8 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
   ! - separer 1ere diffusion et reste
   ! - lom supprime !
 
+  use MRW, only : make_MRW_step, gamma_MRW
+
   integer, intent(in) :: id
   integer, intent(inout) :: lambda, p_lambda
   integer, target, intent(inout) :: icell
@@ -976,9 +978,9 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
   logical, intent(out) :: flag_scatt
 
   real(kind=dp), dimension(4,4) :: M
-  real(kind=dp) :: u1,v1,w1, phi, cospsi, w02, srw02, argmt
+  real(kind=dp) :: u1,v1,w1, phi, cospsi, w02, srw02, argmt, Planck_opacity, rec_Planck_opacity, d, diff_coeff
   integer :: taille_grain, itheta
-  integer :: n_iterations
+  integer :: n_iteractions_in_cell, icell_old
   integer, pointer :: p_icell
   real :: rand, rand2, tau, dvol
 
@@ -1002,7 +1004,7 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
   ! Boucle sur les interactions du paquets:
   ! - on avance le paquet
   ! - on le fait interagir avec la poussiere si besoin
-  n_iterations = 0
+  n_iteractions_in_cell = 0
   infinie : do
 
      ! Longueur de vol
@@ -1020,26 +1022,34 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
      !   if (.not.flag_star) Stokes=0.
      !endif
 
-!     d = 0.
-!     reciprocal_Plank_kappa = 1.
-!     if ((n_iterations > 5) .and. Tdust(icell) > 1) then
-!        d = distance_to_closest_wall()
-!        call diffusion_opacity()
-!     endif
-!
-!     do while  (d > gamma_MRW * reciprocal_Plank_kappa)
-!        call make_MRW_step(id,icell, x,y,z,d, Stokes(1))
-!        d = distance_to_closest_wall()
-!        ! todo : do we update the rec_plack_kappa ???
-!     enddo
-
- !    icell_old = icell
+     ! Do we need to perform a MRW ?
+     !if ((n_iteractions_in_cell > 5)) then
+     !   d = distance_to_closest_wall(icell,x,y,z)
+     !   call compute_Planck_opacities(icell, Planck_opacity,rec_Planck_opacity)
+     !
+     !   !write(*,*) icell, d, rec_Planck_opacity,  d * rec_Planck_opacity, gamma_MRW
+     !   !read(*,*)
+     !
+     !   do while  (d * rec_Planck_opacity > gamma_MRW )
+     !      write(*,*) "MRW"
+     !      call make_MRW_step(id,icell, x,y,z,Stokes(1), d, rec_Planck_opacity)
+     !      d = distance_to_closest_wall(icell,x,y,z)
+     !      ! todo : do we update the rec_plack_opacity ???
+     !      !call compute_Planck_opacities(icell, Planck_opacity,rec_Planck_opacity)
+     !   enddo
+     !
+     !   ! MRW is finished, we choose wl and direction
+     !endif
+
+     ! now that MRW is done, we perform the normal propagation
+     icell_old = icell
      call physical_length(id,lambda,p_lambda,Stokes,icell,x,y,z,u,v,w,flag_star,flag_direct_star,tau,dvol,flag_sortie,lpacket_alive)
- !    if (icell == icell_old) then
- !       n_iterations = n_iterations + 1
- !    else
- !       n_iterations = 0
- !    endif
+     if (icell == icell_old) then
+        n_iteractions_in_cell = n_iteractions_in_cell + 1
+        !write(*,*) "icell=", icell, n_iteractions_in_cell
+     else
+        n_iteractions_in_cell = 0
+     endif
 
      if (flag_sortie) return ! Vie du photon terminee
 
@@ -1074,7 +1084,7 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
      endif ! lmono
 
 
-     if (rand < tab_albedo_pos(p_icell,lambda)) then ! Diffusion
+     if (rand < tab_albedo_pos(p_icell,lambda)) then ! Scatttering event
         flag_scatt=.true.
         flag_direct_star = .false.
 
@@ -1143,7 +1153,7 @@ subroutine propagate_packet(id,lambda,p_lambda,icell,x,y,z,u,v,w,stokes,flag_sta
         ! Mise a jour direction de vol
         u = u1 ; v = v1 ; w = w1
 
-     else ! Absorption
+     else ! Absorption + eventual re-emission
 
         if ((.not.lmono).and.lnRE) then
            ! fraction d'energie absorbee par les grains hors equilibre

src/grid.f90

@@ -19,7 +19,7 @@ module grid
   procedure(indice_cellule_cyl), pointer :: indice_cellule => null()
   procedure(test_exit_grid_cyl), pointer :: test_exit_grid => null()
   procedure(define_cylindrical_grid), pointer :: define_grid => null()
-  procedure(distance_to_closest_wall_Voronoi), pointer :: distance_to_closest_wall => null()
+  procedure(distance_to_closest_wall_cyl), pointer :: distance_to_closest_wall => null()
 
   real(kind=dp) :: v_char, B_char
   logical :: lcalc_ne, lmagnetized
@@ -343,6 +343,7 @@ subroutine setup_grid()
         indice_cellule => indice_cellule_cyl
         test_exit_grid => test_exit_grid_cyl
         define_grid => define_cylindrical_grid
+        distance_to_closest_wall => distance_to_closest_wall_cyl
      else if (grid_type == 2) then
         lcylindrical = .false.
         lspherical = .true.
@@ -352,6 +353,7 @@ subroutine setup_grid()
         indice_cellule => indice_cellule_sph
         test_exit_grid => test_exit_grid_sph
         define_grid => define_cylindrical_grid ! same routine at the moment
+        distance_to_closest_wall => distance_to_closest_wall_sph
      else
         call error("Unknown grid type")
      endif

src/spherical_grid.f90

@@ -10,7 +10,7 @@ module spherical_grid
   implicit none
 
   public :: cross_spherical_cell, pos_em_cellule_sph, indice_cellule_sph, test_exit_grid_sph, &
-       move_to_grid_sph
+       move_to_grid_sph, distance_to_closest_wall_sph
 
 
   private
@@ -448,9 +448,9 @@ end subroutine cross_spherical_cell
 
   !***********************************************************
 
-  real(dp) function distance_to_closest_wall_sph(id,icell,x,y,z) result(s)
+  real(dp) function distance_to_closest_wall_sph(icell,x,y,z) result(s)
 
-    integer, intent(in) :: id, icell
+    integer, intent(in) :: icell
     real(kind=dp), intent(in) :: x,y,z
 
     real(dp) :: r,s1,s2,s3,s4,s5,s6,r2_cyl,rcyl,z0