Output the spatial results for specific layer



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
  • Newest first
  • Active threads
  • 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
    out = toughio.read_output("OUTPUT_ELEME.csv")[-1]

    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]
    • 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.

    • 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
      out = toughio.read_output("OUTPUT_ELEME.csv")[-1]
    • 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.

    • 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).

    • Refaat G Hashish My bad, it should be:

      import matplotlib.pyplot as plt
    • 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. 

    • 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
          pointa=(xmin, z, 0.0),
          pointb=(xmax, z, 0.0),
          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"])
    • Keurfon Luu Thank you.

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

    • 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)
          (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:


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

    • 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?

    • Refaat G Hashish Try removing `cpos` argument:


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

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


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