Pytough mulgrid object
I'm really interested in using Pytough, especially the write_vtk function available for t2listing objects. So I started recently to dive into the world of python and pytough. After going through the neat documentation I'm still not 100% clear about the relationship between the mulgrid and t2data / t2listing objects. Is there any function / process that I can create a mulgrid object from the Tough2 MESH file or the t2data and t2listing object, that I managed to create with Tough2 files I created so far using PetraSim?
Since ELEME and CONNE blocks contain information regarding the position of the blocks in space, I'm wondering if this is not the main information needed to be able to visualize the 3D data set. Or am I missing something here?
And if I understand correctly I need a mulgrid object in order to use the function t2listing.write_vtk.
The TOUGH2 MESH or input file may (optionally) contain block centres but in general only contains the finite volume information (block volumes, connection areas, distances etc.) needed for the TOUGH simulation.
The mulgrid object contains a complete geometrical description of the mesh, i.e. it defines the positions of all the block vertices. This is needed for visualization (and is useful for other things as well).
Hence it usually isn't possible to create a mulgrid from a t2data object- there just isn't enough information to do it (even if the block centres are present in the t2data grid).
When you create a mesh using PetraSim I understand it stores a geometrical description of the mesh somewhere, but unfortunately it's hidden away in a closed format, so you can't easily convert it to a mulgrid or any other format.
Yes, you do need a mulgrid object to use t2listing write_vtk(). But if your mesh is generated in PetraSim your only option at present is to recreate the mulgrid separately. That may be feasible if your mesh is not too complicated. If your mesh is rectangular (regular or irregular) then it's quite easy to create a mulgrid representing it using the mulgrid rectangular() method.
Note that the mulgrid specification is limited to meshes made up of horizontal layers. The columns can be unstructured, and the upper layers are allowed to be incomplete, to represent topography. But if you have e.g. layers with varying thickness then you won't be able to use a mulgrid for your mesh.
I've been thinking some more about this, and have started writing a script to automatically generate a mulgrid object from a t2grid object, assuming the t2grid represents a rectangular grid. I'll let you know how it goes. If it works reliably I'll include it in PyTOUGH.
The other problem I can see is that the block naming convention and block ordering for e.g. PetraSim grids aren't the same as that for mulgrids. I could possibly make PyTOUGH a bit more flexible in this regard by allowing a 'block mapping' parameter in some of the routines (e.g. t2listing write_vtk()), so that you can map block names in a mulgrid into block names generated by other means like PetraSim. That would allow you to do VTK export of PetraSim results fairly easily, if the mulgrid creation routine also generated the block mapping for you.Reply
I have just committed a bunch of new changes to the PyTOUGH testing branch:
Included in these changes is a new t2grid method, rectgeo(), which can reverse-engineer a mulgrid geometry from a rectangular t2grid.
It also produces a block mapping dictionary, which maps the block names in the geometry to those in the t2grid. You can then pass this block mapping in to various PyTOUGH methods (e.g. write_vtk()) which use both a mulgrid and a t2grid. Using the mapping allows you to use a TOUGH2 grid with any block naming convention in conjunction with a mulgrid geometry.
So now you should be able to export your PetraSim results to VTK for visualization like this:
from t2data import * from t2listing import * dat = t2data('INFILE') geo, blockmap = dat.grid.rectgeo() lst = t2listing('OUTFILE') lst.write_vtk(geo, 'OUTFILE.pvd', dat.grid, blockmap = blockmap)Reply
I managed geo.write_vtk so far that looks good. I tried with a model with fixed boundary elements - where PetraSim sets volume of blocks fixed to 'zero' - first and got the error message:
geo, BM = M07s01dat.grid.rectgeo() File "C:\Program Files (x86)\Python278\lib\site-packages\t2grids.py", line 852, in rectgeo blockmap = block_mapping(geo, self, ob, nblks, atmos_volume) File "C:\Program Files (x86)\Python278\lib\site-packages\t2grids.py", line 819, in block_mapping mapping[geoblkname] = blk.name AttributeError: 'NoneType' object has no attribute 'name'
But when switching to a model with all blocks enabled and not fixed, it works fine. Still playing with lst.write_vtk command.Reply
I have a quite similar problem. I wanted to visualize the change in temperature or pressure resulting from the operation of geothermal injection and production well (difference between the initial values and T (orP) after a given time). I thought that Pytough would be perfect to visualise contour of the changes in 3D, but I did not succeed.
I use Petrasim to run my TOUGH2 simulatoins and I want to use the input and output files generated by Petrasim. My mesh is big (about 170 000 cells) and the layers follow the topography. Each layer has a constant thickness but is not horizontal. Is this a problem? Any idea how I could plot the contours in 3D?
I can provide the Petrasim files if needed.
Thank you in advance!Reply
The 'mulgrid' geometry in PyTOUGH is limited to layers that are horizontal and constant thickness. (While I can certainly see the attraction of having layers following stratigraphy etc., I'm not certain that meshes with variable thickness layers are a good idea, because the vertical connections are no longer orthogonal to the block faces, which is a requirement of the TOUGH2 numerical formulation.)
So PyTOUGH won't be able to create a geometry corresponding to your mesh, unfortunately. Can you not do that kind of visualization within Petrasim itself?Reply
In fact each layer has a constant thickness but it varies from one layer to another. I simply use one layer for the top surface which is based on the topography and then the others layers corresponds to this layer shifted downwards by a given value (for example 50 m if I want to have a 50 m thick layer). However, when doing this the layers are indeed not horizontal.
It does not seem to be a problem for the TOUGH2 calculation.
Unfortunately with Petrasim you can visualize the P and T fields at a given time step but you cannot visualize the difference of T or P between two different time steps... which is very important for me.Reply
Though I guess PyTough can do that issue with delta P and T faster/cleaner, here is my little tweak with Petrasim:
Manipulation of mesh.csv (which can be tedious if a lot of time steps are written out):
- Use pressure at different time steps from mesh.csv to calculate delta P/T (using Excel or other)
- Add results as additional column in mesh.csv (or create mesh.csv only containing "timestep") and open 3D output again. (care must be taken with formatting of csv)
From the figure it looks like your mesh isn't just a tilted rectangular mesh- does each layer basically have the same shape as the surface, but translated downwards? (i.e. does the bottom of the mesh have the same shape as the top?)
If so, then this isn't representable via a mulgrid geometry. And the connections in your mesh would almost never be orthogonal to the block faces. Personally I'd be wary of using such a mesh with TOUGH2. This paper explains why.
This is focused on local grid refinement but the issues regarding orthogonality of connections to block faces are the same. In fact, it looks from the picture of your mesh that the vertical block faces wouldn't even be planar, which I don't think TOUGH2 would like either.
We used to use meshes with non-orthogonal connections in TOUGH2 until we found it caused some quite spurious results in certain situations.
If you decide to keep using that mesh though, I suspect what you're doing with creating mesh.csv files and reading them back into PetraSim is probably not a bad approach. Though it could probably be streamlined using a bit of Python scripting. If you find yourself needing to manipulate CSV files or other tabular data I'd recommend looking at the Pandas library. You could easily script the calculation of temperature differences and writing the new mesh.csv file. Much easier than using excel!Reply