Skip to content

hotfix

Leo Michel-Dansac requested to merge hotfix into master

This merge request fixes the problem with negative travel time found previously and the bugs with precision. The code seems to be robust now and does not produce any status==3 photons anymore. The following modifications/corrections have been made:

  1. module_photon.f90: ppos => pcell: use cell center pos now (instead of photon's pos as we wrongly did before)
pcell = cell_corner + 0.5d0*cell_size
cell_fully_in_domain = domain_contains_cell(pcell,cell_size,domaine_calcul)
  1. module_photon.f90: now compute distance to escape of domain border along k. Function domain_distance_to_border_along_k routine added in module_domain.f90.

  2. module_photon.f90 : correct computation of tau in the case of a cell not fully contained in the computational domain. In this case, photon may escape domain before escaping cell => add a test and take the min between distance to cell border and distance to domain border:

! compute distance of photon to border of cell along propagation direction
distance_to_border           = path(ppos_cell,p%k)                   ! in cell units
distance_to_border_cm        = distance_to_border * cell_size_cm     ! cm
distance_to_border_box_units = distance_to_border * cell_size        ! in box units
! if cell not fully in domain, modify distance_to_border to "distance_to_domain_border" if relevant
OutOfDomainBeforeCell = .False.
if (.not.(cell_fully_in_domain)) then
        dborder    = domain_distance_to_border_along_k(ppos,p%k,domaine_calcul)  ! in box units
        dborder_cm = dborder * box_size_cm                                       ! from box units to cm
        ! compare distance to cell border and distance to domain border and take the min
        if (dborder_cm < distance_to_border_cm) then
            OutOfDomainBeforeCell        = .True.
            distance_to_border_cm        = dborder_cm
            distance_to_border_box_units = dborder
            distance_to_border           = distance_to_border_cm / cell_size_cm
        end if
end if
  1. module_domain.f90: domain_distance_to_border debugged (for cube and slab geometries), domain_contains_cell debugged (for all geometries), and make function get_my_new_domain working for any number of overlapping domains.

  2. all (?) real parameters in double precision now (just to be sure...), except the one in the nearest call in ran3.

  3. module_utils.f90: explicitly force normalization of k vectors at numerical precision

  4. module_photon.f90: implement an iterative epsilon pushing method to be sure to move photon outside cells at numerical precision.

  5. module_photon.f90: subroutine propagate has been slightly reshaped due to previous modifications. We now have the following metacode:

[...]
compute scatter_flag
if (scatter_flag == 0) then   ! next scattering event will not occur in the cell or in the domain
        move photon out of cell or domain
        ...
        if (OutOfDomainBeforeCell) then ! photon exits computational domain and is done 
            ...
            exit photon_propagation
        else
            ! photon exits current cell -> find into which new cell it goes
            call whereIsPhotonGoing()
            ...
            exit propag_in_cell
        end if
else
        ! Next event happens inside this cell and in the domain.
        scattering
        ...
end if 
Edited by Leo Michel-Dansac

Merge request reports