PHAREHUB / PHARE

💫 Parallel Hybrid Particle In Cell code with Adaptive mesh REfinement
https://phare.readthedocs.io
GNU General Public License v3.0
71 stars 24 forks source link

add time averaged hierarchies #851

Open nicolasaunai opened 5 months ago

nicolasaunai commented 5 months ago

Possible POC below. Main issue being that hierarchies can change in patch number or positions or shape in time, thus not all points can be averaged blindly. It is possible, since time averaging usually is performed over many times but small time intervals, that hierarchies don't change during the averaging time interval. But there is no general guarantee.

The simplest and most expected behavior then is to return a PatchHierarchy with only patches that are constant over the averaging interval.

def time_avg(hier):
    """
    return a PatchHierarchy where patch datas are averaged over time.
    If there is only one time, the hierarchy is returned as is.
    """

    nbr_times = len(hier.times())

    # no need to average if there is only one time
    if nbr_times == 1:
        return hier

    # we first create the returned hierarchy with the same metadata as the input hierarchy
    # but with the average time and initialized with the first patchlevel
    # the data of the first time is copied into the new hierarchy
    # and the averaging is then done.
    avg_time = hier.times().astype(float).mean()
    first_time_lvls = hier.levels(hier.times()[0])

    new_levels = {}
    for ilvl, lvl in first_time_lvls.items():
        new_levels_patches = []
        for patch in lvl.patches:
            pdatas = {}
            for pdname, pd in patch.patch_datas.items():
                pdatas[pdname] = FieldData(
                    pd.layout, pdname, data=pd.dataset[:] / float(nbr_times)
                )
            new_levels_patches.append(Patch(pdatas, patch.id, attrs=patch.attrs))
        new_levels[ilvl] = PatchLevel(ilvl, new_levels_patches)

    avg_hier = PatchHierarchy(
        new_levels,
        hier.domain_box,
        refinement_ratio,
        time=avg_time,
        data_files=hier.data_files,
        selection_box=hier.selection_box,
    )

    for time in hier.times():
        for ilvl in hier.levels(time):
            lvl = hier.level(ilvl, time)
            avg_lvl = avg_hier.level(ilvl, avg_time)
            for avg_patch, patch in zip(avg_lvl, lvl.patches):
                for avg_pd, pd in zip(
                    avg_patch.patch_datas.values(), patch.patch_datas.values()
                ):
                    avg_pd.dataset[:] += pd.dataset[:] / float(nbr_times)

    return avg_hier