🚀[FEA]: Moving bodies with transient simulation

[michaltakac]

Is this a new feature, an improvement, or a change to existing functionality?

New Feature

How would you describe the priority of this feature request

Critical (currently preventing usage)

Please provide a clear description of problem you would like to solve.

Examples for transient simulations in Modulus Sym only consider geometries to be static. There are many usecases in for example turbomachinery or in modeling flow aorund rotating wind turbines, where some parts of the geometry is moved in some way on every step (or the movement can be modeled by some periodic function).

Either showcasing such feature (how to model moving bodies) or implementing this into Modulus Sym would unlock many usecases.

Describe any alternatives you have considered

No response

Additional context

No response

[martin-muzelak]

Needed feature for transient simulation, vote for it!

[Rulino]

This is a really crucial feature. I was thinking about it like a geometry parametrization (e.g. different rotations of a turbine) but I am not sure if this is the correct approach and have not tested it yet. The idea was to train simulator for multiple different rotations of the given geometry and then on inference rotate it by a fixed angle always taking previous time step output as the input. I am not sure if this approach wouldn't cause the convergence to be extremely slow. An example of moving bodies would certainly help!

[DJopek]

It would be extremely helpful since there are a lot of use cases where transient is necessity. Also there is no clear way how to do that. We would really appreciate an example of moving bodies

[michaltakac]

I started working on this in our internal PoC project at DimensionLab, where I'm modeling a simple wind turbine in a wind tunnel and I want it to be rotating. Preliminary code is available here: https://github.com/DimensionLab/transient-rotating-wind-turbine. My goal is to use PINNs without data.

Let's say, these are time parameters together with angular displacement for the turbine:

# time window parameters
time_window_size = 1.0
t_symbol = Symbol("t")
amplitude = np.pi / 12 # in radians
freq = 10.0 # Frequency in [rad/s] 
w = 2.0 * np.pi * amplitude * cos(freq * t_symbol) # Angular displacement 
time_range = {t_symbol: (0, time_window_size)}
nr_time_windows = 10

What I want to achieve is to, when training the model, to give it slightly rotated boundary condition for the wind turbine blades, but to somehow make the model capture how the rotating blades affect the wind...

I've tried 2 approaches:

1. Created custom moving_body fn where I tried doing

def moving_body():
    blades = blades.rotate(angle=w, axis="y", parameterization=Parameterization({"t": ...iteraded specific time...}))
    s = blades.sample_boundary(nr_points=cfg.batch_size.initial_condition, parameterization=Parameterization({"t": ...iterated specific time...}))
    var_to_polyvtk(s, "outputs/wind_turbine/initial_conditions/constraints/bladesBC")
    time_window_net.move_window()

and then call it in

slv = SequentialSolver(
    cfg,
    [(1, ic_domain), (nr_time_windows, window_domain)],
    custom_update_operation= moving_body,
)

But I guess this is wrong way of doing things, as it only re-generated the bladesBC constraint once when the window is moved, and not for every step of the training

2. Created custom PointwiseRotatingBoundaryCondition:

