data_viewer.py

# /// script
# requires-python = ">=3.11"
# dependencies = [
#     "h5py==3.14.0",
#     "lxml==6.0.1",
#     "numpy==2.3.3",
#     "pyvista==0.46.3",
#     "vtk==9.5.1",
# ]
# ///

import marimo

__generated_with = "0.15.2"
app = marimo.App(width="medium")


@app.cell
def _():
    import marimo as mo
    import numpy as np
    import h5py
    from lxml import etree as ET
    import pyvista as pv
    import vtk
    return ET, h5py, mo, np, pv, vtk


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
    # Data Viewer
    ## Configure setup
    ⚙️ **How it is done**:
    Get the data (`data.h5` and `mesh.h5`) from the [DaRUS repository](https://darus.uni-stuttgart.de/dataset.xhtml?persistentId=doi:10.18419/DARUS-5120) and place them in the same directory as this notebook.
    """
    )
    return


@app.cell(hide_code=True)
def _(mo):
    fn_mesh_input = mo.ui.text(label='Name of mesh file:', value='mesh.h5')
    fn_data_input = mo.ui.text(label='Name of data file:', value='data.h5')
    mo.hstack([fn_mesh_input, fn_data_input])
    return fn_data_input, fn_mesh_input


@app.cell(hide_code=True)
def _(fn_data_input, fn_mesh_input):
    fn_mesh = fn_mesh_input.value
    fn_data = fn_data_input.value
    return fn_data, fn_mesh


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
    ## Select a parameter configuration
    ⚙️ **How it is done**:
    Choose the diffusion coefficient for the bulk $D_\Omega$ as well as the GB thickness $2h$ and the overall size of the RVE $l$ in arbitrary units. The sliders for the relation between the in-plane and normal GB diffusion coefficients $D_\Vert/D_\Omega$ and $D_\perp/D_\Omega$ indicate which simulation results are available within the data base (based on the corresponding 2d representation). Select any configuration of your choice.
    """
    )
    return


@app.cell(hide_code=True)
def _(mo):
    D_bulk_input = mo.ui.text(label=r'$D_\Omega$', value='1.')
    h_input = mo.ui.text(label=r'$h$', value='0.0031622776601683794')
    l_input = mo.ui.text(label=r'$l$', value='1')
    mo.vstack([D_bulk_input,
              h_input,
              l_input])
    return D_bulk_input, h_input, l_input


@app.cell(hide_code=True)
def _(D_bulk_input, h_input, l_input, mo, np):
    # Define reference quantities
    D_0 = float(D_bulk_input.value)
    l_0 = float(l_input.value)

    pi_bulk_log = np.log10(float(D_bulk_input.value) / D_0)
    pi_h_log = np.log10(float(h_input.value)) / l_0
    pi_l_log = np.log10(float(l_input.value)) / l_0

    Pi_para_min = -5
    Pi_para_max = 0
    Pi_perp_min = 0
    Pi_perp_max = 5

    pi_para_min = Pi_para_min + pi_bulk_log + pi_l_log - pi_h_log
    pi_para_max = Pi_para_max + pi_bulk_log + pi_l_log - pi_h_log
    slider_para = mo.ui.slider(start=pi_para_min, stop=pi_para_max, step=0.25, value=0, label=r'$\log_{10}(D_\Vert/D_\Omega)$')

    pi_perp_min = Pi_perp_min + pi_bulk_log + pi_h_log - pi_l_log
    pi_perp_max = Pi_perp_max + pi_bulk_log + pi_h_log - pi_l_log
    slider_perp = mo.ui.slider(start=pi_perp_min, stop=pi_perp_max, step=0.25, value=0, label=r'$\log_{10}(D_\perp/D_\Omega)$')
    return (
        D_0,
        l_0,
        pi_bulk_log,
        pi_h_log,
        pi_l_log,
        pi_para_max,
        pi_para_min,
        pi_perp_max,
        pi_perp_min,
        slider_para,
        slider_perp,
    )


