Creating structured grid with dipping angle

[keurfonluu]

Hello,

I plan to create structured mesh in a dipping structure using meshmaker. When the dipping angle is zero, I use the following to generate the mesh: mesh = toughio.meshmaker.structured_grid(dx, dy, dz)

If I want to assign a dipping angle, what command should I use?

Originally posted by @saharbakhshian in https://github.com/keurfonluu/toughio/issues/127#issuecomment-1684093017

[keurfonluu]

To generate a mesh with a dipping angle, you can either "physically" rotate your mesh, or simply rotate the direction of the gravity vector and provide it when writing the MESH file. In both cases, the BETAX values calculated are the same, however, for the former approach, the ISOT values are all 1 while for the latter, the ISOT values are preserved. For isotropic permeabilities, both methods work great, otherwise, rotating the gravity vector is recommended.

import numpy as np
import toughio

# Rotation matrix
theta = np.deg2rad(-16.0)
Ry = np.array(
    [
        [np.cos(theta), 0.0, np.sin(theta)],
        [0.0, 1.0, 0.0],
        [-np.sin(theta), 0.0, np.cos(theta)],
    ]
)

# Generate structured grid
dx = 10 * [10.0]
dy = [1.0]
dz = 5 * [10.0]
mesh = toughio.meshmaker.structured_grid(dx, dy, dz)

# Rotate gravity vector
gravity = Ry @ [0.0, 0.0, -1.0]
mesh.write_tough("MESH1", gravity=gravity)

# Rotate mesh
mesh.points = mesh.points @ Ry
mesh.write_tough("MESH2")

# Compare calculated angles
mesh1 = toughio.read_mesh("MESH1")
mesh2 = toughio.read_mesh("MESH2")
angles1 = [connection["gravity_cosine_angle"] for connection in mesh1["connections"].values()]
angles2 = [connection["gravity_cosine_angle"] for connection in mesh2["connections"].values()]
isot1 = [connection["permeability_direction"] for connection in mesh1["connections"].values()]
isot2 = [connection["permeability_direction"] for connection in mesh2["connections"].values()]
print("Angles are the same:", np.allclose(angles1, angles2))  # True
print("Permeability directions are the same:", np.allclose(isot1, isot2))  # False

[iunshiya]

hello， I want to know how to asgin the initial conditions cell by cell according to the dip coordinate. Can you show me the toughio code which will read the actual z value after rotating the mesh or gravity vector of each cell from the MESH file,and calculate T=T0+gradTz and P=P0+gradPz, then write an INCON block?

[geowomen]

import pyvista as pv
import toughio
parameters = toughio.read_input("postfermeture_Scellement.in") ## inmput
save_file = toughio.read_input("./Calcul_final/Tough2_operating/save") # a save file from the previous simulaions
# select element and modify the initial conditions
for elem in parameters["elements"].keys():
    if parameters['elements'][elem]['material'] == 1 :
        if parameters['elements'][elem]['volume'] > 0.01 :
            parameters["initial_conditions"][elem]["values"] = [101300.0, 0.0, 0.0, 0.0, 10.2, 20.0]
    if parameters['elements'][elem]['material'] == 2 :
        print("yes i go here",elem)
        parameters["initial_conditions"][elem]["values"] = [101300.0, 0.0, 0.0, 0.0, 10.51, 20.0]
    if parameters['elements'][elem]['material'] == 3 :
        print('yes yes',parameters["initial_conditions"][elem])
        list_10_3 = [elem for elem in parameters["initial_conditions"].keys() 
                     if parameters['elements'][elem]['material'] == 3 and parameters["initial_conditions"][elem]["values"][-2] == 10.3]
        if elem in list_10_3:
            print(save_file["initial_conditions"][elem])
            parameters["initial_conditions"][elem]["values"] = [101300.0, 0.0, 0.0, 0.000, 10.2, 20.0]
    if parameters['elements'][elem]['material'] in {4,5,6,7} :
        pass

write file : toughio.write_input("./Calcul_final/Tough2_postclosure/Nouvelles_conditions.dat", parameters)

in order to give you an idea

[iunshiya]

Thanks for your reply and I'll give it a try.

[keurfonluu]

More generally, if you have a MESH file with at least the Z coordinates, you can assign hydrostatic and geothermal initial conditions as follows (e.g., for a EOS3, EOS5... with 0.0 being the initial mass fraction of component 2):

parameters = tougio.read_input("MESH")
parameters["initial_conditions"] = {
    k: {
        "values": [
            P0 + gradP * v["center"][2],
            0.0,
            T0 + gradT * v["center"][2],
        ],
    }
    for k, v in parameters["elements"].items()
}
toughio.write_input("INCON", parameters, block="incon")