0

# Output the spatial results for specific layer

Dears,

I would like to know if there is a specific way to print out the results (spatial results at a specified time) for a specific layer of grid cells. Since my model is a 2D radial model (350 * 10 grid cells) with 10 layers, the output results (OUTPUT_ELEME) include the results for all grid cells at the specified output times. In this case, I have to filter the results for the remaining layers to get the results for the layer of interest which takes a long time.

15replies Oldest first
• Oldest first
• Popular
• You can use a Python module of mine (https://github.com/keurfonluu/toughio).

To install, you need a Python distribution (e.g. Anaconda) and then run in Windows console:

`pip install toughio[full] --user`

To read the last time step of OUTPUT_ELEME.csv in Python:

```import toughio

```

To extract data (e.g. temperature) at specific depths (e.g. 0.0 <= z < 10.0):

```import numpy

idx = numpy.logical_and(out.data["Z"] >= 0.0, out.data["Z"] < 10.0)
temp = out.data["TEMP"][idx]
```
Like
• Keurfon Luu I have installed anaconda navigator. I have launched the Qt Console and typed the following

`pip install toughio[full] --user`

I got the following result. I would like to know if I follow your steps.

Like
• Keurfon Luu How to output the results versus X at a specific Z ? and What is the rule followed to specify the time step in the following code.

```import toughio

Like
• Refaat G Hashish The function `toughio.read_output` returns a list of all time steps in your file. In Python, `[-1]` refers to the last element in the list ([0] being the first element, [1] the second... a negative integer counts backward).

Let's assume that the targeted layer is at depth z = -5.0 m (meaning the z-component of the elements is equal to -5.0):

```idx = out.data["Z"] == -5.0  # Find indices of elements at z = -5.0 m
x = out.data["X"][idx]       # x-component of target elements
p = out.data["PRES"][idx]    # Pressure of target elements
```

Make sure that these outputs are correctly sorted:

```idx = numpy.argsort(x)
x = x[idx]
p = p[idx]
```

Plotting can be done using the standard plotting library:

```import matplotlib as plt

plt.plot(x, p)
```

But I usually prefer importing the results on a mesh and plot along a line. If you send me your file OUTPUT_ELEME.csv, I can show you how.

Like
• Keurfon Luu Thanks very much for your help. I followed the steps and they work fine, however, I find the following error while plotting. I am using Spyder (Python 3.7). I have attached the excel file for the results (OUTPUT ELEME.csv).

Like
• Refaat G Hashish My bad, it should be:

`import matplotlib.pyplot as plt`
Like
• Keurfon Luu It works fine. Thanks. I would like to know if the stored data can be exported to a text file or excel file so that I can reformat them.

Like
• Refaat G Hashish Here is a script to over a line in a more general fashion.

First, we create a fake mesh from the cloud points and import the output data:

```import numpy
import toughio

# Create a fake mesh from points
points = numpy.c_[out.data["X"], out.data["Z"]]
mesh = toughio.meshmaker.triangulate(points)

for k, v in out.data.items():
if k not in {"X", "Y", "Z"}:
```

Then, we use the module pyvista to post-process the mesh:

```import pyvista

# Convert the mesh to a mesh for pyvista
mesh = mesh.to_pyvista()

# Plot over a line
xmin, xmax = 0.0, 2000.0
z = -10.0
mesh.plot_over_line(
pointa=(xmin, z, 0.0),
pointb=(xmax, z, 0.0),
scalars="SAT_G",
resolution=5000,
title="Saturation profile (Z = {} m)".format(z),
ylabel="Gas saturation",
figsize=(10, 5),
)
```
Like
• Refaat G Hashish The easiest way to export to an Excel speadsheet is to use the module pandas.

```import pandas
import toughio

df = pandas.DataFrame(out.data)
df = df[df["Z"] == -12.5]
df = df.sort_values(by=["X"])
df.to_excel("output.xls")
```
Like
• Keurfon Luu Thank you.

Like
• Keurfon Luu Is it possible to visualize the grid blocks obtained using MESHMAKER in Python?

Like
• Refaat G Hashish Once your fake mesh has been created, you can simply run:

`mesh.plot(notebook=False, show_edges=True)`

Note that because the mesh has been created from the element coordinates in your output file using Delauney triangulation, the fake mesh is made of triangles and not quads as generated by MESHMAKER.

Alternatively, `toughio` can also create a radial mesh from scratch:

```import toughio

dr = [0.3]
while sum(dr) < 100000.0:
dr.append(dr[-1] * 1.1)

dz = 5 * [10.0]
origin_z0 = -sum(dz)
mesh = toughio.meshmaker.cylindric_grid(dr, dz, origin_z0, "sand")
mesh.set_material("well", xlim=(0.0, 0.3), zlim=(origin_z0, 0.0))
```

This should be similar to the mesh you created in your previous posts (I think). To visualize the mesh generated, you can use `pyvista`:

```import pyvista

p = pyvista.Plotter(notebook=False)
p.view_xz()
p.show(cpos=[
(61.18984873150231, -255.23078247814956, -25.601392414026044),
(61.18984873150231, -255.23078247814956, -25.601392414026044),
(0.0, 0.0, 1.0),
])
```

To write the MESH file for TOUGH, simply call:

`mesh.write_tough()`

which will write the MESH file in the current working directory.

Like
• Keurfon Luu Thank you so much for your reply. I have followed the steps, however, I got the graph as follows without showing the discretizations.  Is it possible to use the MESH file generated by MESHMAKER to plot the grid blocks?

Like
• Refaat G Hashish Try removing `cpos` argument:

`p.show()`

MESHMAKER does not output the cell connectivity, so you can only resort to point triangulation to create a fake mesh similarly to output mapping.

Like
• Keurfon Luu

I just wanted to first thank you for developing toughio. Also, I have a little bit of a problem using it for a MINC problem. mesh = toughio.meshmaker.triangulate(points) gives me an error (please see the image). toughio-export command gives me the same error. I really appreciate it if you could let me know what you think about the error.

Best s,

Pedram.

Like
Like Follow
• 6 mths agoLast active
• 15Replies
• 155Views
• 3 Following