Rod geometry not detecting photons

[Qudzu]

Dear Daniel,

Thank you for making this software available. It looks really great to me, and potentially super-useful for my projects.

While I have no problem with LSC (box) geometry, I have encountered an issue with rod (cylinder). I cannot trace photons exiting through its tips (surfaces "cap" and "base"), however, the third one ("hull") seems to work.
Could you give me an advice on how to solve this? I have tried to look into definitions of those surfaces, as well as tracing conditions, but I was not able to identify a solution.

All the best,
Konrad

[danieljfarrell]

Hi Konrad,

I’ll take a look. What are you interesting in simulating? Just a basic rod geometry with a luminescent dye inside?

Can you paste your python script here so I can run it and reproduce the same error?

[Qudzu]

Hi Daniel,

For now, we are interested in both cylindrical and planar LSCs - our target is to change the material (dye and matrix) parameters, to see the impact on efficiency.

At the moment, I use more-or-less a code based on the examples provided with the software.
I changed a bit a definition of the planar soruce.

import numpy as np
import pvtrace
from pvtrace.external import transformations as tf
from pvtrace import *

#Create light source from AM1.5 data, truncate

L=0.02
W=0.05

file = os.path.join(PVTDATA,'sources','AM1.5g-full.txt')
oriel = load_spectrum(file, xbins=np.arange(300,600))
source = PlanarSource(direction=(-1,0,0), spectrum=oriel, length=L, width=W) # Incident light AM1.5g spectrum
source.translate((0.02,-0.01,0))

#Load dye absorption and emission data, and create material
file = os.path.join(PVTDATA, 'dyes', 'lr300.abs.txt')
abs = load_spectrum(file)
file = os.path.join(PVTDATA, 'dyes', 'lr300.ems.txt')
ems = load_spectrum(file)
fluro_red = Material(absorption_data=abs, emission_data=ems, quantum_efficiency=1.00, refractive_index=1.5)

#Give the material a linear background absorption (pmma)
abs = Spectrum([0,1000], [0,0])
ems = Spectrum([0,1000], [0,0])
pmma = Material(absorption_data=abs, emission_data=ems, quantum_efficiency=0.0, refractive_index=1.5)

#Make the LSC

lsc = Rod(radius=0.01,length=0.05)
new_transform= tf.translation_matrix([0.00, 0.0, 0])
lsc.name = "LSC"
lsc.material = CompositeMaterial([pmma, fluro_red], refractive_index=1.5)
lsc.name = "LSC"

scene = Scene()
scene.add_object(lsc)

#Ask python that the directory of this script file is and use it as the location of the database file
pwd = os.getcwd()
dbfile = os.path.join(pwd, 'cylinder.sql') # <--- the name of the database file
if os.path.exists(dbfile):
    os.remove(dbfile)
trace = Tracer(scene=scene, source=source, seed=1, throws=5, database_file=dbfile, use_visualiser=True, show_log=False)
trace.show_lines = True
trace.show_path = True
import time
tic = time.clock()
trace.start()
toc = time.clock()

#Statistics
print("")
print("Run Time: ", toc - tic)
print("")

print("Technical details:")
generated = len(trace.database.uids_generated_photons())
killed = len(trace.database.killed())
thrown = generated - killed
print("\t Generated \t", generated)
print("\t Killed \t", killed)
print("\t Thrown \t", thrown)

print("Summary:")
print("\t Optical efficiency \t", (len(trace.database.uids_out_bound_on_surface('cap', luminescent=True)) + len(trace.database.uids_out_bound_on_surface('base', luminescent=True))) * 100 / thrown, "%")
print("\t Photon efficiency \t", (len(trace.database.uids_out_bound_on_surface('cap')) + len(trace.database.uids_out_bound_on_surface('base')) + len(trace.database.uids_out_bound_on_surface('hull'))) * 100 / thrown, "%")

print("Luminescent photons:")
edges = ['hull']
apertures = ['cap', 'base']

for surface in edges:
    print("\t", surface, "\t", len(trace.database.uids_out_bound_on_surface(surface, luminescent=True))/thrown * 100, "%")

for surface in apertures:
    print("\t", surface, "\t", len(trace.database.uids_out_bound_on_surface(surface, luminescent=True))/thrown * 100, "%")

print("\t", "Losses", len(trace.database.uids_nonradiative_losses())/thrown * 100, "%")

print("Solar photons (transmitted/reflected):")
for surface in edges:
    print("\t", surface, "\t", len(trace.database.uids_out_bound_on_surface(surface, solar=True))/thrown * 100, "%")

for surface in apertures:
    print("\t", surface, "\t", len(trace.database.uids_out_bound_on_surface(surface, solar=True))/thrown * 100, "%")

All the best,
Konrad

[danieljfarrell]

I have reproduced this error. I’m taking a closer look. This use to work, so I’m assumption It has been introduced with the python 3 changes.

By the way, pvtrace2 is almost ready to be released. It is much nicer! I will let you know when I have more news.

[Qudzu]

Dear Daniel,

Thank you! I will wait for your feedback, and if somehow I manage to find a solution, I will also let you know.

Please keep me updated about pvtrace2, I am excited to learn more.

All the best,
Konrad

[danieljfarrell]

Hello Konrad,

I managed to continue taking a look this morning. The problem seems to be here,

pvtrace/pvtrace/Geometry.py

Line 985 in 597f095

def surface_identifier(self, surface_point, assert_on_surface = True):

in the method that assigns an identity to the surface point. It only returns "hull" as you already said. I also opened the database file with https://sqlitebrowser.org/ to check that this is the case. Any ray that hits the cap or base is not getting identified correctly. I'll investigate further.

