Output the spatial results for specific layer

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.

Refaat G Hashish My bad, it should be:

import matplotlib.pyplot as plt

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

# Read output file
out = toughio.read_output("OUTPUT_ELEME.csv")[-1]

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

# Add data to points
for k, v in out.data.items():
    if k not in {"X", "Y", "Z"}:
        mesh.add_point_data(k, v)

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),
)

Refaat G Hashish The easiest way to export to an Excel speadsheet is to use the module pandas.

import pandas
import toughio

out = toughio.read_output("OUTPUT_ELEME.csv")[-1]
df = pandas.DataFrame(out.data)
df = df[df["Z"] == -12.5]
df = df.sort_values(by=["X"])
df.to_excel("output.xls")

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.add_mesh(mesh.to_pyvista(), n_colors=2, show_edges=True)
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.

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.