Unable to run MMC with high photon count

[Viet Tran Ba]

Dear Dr. Fang,

thank you for your work on MMC and iso2mesh, and making them open-source. 

I am trying to setup a simple mesh using pyiso2mesh and using MMC to simulate the tissue model. However, when I try to run MMC with a high photon count (f. ex. 1e6), I always get the following error:

MMC ERROR(36):Invalid command queue in unit mmc_cl_host.c:627

I do not observe such errors with a lower photon count, so I am not sure what the cause of the bug is. Could you help me in finding out on what I did wrong?

---

I use an RTX 3090 in Ubuntu 22.06 and am using the latest MMC commit in the master branch (f86c824b2581fbb6defdb4b18b618918339fca90). I currently run the MMC command using the following command: 

path-to-git-project/mmc/bin/mmc -f path-to-mmc-input-folder/input.json -H 1000000 -O F -b 1 -F jnii -e 1e-3 -x 1

---

Unfortunately, I do not find the option to attach any files, so I will write the code directly into this mail.

Mesh creation and saving: 

To create a mesh, I use pyiso2mesh in Python:

import pathlib
import numpy as np
import iso2mesh as i2m

semantic_volume = np.concatenate([np.ones([32, 32, 2]), np.ones([32, 32, 16]) * 2], axis=2).astype(np.uint8)
nodes, elements, _ = i2m.v2m(semantic_volume, [], 10, 10, "cgalmesh")
elements, _, _ = i2m.meshreorient(nodes, elements)

i2m.plotmesh(nodes, elements)
i2m.plotmesh(nodes, elements[elements[:, -1] != 1])

# Saving the mesh nodes and elements

output_folder_path = pathlib.Path(...)
key = "volume"

node_positions = nodes[:, :3]
row_indices = np.arange(1, node_positions.shape[0] + 1).reshape(-1, 1)
node_data = np.hstack([row_indices, node_positions])
node_file_path = output_folder_path / f"node_{key}.dat"

with open(node_file_path, 'w') as f:
    # Note: Not sure whether the first element in the first line is 1 or the number of labels.
    f.write(f"1\t{node_data.shape[0]}\n")
    np.savetxt(f, node_data, fmt="%d\t%16.8e\t%16.8e\t%16.8e", delimiter='\t')

# Adds the index as first column
row_indices = np.arange(1, elements.shape[0] + 1).reshape(-1, 1)
element_data = np.hstack([row_indices, elements])
element_file_path = output_folder_path / f"elem_{key}.dat"
with open(element_file_path, 'w') as f:
    f.write(f"1\t{element_data.shape[0]}\n")
    np.savetxt(f, element_data, fmt="%d\t%d\t%d\t%d\t%d\t%d", delimiter='\t')

---

input.json

{
"Mesh": {
"MeshID": "volume",
"InitElem": 1
},
"Session": {
"Photons": 1000000.0,
"Seed": 0,
"ID": "volume"
},
"Forward": {
"T0": 0.0,
"T1": 1.0,
"Dt": 1.0
},
"Optode": {
"Source": {
"Type": "planar",
"Pos": [
0.0,
0.0,
0.25
],
"Dir": [
0.0,
0.0,
1.0,
"_NaN_"
],
"Param1": [
33.0,
0.0,
0.0,
0.0
],
"Param2": [
0.0,
33.0,
0.0,
0.0
]
},
"Detector": [
{
"Pos": [
17.0,
17.0,
0.0
],
"R": 32
}
]
}
}

---

prop_volume.dat

1 2
1   1.00000000e-18   1.00000000e-18   1.00000000e+00   1.00000000e+00
2   1.20000000e-02   1.20000000e+00   9.00000000e-01   1.32000000e+00

---

Let me know if you require further information. Thank you for your answer in advance.

Best regards,

Tran Ba, Viet

[Qianqian Fang]

hi Viet,

I see you are using a planar source type - a planar source launches photons at various locations, it is one of what we called a "widefield" sources. If the source-aperture is outside of the mesh domain, it requires an extra "mesh re-tessellation" step, see our 2016 paper

