PMEAL / OpenPNM

A Python package for performing pore network modeling of porous media
http://openpnm.org
MIT License
442 stars 175 forks source link

Visualize Throats in Paraview #507

Closed jgostick closed 4 years ago

jgostick commented 8 years ago

This would be very useful to help ensure things are looking right, but at the moment my Paraview simulations only show lines between pores.

I recall that @jhinebau once sent me a complex recipe for putting throat cylinders in paraview (using a series of paraview commands), but I wonder if there is something we can do to the output VTP files to include exactly what paraview needs to simply add oriented gliphs?

GanttDue: 2018-10-25

jgostick commented 8 years ago

The example on the website already explains how to do this, but it is unreasonably complicated. There MUST be a way to put more information into the VTP files so that throat gliph information is present.

jgostick commented 7 years ago

I really want this! The current process is insane! We need a VTK expert to step up and help here.

TomTranter commented 6 years ago

I don't think there is a better way to do this with vtk and paraview. You could output separate sets of points data for pores and throats and then create separate glyphs but you can only scale by one factor so radius and length have to be set in two steps for the throats (which I haven't managed to figure out how to do properly - probably another set of convoluted filters). Or you could output polydata from OpenPNM leaving the user to control things like number of vertices in each glyph and making the output files much larger. Neither option is better than what we already have. With the GUI we could do all rendering in OpenPNM the way we want it. Paraview works much better with grids of data and images than with polygonal datasets imho

jgostick commented 6 years ago

I guess I'm thinking that the user should just be able to apply a tube filter to get tubes without the 5 step process paraview currently requires. We have throat end points, and throat diameter, so sure there is a shortcut to create the throats. We don't have any vector data in our files, maybe this is the key? If we had vectors we could just put tubes on them?

TomTranter commented 6 years ago

I can use pyevtk to output vector data but the cylinder glyph still needs a two step process. Not sure about tube

On 30 Aug 2017 15:56, "jgostick" notifications@github.com wrote:

I guess I'm thinking that the user should just be able to apply a tube filter to get tubes without the 5 step process paraview currently requires. We have throat end points, and throat diameter, so sure there is a shortcut to create the throats. We don't have any vector data in our files, maybe this is the key? If we had vectors we could just put tubes on them?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/PMEAL/OpenPNM/issues/507#issuecomment-326101102, or mute the thread https://github.com/notifications/unsubscribe-auth/AGvcRRzUgCUzpIKHtPLFzFEGwTKzNjjpks5sdb5RgaJpZM4HYP4d .

jgostick commented 6 years ago

Now that we have throat end points as a normal properties, perhaps we can revisit this and see if it'd doable?

xu-kai-xu commented 5 years ago

this is my problem, now i can export the python file to vtp format, by using this command

import openpnm as op pn = op.network.Cubic(shape=[10, 10, 10], spacing=0.0001) geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts) Hg = op.phases.Mercury(network=pn) phys = op.physics.Standard(network=pn, phase=Hg, geometry=geo) mip = op.algorithms.Porosimetry(network=pn) mip.setup(phase=Hg) mip.set_inlets(pores=pn.pores(['left', 'right', 'top', 'bottom'])) mip.run()

生成vtp格式的文件的命令

op.io.VTK.save(pn, phases=phys, filename='first_examples', delim='|')

did i use the final line code rightly? and i can show pores or throats seperately, but can not show them in the same view window. the first two is what i made.
p1 p2

jgostick commented 5 years ago

One of the "quirks" of openpnm is that you need to return the results of an algorithm back to a phase object. The idea is that you don't necessarily want your phase object to be updated with all kinds of data. Anyway, what this means is that you need to do something like:

Hg['pore.invasion_pressure'] = mip['pore.invasion_pressure']
Hg['throat.invasion_pressure'] = mip['throat.invasion_pressure']

Then when you call the io class, use:

op.io.VTK.save(network=pn, phases=Hg, filename='first_examples')
jgostick commented 5 years ago

Regarding visualizing throats in paraview, there is an example for this because it is not at all easy.

https://github.com/PMEAL/OpenPNM-Examples/blob/master/IO_and_Visualization/view_vtp_in_paraview.md

xu-kai-xu commented 5 years ago

it seems nice, thanks. i did what the example told, and i find the trick is the filter, which is some kind of amazing to me.

xu-kai-xu commented 5 years ago