@app.cell(hide_code=True)
def _(
    D_0,
    mo,
    pi_para_max,
    pi_para_min,
    pi_perp_max,
    pi_perp_min,
    slider_para,
    slider_perp,
):
    mo.vstack([mo.hstack([mo.md(f"Choose value in: [{pi_para_min}, {pi_para_max}]"), slider_para, mo.md(f"Selected: {slider_para.value}"), mo.md(f"$D_\Vert={D_0*10**(slider_para.value):.3f}$")]),
                mo.hstack([mo.md(f"Choose value in: [{pi_perp_min}, {pi_perp_max}]"), slider_perp, mo.md(f"Selected: {slider_perp.value}"), mo.md(f"$D_\perp={D_0*10**(slider_perp.value):.3f}$")])])
    return


@app.cell(hide_code=True)
def _(pi_bulk_log, pi_h_log, pi_l_log, slider_para, slider_perp):
    pi_bulk = 10**pi_bulk_log
    pi_h = 10**pi_h_log
    pi_l = 10**pi_l_log

    pi_para_log = slider_para.value
    pi_para = 10**pi_para_log
    pi_perp_log = slider_perp.value
    pi_perp = 10**pi_perp_log

    Pi_para = pi_para_log - pi_bulk_log + pi_h_log - pi_l_log
    Pi_perp = pi_perp_log - pi_bulk_log + pi_l_log - pi_h_log
    return Pi_para, Pi_perp, pi_bulk, pi_h, pi_l, pi_para, pi_perp


@app.cell(hide_code=True)
def _(D_0, Pi_para, Pi_perp, l_0, mo, pi_bulk, pi_h, pi_l, pi_para, pi_perp):
    mo.md(
        rf"""
    **5d parameter representation** $(D_\Omega, D_\Vert, D_\perp, h, l)$:
    ({pi_bulk*D_0:.2f}, {pi_para*D_0:.2f}, {pi_perp*D_0:.2f}, {pi_h*l_0:.2f}, {pi_l*l_0:.2f})

    **Dimensionless 5d parameter representation** $(\pi_\Omega, \pi_\Vert, \pi_\perp, \pi_h, \pi_l)$:
    ({pi_bulk:.2f}, {pi_para:.2f}, {pi_perp:.2f}, {pi_h:.2f}, {pi_l:.2f})

    **Corresponding dimensionless 2d parameter representation** $(\log_{{10}}(\Pi_\Vert), \log_{{10}}(\Pi_\perp))$:
    ({Pi_para}, {Pi_perp})
    """
    )
    return


@app.cell(hide_code=True)
def _(Pi_para, Pi_perp, np):
    assert Pi_para in np.linspace(-5,0,21), 'No data available for this configuration of Pi_para!'
    assert Pi_perp in np.linspace(0,5,21), 'No data available for this configuration of Pi_perp!'
    return


@app.cell(hide_code=True)
def _(Pi_para, Pi_perp, fn_data, h5py):
    sim_idx_match = -1
    with h5py.File(fn_data, 'r') as f:
        for sim_idx in f.keys():
            if sim_idx != 'general':
                if (f[sim_idx].attrs['Pi_parallel'] == Pi_para) and (f[sim_idx].attrs['Pi_perpendicular'] == Pi_perp):
                    sim_idx_match = int(sim_idx)

    assert sim_idx_match != -1, 'No matching full-field data found!'
    return (sim_idx_match,)


