Yade-batch is being indeterministic (gtsPfacet)

Hello all,

Please use the following code for post-processing

Kind regards,
Rohit

# ------------------------------------------------------------------------------------------------------ python script
import json
import numpy as np
from matplotlib import pyplot as plt
import plotly.graph_objects as go

import sys, os
sys.path.append(".")
current_path = os.getcwd()
# --------------------------------------------------------------------------- setting the directory
os.chdir(current_path)

# --------------------------------------------------------------------------- selecting the data files
files = os.listdir()
files_sel = []
for i in files:
    if i.find('json') > 0:
        files_sel.append(i)

# --------------------------------------------------------------------------- importing data
simulation_results = {}

for i in files_sel:
    with open(i, 'r', encoding="utf-8") as jsonfile:
        data = json.load(jsonfile)
        simulation_results[data['ID']] = data

for i in list(simulation_results.keys()):
    wx = simulation_results[i]['data']['wx']
    plt.plot(wx)

Hello,

Are you *sure* each batch run is running on one core? You don't appear to be setting the core number anywhere.

You can add a column to your batch file titled "!OMP_NUM_THREADS". Set that column equal to 1, re run, and then let us know if you still are seeing indeterminism.

Cheers,

Robert

Dear Robert,

I tried adding the column "!OMP_NUM_THREADS" and I set all the values to one (given below). I am still seeing indeterminism. Please let me know if there is anything that I may have missed.

Kind regards,
Rohit K. John

PC information
Yade version: 2021-08-24.git-d2e8416
Ubuntu 20.04.3 LTS

# -------------------------------------------------------------------------------------------------------------------------------------------- parameters.table
ID !OMP_NUM_THREADS
1 1
2 1
3 1
4 1
5 1
6 1
7 1
8 1
9 1
10 1

What is the output of yade-batch itself?

Hi,
This is intriguing.
You can print this from within the script to be sure:
print("OMP_NUM_THREADS=",os.environ['OMP_NUM_THREADS'] )
Then inspect the log files and see if it's 1.

I can't imagine why only the batch execution would be non deterministic.
Are you sure non-batch on one core is deterministic?

Bruno

Hello,

I have put the content of the .log file at the end of this comment. In a previous, albeit more complicated script, I had found that even running the same script multiple times, even with a single core (I used -j=1), yielded different results.

The simulation monitors the angular velocity of pfacet clump. I tried running the script, the one I have posted above, three times and the results diverge.

I apologise if the title appeared to be misleading.

Kind regards,
Rohit K. John

