GazzolaLab / PyElastica

Python implementation of Elastica, an open-source software for the simulation of assemblies of slender, one-dimensional structures using Cosserat Rod theory.
https://www.cosseratrods.org
MIT License
208 stars 107 forks source link

Restart issue while loading a ring-rod #376

Open armantekinalp opened 2 months ago

armantekinalp commented 2 months ago

For restarting a simulation we load the data of the previous simulation using the following load_state function. This function behaves expectedly for straight rods. However, for ring rods since they have periodic boundaries at the ends these boundaries have to by synched by calling the constrain functions. Thus, after loading the states we need to call the boundary conditions before starting integration.

https://github.com/GazzolaLab/PyElastica/blob/cf862d064b40177588286a828e72d68227a76de5/elastica/restart.py#L56-L103

Following changes will fix the issue.


def load_state(simulator, directory: str = "", verbose: bool = False):
    """ 
    Load the rod-state. Compatibale with 'save_state' method. 
    If the save-file does not exist, it returns error. 
    Call this function after finalize method. 

    Parameters 
    ---------- 
    simulator : object 
        Simulator object. 
    directory : str 
        Directory path name. 
    verbose : boolean 

    Returns 
    ------ 
    time : float 
        Simulation time of systems when they are saved. 
    """
    time_list = []  # Simulation time of rods when they are saved. 
    for idx, rod in enumerate(simulator):
        if isinstance(rod, MemoryBlockCosseratRod) or isinstance(
                rod, MemoryBlockRigidBody
        ):
            continue
        path = os.path.join(directory, "system_{}.npz".format(idx))
        data = np.load(path, allow_pickle=True)
        for key, value in data.items():
            if key == "time":
                time_list.append(value.item())
                continue

            if value.shape != ():
                # Copy data into placeholders 
                getattr(rod, key)[:] = value
            else:
                # For single-value data 
                setattr(rod, key, value)

    if not all_equal(time_list):
        raise ValueError(
            "Restart time of loaded rods are different, check your inputs!"
        )

    # Apply boundary conditions. Ring rods have periodic BC, so we need to update periodic elements in memory block
    self.simulator.constrain_values(time)
    # Damping is also part of constrains, so we need to distinguish BC and damping.
    for feature in self.simulator._feature_group_constrain_rates:
        if feature.__name__ == "_constrain_rates":
            feature(time)

    if verbose:
        print("Load complete: {}".format(directory))

    return time

We also need to add relevant test cases, to test if ring rods are loaded correctly.

skim0119 commented 1 month ago

33

We need to save the entire memory block as part of simulation state.