Morph Drain Algorithm

[kbruning]

At pc = 1.81122 for a [140,30,30] ink bottle, there is a difference in the morph drain solution between my original implementation and the PMMoTo implementation. The difference occurs at Step 3g with the morphology.morph call. I would like to look at this step in more detail but it is difficult for me to read and see exactly what is going on. As a starting point, could you have this function save the new EDT as a csv and add comments describing the steps taken?

[tmweigand]

@kbruning When starting to add that to morpholgy.py, I found a bug with the reservoir condition. Please re-pull and see if I fixed your issue. If not, let me know and I can add the csv write but we will need to chat about how you want that csv data because I handle morph operations in a slightly weird way.

[kbruning]

Unfortunately made the difference larger. I'll email you a zoom link.

[kbruning]

I have narrowed down the difference to the edt calculation within morphadd, see the plot below (everything above ~x=70 is being ignored). My serial implementation is using distance_transform_edt from scipy.ndimage. If I use edt3d outside of PMMoTo, the error between the plots below goes to e-6 as shown in the second set of plots below. We are matching at step 3f, and the error between our approaches is e-8 for the first edt calculation in calcdrainage. I'm not sure where to go from here.

[tmweigand]

Okay. For clarity, the top three plots are comparing pmmoto and the bottom three are comparing edt? I think it must be an issue with the reservoirs and halo communication. I will look into it.

[kbruning]

The top plots are the data output straight from PMMoTo vs data output from my own implementation. The bottom plots are comparing the EDT methods in a totally new python script. Thanks!

[tmweigand]

@kbruning I just pushed a commit to timMA that fixes this issue.

[kbruning]

I pulled and confirmed that I am on the correct commit, but am getting the same issue. You are getting a new EDT and a different result?

[kbruning]

We discovered a few days ago that this difference is due to what is happening in the nonwetting phase reservoir that was not being visualized above. The difference between implementations actually starts at the beginning of calcDrainage with the first calcEDT function. The PMMoTo calcEDT gives the result below. It appears that there is a partial solid wall at the first voxel which should be removed, and I'm not sure what is happening at the reservoir/domain boundary. Step 2 seems to be working just fine, but step 3c is being impacted by the disconnect at the reservoir/domain boundary.

[tmweigand]

@kbruning I just started looking into this. I increased the reservoir size to 20 and the poreSpaceDist map I am getting seems correct. How did you generate that image? Here is what I have.

[kbruning]

You are correct, the error was mine, thanks! I was not reading the csv reservoir data correctly. I now have an identical match including the reservoirs through Step 3f, these are 10 voxel top and bottom reservoirs:

I believe I have the difference in the final result figured out. We have the same input into the second EDT, but it looks like the PMMoTo version might not be considering any positive values of the indicator function within the reservoir:

Also if you get a moment, could you show me how to plot the reservoirs in paraview? I haven't been able to figure it out.

[tmweigand]

To view the reservoirs in paraview, you have to open up the output for each processor. For example, open grid/gridProc.0.vtr. The package I am using to generate the files forces the complete domain to start at (0,0,0) so the reservoirs for this problem are cut off. I raised Issue #32 but haven't been able to fix it yet.

[kbruning]

I found this issue and how to fix it. The checks in Step 3f should be reduced to only the second one, if nwCheck..., and the other two should be removed. I should have taken your recommendation to look closely at this step earlier!

This gets us identical answers for every ink bottle drainage case except for the first few pc points from the analytic solution. To fix this part, Step 3 eqDist.checkPoints should also consider points within the reservoirs, but the other eqDist.checkPoints calls within calcDrainage (nwCheck and wCheck) should remain as is and not consider points within the reservoirs.

[tmweigand]

Great! I will make those changes as soon as possible. I am still working on the medial extraction with some urgency.

[kbruning]

I have a further correction to this issue. The above works for the ink bottle, but to work for other media as well it should actually be changed to what I have below, thanks!

Steb 3c and 3d - Already checked at Step 3 so Collect Sets with ID = 1

indSets,indSetCount = sets.collectSets(ind,1,mP.inlet[mP.nwID],mP.outlet[mP.nwID],mP.loopInfo[mP.nwID],mP.subDomain)

ind = eqDist.getInletConnectedNodes(indSets,1)

ind2 = eqDist.getInletConnectedNodes(indSets,1)

Step 3f -- Unsure about these checks!

if wCheck and nwCheck:

ind = np.where( (ind == 1) & (nwGrid == 1) & (wGrid == 0),1,0).astype(np.uint8)

if nwCheck: ind = np.where( (ind2 != 1) & (nwGrid != 1),0,ind).astype(np.uint8)

elif wCheck:

ind = np.where( (ind == 1) & (wGrid == 1),1,0).astype(np.uint8)

[tmweigand]

@kbruning are the edits above the only edits necessary or do you still need checkPoints to be able to include/exclude reservoirs?

[kbruning]

Yes, also the checkPoints piece on including/excluding reservoirs.

[tmweigand]

I just pushed some code to main that allows the option to examine the reservoir with checkPoints. You set includeInlet = True to include the reservoirs. I also modified those if statements so that we only have the one. I havent tested this code so please let me know if it is not working.

Also, I did a major refactor to the code to include the new medial extraction algorithm. I think all of your code should still be fine but I did modify some of communications. It would be great if we could develop a unit test for you algorithms and my code so we know if things break.

[kbruning]

Thanks for the modifications. I have been unable to switch to includeInlet = True. I get the error "AttributeError: 'equilibriumDistribution' object has no attribute 'ownNodesIndex'" from step 3 when changing to "def checkPoints(self,grid,ID,includeInlet=True):" and leaving everything else as is in the commit other than the KDTree import in distance.pyx.

Once we get this drainage function working again I will start on that unit test.

[tmweigand]

Sorry about that. That is fixed now.

[kbruning]

Thanks! I was able to use this to get the ink bottle result I was looking for. Step 3 includes the inlet, steps 3a and 3b (nwCheck and wCheck) do not.

[tmweigand]

Great!

[kbruning]

The KB_morphtools branch can match my drainage, imbibition, and scanning curve implementations.