hotfix
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:
- 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)
-
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.
-
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
-
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.
-
all (?) real parameters in double precision now (just to be sure...), except the one in the nearest call in ran3.
-
module_utils.f90: explicitly force normalization of k vectors at numerical precision
-
module_photon.f90: implement an iterative epsilon pushing method to be sure to move photon outside cells at numerical precision.
-
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