The weird thing I noticed is that the local coordinate system of the cylinder object should have the base at z=0 and cap at z=L (cylinder length). However, it seems rays are intersecting with the cap at L/2 and I cannot see any intersections with base. The surface identifier method labels points based on these values in the local coordinate system of the shape. If that has (somehow?) changed then this would cause the bug.

I don't know how urgent this is, but there are some ideas there if you wish to investigate yourself. I would start by using the SimpleSource and moving it to different locations and orientations, trace a single ray, and verify the intersection is correct and the surface has the correct identifier.

Yes, pvtrace2, is much nicer, it has unit tests so this type of bug would be caught. Also, it can trace meshes so can trace any closed shape. You could make a cylinder shape in Blender, import that and trace it. I am thinking about making a beta release soon.

[danieljfarrell]

Hi Konrad,

I have not returned to this issue yet. But I did push pvtrace2 to a branch. In principle, you could use pvtrace2 and a mesh of a cylinder to achieve what you want.

Pvtrace2 is not yet officially released because I’m still preparing things like documentation, pypi pages etc. But the code is frozen for the release now.

If are still interested I can help you get pvtrace2 working with your cylinder geometry.

[danieljfarrell]

Hello Konrad,

I added cylinder geometry to pvtrace 2.0 this weekend and below is a script that will get you started using the programme.

You should be able to install pvtrace (better in clean python environment) with the command pip install pvtrace then save this script below as a cyl.py and run it as python cyl.py.

Read the docs here which describe how to go about processing ray tracer statistical information to generate optical efficiencies for luminescence solar concentrators. You can apply the same concepts for your geometry.

""" This example simulates a thin and long cylinder which is sometime studied as a
    luminescent concentrator geometry.
"""
from pvtrace.geometry.cylinder import Cylinder
from pvtrace.geometry.sphere import Sphere
from pvtrace.scene.renderer import MeshcatRenderer
from pvtrace.scene.scene import Scene
from pvtrace.scene.node import Node
from pvtrace.light.light import Light
from pvtrace.trace.tracer import PhotonTracer
from pvtrace.material.material import Dielectric, Lumophore, Host
import numpy as np
import functools
import sys
import time

# Make a world node filled with air

world = Node(
    name="world (air)",
    geometry=Sphere(
        radius=10.0,
        material=Dielectric.air()
    )
)

# Make a very simple (and probably unphysical) absorption and emission.

def make_absorprtion_coefficient(x_range, wavelengths, absorption_coefficient, cutoff_range, min_alpha=0):
    wavelength1, wavelength2 = cutoff_range
    alpha = absorption_coefficient
    halfway = wavelength1 + 0.5 * (wavelength2 - wavelength1)
    x = [x_range[0], wavelength1, halfway, wavelength2, x_range[1]]
    y = [alpha, alpha, 0.5 * alpha, min_alpha, min_alpha]
    abs_coeff = np.interp(wavelengths, x, y)
    return abs_coeff

def make_emission_spectrum(x_range, wavelengths, cutoff_range, min_ems=0):
    wavelength1, wavelength2 = cutoff_range
    halfway = wavelength1 + 0.5 * (wavelength2 - wavelength1)
    x = [x_range[0], wavelength1, halfway, wavelength2, x_range[1]]
    y = [min_ems, min_ems, 1.0, min_ems, min_ems]
    abs_coeff = np.interp(wavelengths, x, y)
    return abs_coeff

# Use the functions above to make the spectrum needed by the Lumophore. At this point
# you may want to customise this script and import your own experimental data.

x_range = (300, 1000)
wavelength = np.linspace(*x_range)
abs_coef = make_absorprtion_coefficient(x_range, wavelength, 1.0, (700, 800))
ems_spec = make_emission_spectrum(x_range, wavelength, (600, 700))
lumophore = Lumophore(
    np.column_stack((wavelength, abs_coef)),  # abs. coef. spectrum
    np.column_stack((wavelength, ems_spec)),  # emission spectrum
    1.0  # quantum yield
)

# Make a host material which has a refractive index of 1.5 and also contains the 
# lumophore we just created.

material = Host(
    np.column_stack( # refractive index spectrum
        (wavelength,
         np.ones(wavelength.size) * 1.5)
    ), 
    [lumophore],  # list of lumophores, reuse the one we already have.
)

# Make the cylinder node with length 5cm and radius 0.02cm and give the material we 
# just made.

cylinder = Node(
    name="cylinder (glass)",
    geometry=Cylinder(
        length=5.0,
        radius=0.02,
        material=material
    ),
    parent=world
)

# Make a light source. This is a laser emitting light along the whole length of the
# cylinder. We need to translate and rotate the light source to get it to fire along
# the axis. We use the position delegate to generate photons along the same length
# as the cylinder.

light = Node(
    name="light (555nm laser)",
    light=Light(position_delegate=lambda : (np.random.uniform(-2.5, 2.5), 0.0, 0.0)),
    parent=world
)
cylinder.rotate(np.radians(90), (0.0, 1.0, 0.0))
light.translate((0.0, 0.0, -1.0))

# Setup renderer and the scene for tracing
rend = MeshcatRenderer()
scene = Scene(world)
tracer = PhotonTracer(scene)
rend.render(scene)

# Trace 100 photons and visualise
for ray in light.emit(100):
    path = tracer.follow(ray)
    print(path)
    rend.add_ray_path(path)

# You can kill the simulation any time by pressing ctrl-c.
# Do something with the data. Read the pvtrace documentation on how to process data
# for luminescence solar concentrators.
while True:
    try:
        time.sleep(.3)
    except KeyboardInterrupt:
        print("Bye")
        sys.exit()

[danieljfarrell]

Going to close this now. Anybody wanting to do cylinder/rod LSC simulations should use pvtrace 2