Contour plots

[jousheinfo]

Hi there!

I am trying to plot a contour figure typical from the oil and gas industry, let's say depth(y-axis)-time(x-axis)-saturation. Please find attached the code below. I would very much appreciate your support. I am mixing 2 python libraries: toughio and pytough. I am attaching the MESH file. Thank you in advance for your time and consideration!

import toughio import numpy as np from t2listing import *

#Reading the mesh file mesh_inp = toughio.read_mesh("MESH", file_format="tough")

#Reading the output file lst_gas = t2listing('injection_Pwh_50_gas_CO2_salinity_pytough.listing')

#time = [] pressure_gas = [] gas_saturation_gas = [] blocks_gas = [] time_gas = [] #Getting the blocks name for block_index in mesh_inp["elements"].keys(): blocks_gas.append(block_index) #print(blocks)

#Getting the properties for the blocks
[(t_gas,p_gas), (t_gas, temp_gas), (t_gas, d_gas), (t_gas, sg_gas)] = lst_gas.history([('e', block_index, 'P'), ('e', block_index, 'T'), ('e', block_index, 'DG'), ('e', block_index, 'SG')]) #e stands for 'element'

#Getting the properties' values for all blocks for the whole simulation time
pressure_gas.append(p_gas)
gas_saturation_gas.append(sg_gas)

#Getting time values for all blocks for the whole simulation time
time_gas.append(t_gas)

#Getting time values for the whole simulation t_gas = t_gas

#Getting the properties' values for wellbore cells only for the whole simulation time P_gas = np.concatenate((pressure_gas[:50], pressure_gas[100:101])) SG_gas = np.concatenate((gas_saturation_gas[:50], gas_saturation_gas[100:101]))

#Getting the Z-coordinate values for wellbore cells only z_gas_wellbore = np.concatenate((z_array[:50], z_array[100:101]))

#Getting the time values for the first time steps up to n (according to the Python file CO2_injection_model_comparion_results.py) t_gas_specified = t_gas[:24]

#Getting the properties' values for the first time steps up to n (according to the Python file CO2_injection_model_comparion_results.py) P_gas_time_specified = [arr[:24] for arr in P_gas] SG_gas_time_specified = [arr[:24] for arr in SG_gas]

#Plotting plt.tricontourf(t_gas_specified, np.linspace(0, 2500, 24), SG_gas_time_specified) #plt.tricontour(x1_95_reservoir.values.flatten(), T1_95_reservoir.values.flatten(), SG1_95_reservoir.values.flatten()) plt.colorbar() plt.show() MESH.txt

[acroucher]

hi, can you tell me what the problem is with your script? I can't run it without the listing file (and I actually don't have toughio installed either).

However I can see a number of problems with your script just from reading it:

• in the lst_gas.history() call, you are using block_index to identify the block, but this is not defined except as a loop variable in the for loop above it.
• also in the history() call you are using the same variable t_gas as the returned time array for all the different element properties (P, T, etc.). This is OK if you are sure all those variables are present at all times in the listing (they might not be if e.g. you are using 'short' output in an AUTOUGH2 model).
• I don't understand what you're doing with the np.concatenate() calls, which seem to be selecting the first 50 results and then adding the 100th one on the end

[jousheinfo]

The problem was with the dimensions when plotting plt.tricontourf. But, I fixed it. Thanks! However, I have another question.

In the code below, my wellbore radius is 0.3 m and the reservoir radius is 6000 m. When running the code, I get in block.centre[0] much higher than this (2.747e+04). I still cannot understand it. Morever, if I write well_radius = 0.3, in the input file, it will give me a value equal to 0.15. As far as I know, the block.centre[0] in the input file is the centre of the cell. That's why I changed the value to 0.6. Is this approach correct? I am confused. Thanks!

well_radius = 0.6 reservoir_radius = 6000 total_num_blocks = 53

values = np.logspace(np.log10(well_radius), np.log10(reservoir_radius), total_num_blocks, endpoint=True) dat = t2data() dat.title = 'Cold_CO2_liquid_injection' dat.grid = t2grid().radial(dr, dz, origin = [0,0])

[jousheinfo]

I forgot to ask, how to change the code in order to actually have a radius reservoir of 6000 m? I very much appreciate your help.

will this help with the desired goal?

Calculate the logarithmic increment factor

increment_factor = np.log10(reservoir_radius / well_radius) / (total_num_blocks - 1)

Generate the dr values on a logarithmic scale

dr = well_radius * np.power(10, increment_factor * np.arange(total_num_blocks))

create TOUGH2 input data:

dat = t2data() dat.title = 'Cold_CO2_liquid_injection' dat.grid = t2grid().radial(dr, dz, origin = [0,0])

[acroucher]

On the first question, you're creating an array called values that looks like the radial block sizes, but then passing in another array called dr to t2grid.radial(). Could that be the problem?

Regarding the well radius, the first value in your dr array is the radial thickness of the first block. The centre of that first block will by default be at the middle of the block, so at half the specified radial thickness. If you aren't including the well in the model you usually use the origin parameter in t2grid.radial() to shift the grid origin to the well radius (as in the example on the PyTOUGH wiki).

If you do want to include the well in the model then you don't need to shift the grid origin but you may need to adjust the centre of the first block in your t2grid so it's at zero radius and you get the connection distance to the second block equal to the well radius.

If you want the grid radius of 6000 m you can either compute the expansion factor of the grid sizes so that you end up with the correct total radius, or (more easily) just truncate or extend the last radial block thickness.