@app.cell(hide_code=True)
def _(mo, sim_idx_match):
    mo.md(rf"""**Index of matching full-field data in the database**: {sim_idx_match:03}""")
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
    ## Visualize full-field solutions
    🎯 **Goal**:
    Visualize the results for the 2d-parameter configuration selected above. Note that data for the bulk and GB domain are provided on separate meshes. The number at the beginning of the dataset names (0,1,2) indicates the direction of the imposed gradient.

    ⚙️ **How it is done**:

    - After selecting a parameter configuration above, mesh data and respective simulation data are linked via an `xdmf`-file. Simply select a file name and press the button to generate the required file.
    - For a thorough investigation of all data this file can be directly opened in ParaView.
    - Alternatively, select a mesh (bulk or GB) and one of the point or cell data fields within the notebook. A separate window will open for an interactive 3d visualization. Note that the new window needs to be closed in order to proceed within the notebook.

    🧮 **Rescaling of data**:
    Note that some of the displayable data fields depend on the 5d parameter set. Therefore, a rescaling of the field magnitude needs to be applied for certain configurations (namely changes of $h$ and $D_\perp$). If this applies the rescaling factors are displayed below next to the unscaled fields.

    🦺 **Important**:

    - The `xdmf`-file needs to be regenerated after every update in the parametric setup (click the button).
    - Note that the `xdmf`-contains relative links to the `hdf5`-files.
    """
    )
    return


@app.cell(hide_code=True)
def _(mo):
    fn_xdmf_input = mo.ui.text(label='Name of xdmf file:', value='visu.xdmf')
    fn_xdmf_input
    return (fn_xdmf_input,)


@app.cell(hide_code=True)
def _(ET, fn_data, fn_mesh, h5py):
    def ExportGeometryAndMesh(domain, name_domain, name_vx, name_el):
        """Link geometry and mesh to root-ET.Element for generating xdmf-file.
        For the individual layers of cohesive interface elements nodal information from other layers (master) can be reused. Only cell data
        differs here.

        Args:
            domain (ET.SubElement): Domain that meshes and data are built on.
            name_domain (string): Name of the domain. Possible values: 'bulk', 'gb'.
            name_vx (string): Name of the mesh (nodal part). Possible values: 'node_positions'.
            name_el (string): Name of the mesh (element part). Possible values: 'elements', 'elements_-1', 'elements_0', 'elements_+1'.

        Returns:
            grid (ET.SubElement): Grid object of the added mesh.
            n_vx (int): Number of nodes in the added mesh.
            n_el (int): Number of elements in the added mesh.
        """

        with h5py.File(fn_mesh, 'r') as F:
            n_vx = F[f'{name_domain}/{name_vx}'][:].shape[0]
            el = F[f'{name_domain}/{name_el}']
            n_el = el[:].shape[0]
            n_vx_per_el = el.attrs['nodes_per_element']
            el_type =  el.attrs['element_type']

        # Define grid with name
        grid = ET.SubElement(domain, "Grid", Name="mesh_"+name_domain, GridType="Uniform")

        # Add the topology (elements)
        topology = ET.SubElement(grid, "Topology", TopologyType=el_type, NumberOfElements=str(n_el))
        # Reference to the HDF5 file and the mesh data
        ET.SubElement(topology, "DataItem", DataType="Int", Dimensions=f'{n_el} {n_vx_per_el}', Format="HDF").text = f'{fn_mesh}:/{name_domain}/{name_el}'

        # Add the geometry (nodes)
        geometry = ET.SubElement(grid, "Geometry", GeometryType="XYZ")
        # Reference to the HDF5 file and the coordinates data
        ET.SubElement(geometry, "DataItem", DataType="Float", Dimensions=f'{n_vx} 3', Format="HDF").text = f'{fn_mesh}:/{name_domain}/{name_vx}'

        return grid, n_vx, n_el


    def ExportData(grp, grid, n_vx, n_cells):
        """Link nodal and cell data from data_file to grid from root-ET.Element for generating xdmf-file.
        Names of data sets are assumed as 'vx_data' and 'cell_data'.

        Args:
            grp (string): Name of group in which data is located in HDF5-file. Possible values are: 'bulk', 'gb'.
            grid (ET.SubElement): Grid object of the added mesh.
            n_vx (int): Number of nodes in the mesh (corresponds to number of nodal data).
            n_cells (int): Number of elements in the mesh (corresponds to number of integration point data).
        """

        path_vx = grp + '/vx_data'
        path_cell = grp + '/cell_data'

        with h5py.File(fn_data, 'r') as F:

            # Link nodal data
            if path_vx in F:
                vx_data_keys = list(F[path_vx].keys())

                # print('Nodal visualization data in ' + grp + ' available for the following data sets: ' + ', '.join(str(key) for key in vx_data_keys))

                center = "Node"

                for key in vx_data_keys:
                    dims = F[f'{path_vx}/{key}'][:].shape

                    if len(dims) == 1 or dims[1] == 1:
                        attribute_type = "Scalar"
                        dim = 1
                    else:
                        attribute_type = "Vector"
                        dim = dims[1]

                    # Add the attribute (assuming scalar data)
                    attribute = ET.SubElement(grid, "Attribute", Name=key, AttributeType=attribute_type, Center=center)
                    # Reference to the HDF5 file and the results data
                    ET.SubElement(attribute, "DataItem", DataType="Float", Dimensions=f'{n_vx}{(dim>1)*(f" {dim}")}', Format="HDF").text = fn_data+':/'+path_vx+'/'+key

            # Link cell data
            if path_cell in F:
                cell_data_keys = list(F[path_cell].keys())

                # print('Cell visualization data in ' + grp + ' available for the following data sets: ' + ', '.join(str(key) for key in cell_data_keys))

                center = "Cell"

                for key in cell_data_keys:
                    dims = F[f'{path_cell}/{key}'][:].shape

                    if len(dims) == 1 or dims[1] == 1:
                        attribute_type = "Scalar"
                        dim = 1
                    else:
                        attribute_type = "Vector"
                        dim = dims[1]

                    # Add the attribute (assuming scalar data)
                    attribute = ET.SubElement(grid, "Attribute", Name=key, AttributeType=attribute_type, Center=center)
                    # Reference to the HDF5 file and the results data
                    ET.SubElement(attribute, "DataItem", DataType="Float", Dimensions=f'{n_cells}{(dim>1)*(f" {dim}")}', Format="HDF").text = fn_data+':/'+path_cell+'/'+key
    return ExportData, ExportGeometryAndMesh


@app.cell(hide_code=True)
def _(mo):
    export_visu = mo.ui.run_button(label="Create visualization")
    export_visu
    return (export_visu,)


@app.cell(hide_code=True)
def _(
    ET,
    ExportData,
    ExportGeometryAndMesh,
    export_visu,
    fn_xdmf_input,
    mo,
    sim_idx_match,
):
    fn_xdmf = fn_xdmf_input.value

    if export_visu.value:
        # Create the XDMF root element
        xdmf_root = ET.Element("Xdmf", Version="2.0")

        # Create the domain
        domain = ET.SubElement(xdmf_root, "Domain")

        # Domain and mesh information for bulk
        _grid_bulk, _n_vx_bulk, _n_el_bulk = ExportGeometryAndMesh(domain, 'bulk', 'node_positions', 'elements')
        # Nodal and cell data for bulk
        ExportData(f'{sim_idx_match:03}/bulk', _grid_bulk, _n_vx_bulk, _n_el_bulk)
        ExportData('general/bulk', _grid_bulk, _n_vx_bulk, _n_el_bulk)

        # Domain and mesh information for GB
        _grid_GB, _n_vx_GB, _n_el_GB = ExportGeometryAndMesh(domain, 'gb', 'node_positions', 'elements_-1')
        # Nodal and cell data for GB
        ExportData(f'{sim_idx_match:03}/gb', _grid_GB, _n_vx_GB, _n_el_GB)
        ExportData('general/gb', _grid_GB, _n_vx_GB, _n_el_GB)

        # Write the XDMF file
        with open(fn_xdmf, "wb") as _xdmf_f:
            _xdmf_f.write(ET.tostring(xdmf_root, pretty_print=True))

        mo.output.replace('Visualization updated!')
    return (fn_xdmf,)


@app.cell(hide_code=True)
def _(D_0, mo, np, pi_h, pi_perp):
    mo.output.replace(mo.md('**Rescaling required**:'))

    if not np.isclose(pi_h, 10**(-2.5)):
        scaling_factor = pi_h / (10**(-2.5))
        mo.output.append(mo.md(f'''
            `X_GRAD_PERP_Y` with {scaling_factor:.4f}\n
            `X_NET_INFLUX` with {scaling_factor:.4f}\n
            `X_FIRST_ORDER_MODE_ABS` with {scaling_factor:.4f}\n
            `X_SECOND_ORDER_MODE` with {scaling_factor**2:.4f}'''))
    elif not np.isclose(pi_perp*D_0, 1):
        scaling_factor = pi_perp*D_0
        mo.output.append(mo.md(f'''
            `X_NET_INFLUX` with {scaling_factor:.4f}'''))
    else:
        mo.output.append(mo.md('None.'))
    return


@app.cell(hide_code=True)
def _(export_visu, fn_xdmf, pv, vtk):
    if export_visu.value:
        reader = vtk.vtkXdmfReader()
        reader.SetFileName(fn_xdmf)
        reader.Update()

        vtk_data = reader.GetOutputDataObject(0)
        multi_block = pv.wrap(vtk_data)
    return (multi_block,)


@app.cell(hide_code=True)
def _(export_visu, mo, multi_block):
    dropdown_mesh = None
    if export_visu.value:
        dropdown_mesh = mo.ui.dropdown(options=multi_block.keys())
    dropdown_mesh
    return (dropdown_mesh,)


@app.cell(hide_code=True)
def _(dropdown_mesh, mo, multi_block):
    tabs = None
    if dropdown_mesh != None:
        if dropdown_mesh.value != None:
            mesh = multi_block[dropdown_mesh.value]

            vx_dropdown = mo.ui.dropdown(options=multi_block[dropdown_mesh.value].point_data.keys())
            cell_dropdown = mo.ui.dropdown(options=multi_block[dropdown_mesh.value].cell_data.keys())

            tabs = mo.ui.tabs({
                "point data": vx_dropdown,
                "cell data": cell_dropdown
            })
    tabs
    return cell_dropdown, mesh, tabs, vx_dropdown


@app.cell(hide_code=True)
def _(cell_dropdown, dropdown_mesh, mesh, pv, tabs, vx_dropdown):
    if tabs != None:
        if dropdown_mesh.value != None:
            match tabs.value:
                case "point data":
                    field = vx_dropdown.value
                case "cell data":
                    field = cell_dropdown.value
            if field != None:
                plotter = pv.Plotter(notebook=False)
                plotter.add_mesh(mesh, scalars=field, cmap='turbo')
                plotter.show_axes()
                plotter.show()
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
    ## Compute effective diffusivity tensor
    🎯 **Goal**:
    Compute the effective diffusivity tensor for the 5d-parameter configuration selected above.

    ⚙️ **How it is done**:
    Results are automatically updated in real-time based on the selected setup above.
    """
    )
    return