# ----------------------------------------------------------------------------------------------------------- YADE LOG file for batch
Welcome to Yade 2021-08-24.git-d2e8416
Using python version: 3.8.10 (default, Jun 2 2021, 10:49:15)
[GCC 9.4.0]
Warning: no X rendering available (see https://bbs.archlinux.org/viewtopic.php?id=13189)
XMLRPC info provider on http://localhost:21000
TCP python prompt on localhost:9000, auth cookie `ksdsyu'
Running script to_LP_main.py

=================== JOB SUMMARY ================
id : ID=1,OMP_NUM_THREADS=1
status : 0 (OK)
duration: 00:03:28
command : YADE_BATCH=parametric.table:2 DISPLAY= /usr/local/bin/yade-2021-08-24.git-d2e8416 --threads=1 --nice=10 -x to_LP_main.py> ./logs/to_LP_main.py.2.log 2>&1
started : Wed Sep 22 15:24:21 2021
finished: Wed Sep 22 15:27:49 2021

Hi,

So it is the script performance that is being indeterministic, and the question is why. Maybe you should split your simulation into a couple of steps to determine which part introduces randomness. I can't tell which part of the code is responsible for that. But if you are able to isolate this, you could:
a) solve the problem at the source,
b) find a workaround (e.g. run the indeterministic part just once, and then load and modify the simulation within the batch mode).

I bet on the gtsPFacet() function. I ran your script a couple of times in a row with option -j1 (only the script, no simulation running). Then I inspected pfacets by simply checking one of them:

***
O.bodies[6].state.se3
***
And once in two to five simulations, the result changes. This change is probably due to the different order of facets processing, but still indeterministic.

 I am not very familiar with gts surfaces, but I think that inside the gts.read()* function is some kind of algorithm at the end to 'patch' all the numerical inaccuracies.
*which is part of gtsPFacet()[1]

Best wishes,
Karol

[1] https://yade-dem.org/doc/_modules/yade/gridpfacet.html#gtsPFacet

Dear Karol,

Thanks for spending your time to find the probable cause and proposing a viable solution. I shall look into this and get back as soon as possible.

Kind regards,
Rohit K. John

Dear Karol,

You are correct. The problem lies with the gtsPfacet() function. It is introducing randomness somehow.

I tried writing my own function to read a GTS file and create a pfacet object. When I use this function (create_GTSPfacet()) instead of gtsPfacet(), all simulations in batch-mode produce the same result. I will paste the code below for future reference. If I find any other issues regarding this I shall post them here.

Kind regards,
Rohit K. John

# ---------------------------------------------------------------------------------------- Python
# --------------------------------------------------------------------------------------------- split_bySpaces
def split_bySpaces(input_string):
    """
    Uses spaces (' ') to split a string into a list of substrings

    Args:
        input_string (str): input string

    Returns:
        list: list of substrings
    """
    string = ['']
    for i in input_string:
        if i == ' ':
            string = string + ['']

        else:
            string[-1] = string[-1] + i

    if string[-1] == '':
        string = string[:-1]

    return string

# --------------------------------------------------------------------------------------------- GTSParser
class GTSParser:
    """
    Parses the GTS file and extracts the data pertaining to vertices, edges and faces
    """
    # ------------------------------------------------------------------
    def __init__(self, input_str):
        """
        Initialises and parses the dataset

        Args:
            input_str (str): The string data from GTS file
        """
        self.raw_data = input_str

        self._split_rawByLines()
        self._get_3DMetadata()
        self._get_vertices()
        self._get_edges()
        self._get_faces()

    # ------------------------------------------------------------------
    def _split_rawByLines(self):
        """
        Splits the GTS string into lines based on '\n' character. This is stored
        in list
        """
        self._data_line_splitted = ['']
        for i in self.raw_data:
            if i == '\n':
                self._data_line_splitted = self._data_line_splitted + ['']
            else:
                self._data_line_splitted[-1] = self._data_line_splitted[-1] + i

    # ------------------------------------------------------------------
    def _get_3DMetadata(self):
        """
        Extracts the metadata (number of vertices, edges and faces) from the first line
        """
        self._metadata_string = self._data_line_splitted[0]
        metadata_string_splitted = split_bySpaces(self._metadata_string)
        self._metadata = {}
        for i in range(3):
            self._metadata[metadata_string_splitted[-i-1]] = int(metadata_string_splitted[i])

    # ------------------------------------------------------------------
    def _get_vertices(self):
        """
        Extracts the vertices (Their location) and stores them in a variable
        """
        no_vertices = self._metadata['GtsVertex']
        start_idx = 1
        vertex_data = self._data_line_splitted[start_idx:start_idx + no_vertices]
        vertices = []
        for i in vertex_data:
            vertices = vertices + [[float(j) for j in split_bySpaces(i)]]

        self._vertices = vertices

    # ------------------------------------------------------------------
    def _get_edges(self):
        """
        Extracts the edges (the vertices which form them) and stores them in a variable
        """
        no_vertices = self._metadata['GtsVertex']
        no_edges = self._metadata['GtsEdge']
        start_idx = 1 + no_vertices
        edge_data = self._data_line_splitted[start_idx:start_idx + no_edges]
        edges = []

        for i in edge_data:
            edges = edges + [[int(j) for j in split_bySpaces(i)]]

        self._edges = edges

    # ------------------------------------------------------------------
    def _get_faces(self):
        """
        Extracts the faces (the edges which form them) and stores them in a variable
        """
        no_vertices = self._metadata['GtsVertex']
        no_edges = self._metadata['GtsEdge']
        np_faces = self._metadata['GtsFace']
        start_idx = 1 + no_vertices + no_edges
        face_data = self._data_line_splitted[start_idx:start_idx + np_faces]
        faces = []
        for i in face_data:
            faces = faces + [[int(j) for j in split_bySpaces(i)]]

        self._faces = faces

    # ------------------------------------------------------------------
    @property
    def vertices(self):
        return self._vertices

    @property
    def edges(self):
        return self._edges

    @property
    def faces(self):
        return self._faces

# --------------------------------------------------------------------------------------------- GTSPfacet2
class GTSPfacet2:
    '''
    Creates a pfacet mesh from a GTS file
    '''
    def __init__(
        self,
        gts_file_path,
        radius,
        shift,
        scale,
        wire,
        fixed,
        color,
        material_nodes,
        material_ext
        ):
        """
            Initialises and creates the pfacet
        """

        self._gts_file_path = gts_file_path
        self._radius = radius
        self._shift = Vector3(shift)
        self._scale = scale
        self._wire_Q = wire
        self._fixed_Q = fixed
        self._color = color
        self._material_nodes = material_nodes
        self._material_ext = material_ext

        self._read_file()
        self._parsed_gts = GTSParser(self._gts_data)
        self._generate_Nodes()
        self._generate_Edges()
        self._generate_Faces()

    # ------------------------------------------------------------------
    def _read_file(self):
        """
        Reads the gts into a variable
        """
        with open(self._gts_file_path, 'r') as file:
            self._gts_data = file.read()

    # ------------------------------------------------------------------
    def _generate_Nodes(self):
        """
        Generates the nodes from the vertices defined in the GTS file
        """
        self._node_ids = []
        self._node_id_to_gts_node_idx = {}

        # creating nodes
        for idx,i in enumerate(self._parsed_gts.vertices):
            self._node_ids.append(
                O.bodies.append(
                    gridNode(
                        self._scale*Vector3(i) + self._shift,
                        self._radius,
                        # dynamic = None,
                        wire = self._wire_Q,
                        fixed = self._fixed_Q,
                        color = self._color,
                        highlight = False,
                        material = self._material_nodes
                    )))

            # GTS vertex index -> ID
            self._node_id_to_gts_node_idx[idx + 1] = self._node_ids[-1]

    # ------------------------------------------------------------------
    def _generate_Edges(self):
        """
        Generates the edges/GridConnections using the GTS file
        """
        self._cyl_ids = []
        self._edge_id_to_gts_edge_idx = {}

        # creating Edges
        for idx,k in enumerate(self._parsed_gts.edges):
            id1 = self._node_id_to_gts_node_idx[k[0]]
            id2 = self._node_id_to_gts_node_idx[k[1]]
            self._cyl_ids.append(
                O.bodies.append(
                    gridConnection(
                        id1,id2,
                        self._radius,
                        wire = self._wire_Q,
                        color = self._color,
                        highlight = False,
                        material = self._material_ext,
                        mask = -1
                )))
            # GTS edge index -> ID
            self._edge_id_to_gts_edge_idx[idx + 1] = self._cyl_ids[-1]

    # ------------------------------------------------------------------
    def _generate_Faces(self):
        """
        Generates the pfacet faces defined the GTS file
        """
        self._pfacet_ids = []
        self._face_id_to_gts_face_idx = {}

        # creating faces
        for idx, j in enumerate(self._parsed_gts.faces):
            id1 = self._edge_id_to_gts_edge_idx[j[0]]
            id2 = self._edge_id_to_gts_edge_idx[j[1]]
            id3 = self._edge_id_to_gts_edge_idx[j[2]]
            pfacetCreator4(
                id1,id2,id3,
                pfIds = self._pfacet_ids,
                wire = self._wire_Q,
                fixed = self._fixed_Q,
                color = self._color,
                material = self._material_ext,
                mask = -1
                )

            # GTS face index -> ID
            self._face_id_to_gts_face_idx[idx + 1] = self._pfacet_ids[-1]

    # ------------------------------------------------------------------
    @property
    def node_ids(self):
        return self._node_ids

    @property
    def cyl_ids(self):
        return self._cyl_ids

    @property
    def face_ids(self):
        return self._pfacet_ids

# --------------------------------------------------------------------------------------------- create_GTSPfacet
def create_GTSPfacet(
        gts_file_path,
        radius,
        shift,
        scale,
        wire,
        fixed,
        color,
        materialNodes,
        material
        ):
    """
    Creates a pfacet object from the input GTS file.

    Args:
        gts_file_path (str): The path to the GTS file
        radius (float): Radius of the pfacet nodes.
        shift (Vector3): The location of the pfacet object
        scale (Float): The amout the object is to be scaled
        wire (bool) :
        fixed (bool): Is it fixed
        color (list): Color of the object
        materialNodes (str): The material name of the nodes
        material (str): The material name of the pfacet and the grids

    Returns:
        (list, list, list): ids of the nodes, cylinder, pfacets
    """
    object = GTSPfacet2(
        gts_file_path,
        radius,
        shift,
        scale,
        wire,
        fixed,
        color,
        materialNodes,
        material
        )

    return object.node_ids, object.cyl_ids, object.face_ids