Benchmark 1: Monodisperse Cube
==============================

The script that performs the **Pysammos** Coarse-graining of this benchmark simulation is located at:

``pysammos/benchmarks/monodisperse_cube/run_sweep_pysammos.ipynb``

The plotting of the results against other benchmarks is located at:

``pysammos/benchmarks/monodisperse_cube/plot.ipynb``


Before running Pysammos (Important!)
------------------------------------


1. Contact Data Source 
~~~~~~~~~~~~~~~~~~~~~~~

For benchmarking, the point data reader must use JP’s contact data reader.

**Required code change:**

- Go to the module: ``pysammos.data_read.mfix.point_data``
- Open the function: ``contacts``
- Comment out the default contact reader code and uncomment the code block
  immediately below the comment:

  ``# only to benchmark with JP's contact data reader``

Such that it looks as below:

.. code-block:: python

    # poly_output = InputConnection.GetOutput()
    # F_ij = get_point_data_variable(Force_ij_string, poly_output).astype(np.float32) if Force_ij_string else None
    # Particle_i = get_point_data_variable(Particle_i_string, poly_output).astype(np.float32) if Particle_i_string else None
    # Particle_j = get_point_data_variable(Particle_j_string, poly_output).astype(np.float32) if Particle_j_string else None
    # Contact_ij = get_point_data_variable(Contact_ij_string, poly_output).astype(np.float32) if Contact_ij_string else None

    # only to benchmark with JP's contact data reader
    poly_output = InputConnection.GetOutput().GetCellData(); print("Contact Data loaded as Cell Data")
    F_ij = vtk_to_numpy(poly_output.GetArray(Force_ij_string)).astype(np.float32)
    contact_ids = vtk_to_numpy(poly_output.GetArray(Particle_i_string)).astype(np.float32)
    Particle_i = contact_ids[:, 0]
    Particle_j = contact_ids[:, 1]
    Contact_ij = vtk_to_numpy(poly_output.GetArray(Contact_ij_string)).astype(np.float32)


2. Smoothing Functions c/w 
~~~~~~~~~~~~~~~~~~~~~~~~~~~

The benchmarks use a specific range-to-width ratio.

For consistency, in the function ``calc_cutoff`` of the module
``pysammos.spatial_weights_resolution``, change them to the following:

.. code-block:: python

    if function == 'Lucy':
        c = w
    elif function == 'Gaussian':
        c = 3*w
    elif function == 'HeavySide':
        c = w


4. Kernel Restart
~~~~~~~~~~~~~~~~~

If any code changes are made before re-running the benchmark, you must restart
the notebook.

- Select **Restart & Run All**

This ensures no cached variables are used and results are fully reproducible.


Running the Benchmark
---------------------

1. Sensitivity Setting
~~~~~~~~~~~~~~~~~~~~~~

- The sensitivity parameter must be increased to at least:

  ``sensitivity ≥ 10000``

- This is required to obtain a smooth match across all ``w/d`` values.


2. The Tested CG Widths
~~~~~~~~~~~~~~~~~~~~~~~

The CG widths for which the sweep is carried out must be the following:

.. code-block:: python

    # cg values used by JP:
    cgWidths = np.array([
      0.41666667, 0.5, 0.58333333, 0.66666667, 0.75,
      0.83333333, 0.91666667, 1.0, 1.08333333, 1.16666667,
      1.25, 1.33333333, 1.5, 1.75, 2.0,
      2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 4.0
    ])

Each smoothing function uses a different multiple of the baseline CG widths.

.. code-block:: python

    dp = 0.02  # particle diameter in m

    Lucy_widths = cgWidths * dp * 3
    Gaussian_widths = cgWidths * dp
    HeavySide_widths = cgWidths * dp * 3


3. Reference Point Location
~~~~~~~~~~~~~~~~~~~~~~~~~~~

All benchmark comparisons must be performed at the following reference point:

``P(x,y,z) = [0.0, 0.0, 0.42725]``

This point must be used as the center of a **3×3×3 coarse-graining grid**
by entering it in the ``CG.sweep_CG_widths()`` function.

The index of the point within the grid is printed when running that function.
It tends to be ``index = 13``.

.. code-block:: python

    CG.sweep_CG_widths(
      w_d=eval(f"{cg_func}_widths") / 0.02,
      center=np.array([0.0, 0.0, 0.42725])
    )