writes an Fluent ip file from an MCNP meshtal file

mesh_tal_to_fluent.py

# writes an Fluent ip file from an MCNP meshtal file

import typing
import numpy as np
import pprint


def get_mesh_info(input_filename: str) -> dict:
    """Converts an MCNP meshtal file to a dictionary containing the mesh information for each mesh tally in the file

    Args:
        input_filename (str): input input_filename of MCNP meshtal file

    Returns:
        dict: a dictionary containing the mesh information for each mesh tally in the file
    """
    mesh_info = {}
    found_mesh_tally = False
    lines_to_capture = 0
    current_mesh_tally_number = None

    print(f"Reading MCNP meshtal info from file {input_filename}")
    with open(input_filename, "r") as file:
        for i, line in enumerate(file, start=1):
            if found_mesh_tally:
                mesh_info[current_mesh_tally_number]["description"] = line.strip()
                found_mesh_tally = False

            if lines_to_capture > 0:
                if lines_to_capture == 3:
                    mesh_info[current_mesh_tally_number][
                        "x_boundaries"
                    ] = line.strip().split()[2:]
                elif lines_to_capture == 2:
                    mesh_info[current_mesh_tally_number][
                        "y_boundaries"
                    ] = line.strip().split()[2:]
                elif lines_to_capture == 1:
                    mesh_info[current_mesh_tally_number][
                        "z_boundaries"
                    ] = line.strip().split()[2:]
                    mesh_info[current_mesh_tally_number]["number_of_mesh_elements"] = (
                        (len(mesh_info[current_mesh_tally_number]["x_boundaries"]) - 1)
                        * (
                            len(mesh_info[current_mesh_tally_number]["y_boundaries"])
                            - 1
                        )
                        * (
                            len(mesh_info[current_mesh_tally_number]["z_boundaries"])
                            - 1
                        )
                    )
                lines_to_capture -= 1

            if "Mesh Tally Number" in line:
                current_mesh_tally_number = int(line.split()[-1])
                mesh_info[current_mesh_tally_number] = {
                    "description": "",
                    "x_boundaries": [],
                    "y_boundaries": [],
                    "z_boundaries": [],
                }
                found_mesh_tally = True

            if "Tally bin boundaries:" in line:
                lines_to_capture = 3

            # Check for the specific line and store its line number
            if line.strip() == "X         Y         Z     Result     Rel Error":
                if (
                    current_mesh_tally_number
                ):  # Ensure we are within a mesh tally section
                    mesh_info[current_mesh_tally_number]["data_start_line"] = i

    # Read the data for each mesh tally
    for mesh_id in mesh_info.keys():

        content = np.genfromtxt(
            "3E9_global_n_LOST1E4_1cm.t",
            skip_header=mesh_info[mesh_id]["data_start_line"],
            max_rows=mesh_info[mesh_id]["number_of_mesh_elements"],
        )

        mesh_info[mesh_id]["x_coords"] = content[:, 0]
        mesh_info[mesh_id]["y_coords"] = content[:, 1]
        mesh_info[mesh_id]["z_coords"] = content[:, 2]
        mesh_info[mesh_id]["result"] = content[:, 3]
        mesh_info[mesh_id]["result_error"] = content[:, 4]

    return mesh_info


def write_fluent_ip_file(
    mesh_info: dict,
    filename: str = "fluent_ip_from_mcnp_mctal.ip",
    scaling: float = 100.0,
    uds_numbers: typing.Union[None, typing.Iterable] = None,
):
    """Converts a dictionary containing the mesh information for each mesh tally in the file to a Fluent IP file

    Args:
        mesh_info (dict): A dictionary containing the mesh information for each mesh tally in the file
        filename (str, optional): The filename to write. Defaults to "fluent_ip_from_mcnp_mctal.ip".
        scaling (float, optional): A scaling factor to apply to the coordinates. Defaults to 100.0.
        uds_numbers (typing.Union[None, typing.Iterable], optional): A list of user defined tags to assign to each mesh, can be an iterable of ints, strings, floats. Defaults to None.
    """

    with open(filename, mode="w") as file:

        first_mesh = mesh_info[list(mesh_info.keys())[0]]
        file.write("3\n")  # first line is the interpolation_file_version
        file.write("3\n")  # second line is the dimension
        # number_points (the ip format assumes all the fields are on the same mesh)
        file.write(f"{first_mesh['number_of_mesh_elements']}\n")
        file.write(f"{len(mesh_info.keys())}\n")

        for i, mesh_id in enumerate(mesh_info.keys()):
            if uds_numbers:
                file.write(f"uds-{uds_numbers[i]}\n")
            else:
                file.write(f"uds-{mesh_id}\n")
        # writing the coordinates, special case for first coordinate to include opening bracket
        for coord_name in ["x_coords", "y_coords", "z_coords"]:
            for i, coord in enumerate(first_mesh[coord_name]):
                if i == 0:
                    file.write(f"({coord/scaling}\n")
                else:
                    file.write(f"{coord/scaling}\n")
            file.write(")\n")
        # writing the values, special case for first value to include opening bracket
        for mesh_id, mesh_values in mesh_info.items():
            for i, field_value in enumerate(mesh_values["result"]):
                if i == 0:
                    file.write(f"({field_value}\n")
                else:
                    file.write(f"{field_value}\n")
            file.write(")\n")
    print(f"Fluent IP file written to {filename}")


# Example minimal usage, write all the meshes to a file with defaults
mesh_info = get_mesh_info(
    input_filename="3E9_global_n_LOST1E4_1cm.t",
)
write_fluent_ip_file(
    mesh_info=mesh_info,
    filename="fluent_ip_from_mcnp_mctal_default.ip",
)


# Example usage, write all the meshes to a file
mesh_info = get_mesh_info(
    input_filename="3E9_global_n_LOST1E4_1cm.t",
)
write_fluent_ip_file(
    mesh_info=mesh_info,
    filename="fluent_ip_from_mcnp_mctal.ip",
    scaling=100,
    uds_numbers=[1, 2, 3],
)

# Example usage, how to write just two of the three of the mesh tallies
mesh_info = get_mesh_info(
    input_filename="3E9_global_n_LOST1E4_1cm.t",
)
first_mesh_key = list(mesh_info.keys())[0]
first_mesh_info = {1: mesh_info[first_mesh_key]}
write_fluent_ip_file(
    mesh_info=first_mesh_info,
    filename="fluent_ip_from_mcnp_mctal_subset.ip",
    scaling=100,
    uds_numbers=[1, 2, 3],
)