@app.cell(hide_code=True)
def _(
    D_0,
    fn_data,
    fn_mesh,
    h5py,
    l_0,
    mo,
    np,
    pi_bulk,
    pi_h,
    pi_l,
    pi_para,
    pi_perp,
    sim_idx_match,
):
    # Load all relevant data
    with h5py.File(fn_mesh, 'r') as _f:
        v_bulk_0 = _f['bulk/volume_ref'][()]
        surf_GB_0 = _f['gb/surface_area_ref'][()]
        struct_tensor_2 = _f['gb/structural_tensor_2'][:]

    with h5py.File(fn_data, 'r') as _f:
        grp = _f[f'{sim_idx_match:03}/integrated_grad']
        G_int_bulk = np.hstack((grp['0_integrated_grad_bulk'][:][:,None], grp['1_integrated_grad_bulk'][:][:,None], grp['2_integrated_grad_bulk'][:][:,None]))
        G_int_para = np.hstack((grp['0_integrated_grad_gb_para'][:][:,None], grp['1_integrated_grad_gb_para'][:][:,None], grp['2_integrated_grad_gb_para'][:][:,None]))
        G_int_perp = np.hstack((grp['0_integrated_grad_gb_perp'][:][:,None], grp['1_integrated_grad_gb_perp'][:][:,None], grp['2_integrated_grad_gb_perp'][:][:,None]))

    G_bar = np.eye(3)

    v_GB_0 = 2 * surf_GB_0 * l_0
    f_GB = v_GB_0 / (pi_l/pi_h * v_bulk_0 + v_GB_0)

    G = ((1. - f_GB) * (G_bar + (1/(pi_l*v_bulk_0) * G_int_bulk))
            + f_GB * ((np.eye(3) - 1/surf_GB_0 * struct_tensor_2) @ G_bar + 1/(2*pi_l*surf_GB_0) * G_int_para)
            + f_GB * 1/(2*pi_h*surf_GB_0) * G_int_perp)

    Q = D_0 * (((1. - f_GB) * pi_bulk * (G_bar + (1/(pi_l*v_bulk_0) * G_int_bulk))
        + f_GB * pi_para * ((np.eye(3) - 1/surf_GB_0 * struct_tensor_2) @ G_bar + 1/(2*pi_l*l_0*surf_GB_0) * G_int_para)
        + f_GB * pi_perp * 1/(2*pi_h*l_0*surf_GB_0) * G_int_perp))

    D_eff = Q @ np.linalg.inv(G)

    mo.output.replace(D_eff)
    return


@app.cell(hide_code=True)
def _():
    return


if __name__ == "__main__":
    app.run()