https://opg.optica.org/boe/fulltext.cfm?uri=boe-7-1-171&id=333806

you can also see how this is done in our built-in mmclab examples:

https://github.com/fangq/mmc/blob/master/mmclab/example/demo_sfdi_2layer.m#L55-L60

if you don't have this step, you can not set the "InitElem" value (the index of the tetrahedron that encloses all launched photon position) in the input file because the photon launch position is outside of the mesh. I see you set InitElem to 1, but if a photo is launched outside of element#1, then mmc could crash and cause memory errors.

to launch a widefield source located inside the mesh is a bit tricky.

if your mesh is generated so that there is a single tetrahedron element completely enclosing your source aperture, then, in theory, you should be able to set InitElem (or cfg.e0) to the index of this element, just like a pencil beam source. I used "in theory" here because this did not actually work, until last week after this fix

https://github.com/fangq/mmc/commit/e68a427e8e9f3d962eaa5b6e97a9a4cc91e34cf4

if a widefield source embedded inside a mesh covers more than a single element, it is still possible to model it - in this case, you should find all elements that are potential photon launch sites, and label their element type (cfg.elemprop) to -1. This is a special label indicating that these are elements to be tested for photon launching. When running the simulation, mmc sequentially test each -1 labeled elements and find the initial element ID dynamically. The only down-side of this approach is that mmc treats -1 labeled elements as air (use cfg.prop(1,:) as optical property), so it is quite limited.

Qianqian

--
You received this message because you are subscribed to the Google Groups "mmc-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mmc-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/mmc-users/2dfe4db7-e73f-432d-bc3f-220d0ed80fcdn%40googlegroups.com.

[Qianqian Fang]

by the way, as I announced recently, we also have released pmmc - the Python version of mmc - on pypi

https://pypi.org/project/pmmc/

you don't really need to save mesh/input to files and run mmc binary for this simulation. You can just define a cfg dict and call pmmc.run(cfg) in Python, see examples

https://github.com/fangq/mmc/tree/master/pmmc#how-to-use

To view this discussion visit https://groups.google.com/d/msgid/mmc-users/d304319f-6237-4b80-8610-61363c1c7b45%40neu.edu.

[Viet Tran Ba]

Dear Dr. Fang, 

thank you for your answer! Since I do indeed want the light source to be within the mesh which represents air, I replaced the air label of the tetrahedrons and the init elem to be -1 (which was originally 1).

After removing the entry of the air label in the prop file, MMC successfully runs.

Regarding pmmc, I prefer using MMC directly in case I want to take a look or try out some debugging at the source code.

Best regards, 

Viet 

[Viet Tran Ba]

Dear Dr. Fang,

unfortunately, I have to reopen my issue as MMC throws the same error again with a different mesh.

I did a little debugging and found out that it appears somewhere within the infinite loop in the function onephoton() in mmc_core.cl after a small amount of iterations. But, I can't tell which line exactly causes that issue.

Could you perhaps take another look at my problem? I feel like something about the mesh setup causes that issue.

---

Updated mesh code:

semantic_volume = np.concatenate([np.ones([32, 32, 2]), np.ones([32, 32, 4]) * 2, np.ones([32, 32, 14]) * 3], axis=2).astype(np.uint8)

nodes, elements, _ = i2m.v2m(semantic_volume, [], 2, 10, "cgalmesh")

elements, _, _ = i2m.meshreorient(nodes, elements)

elements[elements[:, -1] == 1, -1] = -1  # All air labels are replaced with -1

nodes_copy = np.copy(nodes)
nodes_copy[:, 2] = -nodes_copy[:, 2]  # Invert the z-Axis to show the tissue correctly.
i2m.plotmesh(nodes_copy, elements)
i2m.plotmesh(nodes_copy, elements[elements[:, -1] != -1])

Updated prop_volume.dat (with one more entry for label 3)

1 2

2   1.20000000e-02   1.20000000e+00   9.00000000e-01   1.32000000e+00

3   1.40000000e-01   1.20000000e+00   9.00000000e-01   1.40000000e+00

Updated input.json

"InitElem": -1  (previously: 1)

---

Thank you for your answer in advance.

Best regards,

Viet