the two line code you suggest maybe means to specify that the phases is the invasion_pressure. at first i just copy them before the last line, and it did produce the vtp files, but the import to paraview is fail, which tells "Cannot read point data array "phase|phase_01|properties|pore.invasion_pressure" from PointData in piece 0. The data array in the element may be too short." and i check the new vtp files with the old one, they show difference after the 96489 index (i write some codes to do this, do you think it is right?)

" example_1 = open('first_examples.vtp', 'rt') e_1 = example_1.read() example_1.close()

example_2 = open('examples_2.vtp', 'rt') e_2 = example_2.read() example_2.close()

i = 0 while True: i += 1 if e_1[i] != e_2[i]: print(e_1[i], e_2[i], i) break " BTY, i really want to know how do you write the code in the font, which looks like you write in the compiler. thanks for all the help.

jgostick commented 5 years ago

Hi @xu-kai-xu

  1. I'm not sure why paraview is giving you an error...we use paraview everyday and it works fine.
  2. To get the syntax highlighting put the text inside single back ticks (`) or do multi-lines with triple back ticks (```)
TomTranter commented 5 years ago

Hi @xu-kai-xu Maybe you will find this useful https://www.youtube.com/watch?v=g9yMuU5a5RA&feature=youtu.be Code is in the comment

TomTranter commented 5 years ago

@xu-kai-xu @jgostick I think I know why you are getting this error with the data array may be too short. We have a few percolation algorithms but some assign -infinity to the invasion pressure for boundary pores. It appears that infinity and nan do not get read properly in paraview. I think we should address this issue by specifying a minimum inlet pressure. It is useful to be able to filter out the inlets which is why I think infinity or nan might have been used whereas a value match some of the genuine invasion pressures for other pores. The alternative is to catch infs and nans in the output and change them to some arbitrary very high or negative value

xu-kai-xu commented 5 years ago

@TomTranter @jgostick thanks for the video and code, i will try figuring the problem, and it seems a lot to learn.

TomTranter commented 5 years ago

You just need to replace the bad values like phase[property][numpy.isinf(phase[property])]=9999

xu-kai-xu commented 5 years ago

@TomTranter emm, bad news. i just copy your code in my spyder, but it gives me a syntaxerror. the mysterious code is print('Max Pressure', mip['throat.invasion_pressure'].max()) , the last line, and the feedback is " File "D:\program_files_e\anaconda\lib\site-packages\IPython\core\interactiveshell.py", line 2961, in run_code exec(code_obj, self.user_global_ns, self.user_ns)

File "", line 1, in runfile('E:/xu/python_files/examples/tom\'s advice for openpnm.py', wdir='E:/xu/python_files/examples')

File "D:\program_files_e\anaconda\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 668, in runfile execfile(filename, namespace)

File "D:\program_files_e\anaconda\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile exec(compile(f.read(), filename, 'exec'), namespace)

File "E:/xu/python_files/examples/tom's advice for openpnm.py", line 42 print('Max Pressure', mip['throat.invasion_pressure'].max()) ^ SyntaxError: invalid character in identifier", i have no idea what happened.

TomTranter commented 5 years ago

So it works without that line? I would rename the file something without spaces and apostrophe to see if that helps

xu-kai-xu commented 5 years ago

yes, it works without the last line. the video is so detailed. i feel i understand a little about how the paraview works in visualization. and we use the different version, your is 5.2.0 and mine is 5.6.0-RC2. and i can't find the scaler factor and scale mode in the properties panel, but it seems does not influence the result. thanks very much.

default

default

xu-kai-xu commented 5 years ago

You just need to replace the bad values like phase[property][numpy.isinf(phase[property])]=9999

i think i know what you mean. when i read the openpnm.io.vtk.py files, i find that def save(cls, network, phases=[], filename='', delim=' | ', fill_nans=None): line, tht last Parameter is fill_nans=None, and the documentation about it is that

fill_nans : scalar The value to use to replace NaNs with. The VTK file format does not work with NaNs, so they must be dealt with. The default is None which means property arrays with NaNs are not written to the file. Other useful options might be 0 or -1, but the user must be aware that these are not real values, only place holders.

then i change it to fill_nans=9999, just as you said, then the output works.

xu-kai-xu commented 5 years ago

it should be fill_nans=0, the '9999' value receive that"scalar range is zero", and the '0' value works.

xu-kai-xu commented 4 years ago