class PointwiseRotatingBoundaryConstraint(PointwiseConstraint):
    """
    Pointwise Constraint applied to boundary/perimeter/surface of rotating geometry.
    For example, in 3D this will create a constraint on the surface of the
    given geometry.

    Parameters
    ----------
    nodes : List[Node]
        List of Modulus Nodes to unroll graph with.
    geometry : Geometry
        Modulus `Geometry` to apply the constraint with.
    angular_displacement : Union[float, sp.Basic]
        The angular displacement of the geometry.
    axis : str
        The axis to rotate the geometry around.
    outvar : Dict[str, Union[int, float, sp.Basic]]
        A dictionary of SymPy Symbols/Expr, floats or int.
        This is used to describe the constraint. For example,
        `outvar={'u': 0}` would specify `'u'` to be zero everywhere
        on the constraint.
    batch_size : int
        Batch size used in training.
    criteria : Union[sp.Basic, True]
        SymPy criteria function specifies to only apply constraint to areas
        that satisfy this criteria. For example, if
        `criteria=sympy.Symbol('x')>0` then only areas that have positive
        `'x'` values will have the constraint applied to them.
    lambda_weighting :  Dict[str, Union[int, float, sp.Basic]] = None
        The spatial pointwise weighting of the constraint. For example,
        `lambda_weighting={'lambda_u': 2.0*sympy.Symbol('x')}` would
        apply a pointwise weighting to the loss of `2.0 * x`.
    parameterization : Union[Parameterization, None], optional
        This allows adding parameterization or additional inputs.
    compute_sdf_derivatives: bool, optional
        Compute SDF derivatives when sampling geometery
    batch_per_epoch : int = 1000
        If `fixed_dataset=True` then the total number of points generated
        to apply constraint on is `total_nr_points=batch_per_epoch*batch_size`.
    quasirandom : bool = False
        If true then sample the points using the Halton sequence.
    num_workers : int
        Number of worker used in fetching data.
    loss : Loss
        Modulus `Loss` module that defines the loss type, (e.g. L2, L1, ...).
    shuffle : bool, optional
        Randomly shuffle examples in dataset every epoch, by default True
    """

    def __init__(
        self,
        nodes: List[Node],
        geometry: Geometry,
        outvar: Dict[str, Union[int, float, sp.Basic]],
        batch_size: int,
        time_window_net: Arch,
        angular_displacement: Union[float, sp.Basic] = 0.0,
        axis: str = "z",
        criteria: Union[sp.Basic, Callable, None] = None,
        lambda_weighting: Dict[str, Union[int, float, sp.Basic]] = None,
        parameterization: Union[Parameterization, None] = None,
        batch_per_epoch: int = 1000,
        quasirandom: bool = False,
        num_workers: int = 0,
        loss: Loss = PointwiseLossNorm(),
        shuffle: bool = True,
    ):
        self.geometry = geometry
        # invar function
        def invar_fn():
            self.geometry = self.geometry.rotate(angle=angular_displacement, axis=axis, parameterization=parameterization)
            return self.geometry.sample_boundary(
                batch_size,
                criteria=criteria,
                parameterization=parameterization,
                quasirandom=quasirandom,
            )

        # outvar function
        outvar_fn = lambda invar: _compute_outvar(invar, outvar)

        # lambda weighting function
        lambda_weighting_fn = lambda invar, outvar: _compute_lambda_weighting(
            invar, outvar, lambda_weighting
        )

        # make point dataloader
        dataset = ContinuousPointwiseIterableDataset(
            invar_fn=invar_fn,
            outvar_fn=outvar_fn,
            lambda_weighting_fn=lambda_weighting_fn,
        )

        # initialize constraint
        super().__init__(
            nodes=nodes,
            dataset=dataset,
            loss=loss,
            batch_size=batch_size,
            shuffle=shuffle,
            drop_last=True,
            num_workers=num_workers,
        )

This, when used looks like this:

bladesBC = PointwiseRotatingBoundaryConstraint(
        nodes=nodes,
        geometry=blades,
        angular_displacement=w,
        axis="y",
        outvar={"u": 0, "v": 0, "w": 0},
        batch_size=cfg.batch_size.initial_condition,
        lambda_weighting={"u": 100, "v": 100, "w": 100},
        parameterization=OrderedParameterization(time_range, key=t_symbol),
    )
    ic_domain.add_constraint(bladesBC, "bladesBC")
    window_domain.add_constraint(bladesBC, "bladesBC")

Although it generates point cloud for boundary condition across whole parameterization space (turbine rotation), which looks okay (white dots in image)... although it does that once and I don't know how to make it so that only one "state" of the rotation is used for particular training step... and also how to make it so that wind is affected by rotating blades?