@TomTranter @jgostick Hi, these days I am studying and working with openpnm. I used the method Tom told me nearly one year ago. I want to create a simple network with only one same direction throats, like K0206I YFS__)7L066GJRY5

Then when I visualize it with ParaView, I found something confusing me. Tom told me to extract the throats using three filters: cell data to point data, extract surface, tube, which looks like image when I use only two filters: extract surface, tube, the result looks like another way: image The latter is much closer to what I expected.

So I just want to ask the function of the filter "cell data to point data". Which of the above two visualization is better and close to my expected.

I attached the vtp file for your inspection.

Thanks for your help. water_alg.vtp.txt

PS: May you recommend some nice Paraview Tutorials? The official tutorial is too boring.

TomTranter commented 4 years ago

I think it is the celldatatopointdata function doing some rounding errors you have very small throat diameter and maybe this is causing trouble. If you don't do this step you can't scale the throat glyphs. Instead of making the diameter tiny why not use the topotools trim function to actually remove the ones you don't want.

xu-kai-xu commented 4 years ago

@TomTranter , thanks for your advice. I repeat the unidirection model and use openpnm.topotools,trim method to process the non-z-direction throats. The result is very nice. (Sorry for the laaaaaaate reply, and thank you again) image

xu-kai-xu commented 4 years ago

hi, @TomTranter and @jgostick, when I use paraview to visualize throats, I found that the applied filters (cell data to point data --> extract surface --> tube) suggested by Tom is with one little problem. The problem is that the radius of the throats are not visualized to be their real radius, which means the throat's visual result depends on the two pores this throat is connected to.

e.g., when I created a unidirectional network, in which pores and throats shown as follow:

image

this network only has two kinds of pore diameters: 2.5e-5 and 5e-5 m. And only the throats connected by two small pore diameter values are assigned with the value 2.5e-5, which means only one group throats are 2.5e-5(number of 100), and other 800 are 5e-5. The throats above are visualized by only using tube filter directly on "cubic-5.vtp" file. Then, I show the results by applying all three filters on cubic-5.vtp:

image

After applying "cell data to point data" and "extract surface" filters, the throats seem to be changed for visualization. The throat diameter changes too much than the half of size.

For now, I don't know how to visualize it properly. So I just show the throat diameter distribution by color. Do you know how to show them correctly?

jgostick commented 4 years ago

I think you're doing it correctly, but the throat scaling in paraview is not perfect. It's only for visualization, so you manually adjust the the scale factor until it looks right. If you want to do real investigation of pore and throat sizes you should use real data in python, not visualization.

xu-kai-xu commented 4 years ago

Thank you! So you mean, just focus the data processed by python. Paraview is just an assistant work. The real pore and throat relations are reflected in the results of openpnm. Do I understand right?

jgostick commented 4 years ago

In paraview you are just making pictures to help visualize, but the pictures are not 'quantitative'. You basically have a bit of artistic license to make the picture look right. If you want to be sure that throats are correct sizes, a picture is not the right way to do this...you should analyze the numerical arrays instead. So, basically don't stress about getting throat scaling perfect, just approximate. If you are confident in your throat sizes, via numerical analysis, then this is not a trick.

xu-kai-xu commented 4 years ago

@jgostick Thank you!

ma-sadeghi commented 4 years ago

@xu-kai-xu I recently found a way to get throat scaling perfect. Here's what you need to do:

  1. Create a new label in OpenPNM to store throat radius data via

    net["throat.radius"] = net["throat.diameter"] / 2
  2. Apply Shrink filter and set Shrink Factor to 1

  3. Apply CellDatatoPointData filter

  4. Apply ExtractSurface filter

  5. Apply Tube filter

Up to this step, pretty much the same as what you can find in our documentation (except for step 0), here's the extra steps:

  1. In the Tube filter window, select throat.radius for your Scalars
  2. [Optional] Set Number of Sides to 32 for better visualization
  3. Set Radius parameter to 1.
  4. Set Vary Radius parameter to By Absolute Scalar
  5. Set Radius Factor to 1.
  6. Apply!

Here's a screen shot of a network with spacing = 1 and throat diameter = 0.5. image

xu-kai-xu commented 4 years ago

@ma-sadeghi, Thanks a lot! I tried it. The result is fascinating. In the following pictures, I change throats' diameters for a group of tubes crossing the tube direction. Very nice. image