large values for the Pinball env

[ReHoss]

Hello, @jcallaham

By running the Pinball Re=30 or Re=130 with the default initial field, I observe very large values of the observation at the beginning.
It looks like it is independent of the mesh (I tried with a coarse and the fine mesh), see above for 0.2s:

[jcallaham]

By default initial field do you mean all zeros? I would say that's likely expected behavior - I would recommend initializing with a nonzero field if at all possible, whether a snapshot from a previous transient run, a perturbed steady-state solution, mean field, etc.

[ReHoss]

@jcallaham Hello,

The initial field is the default field when you do not provide restart or apply load_checkpoint to the flow. So I guess it is zero.

Thanks to your example on the Cavity, I also tried with a perturbed steady-state flow, and I observed the same behaviour. I will provide you with some evidence.

[jcallaham]

Can you share the code to reproduce in the case of the perturbed steady-state flow? I haven't seen anything that looks like the screenshot you shared when I've done that.

…

On Wed, Sep 11, 2024, 10:45 AM Hosseinkhan Rémy ***@***.***> wrote: @jcallaham <https://github.com/jcallaham> Hello, The initial field is the default field when you do not provide restart or apply load_checkpoint to the flow. So I guess it is zero. I also tried with a perturbated steady-state flow and I observe the same behaviour. — Reply to this email directly, view it on GitHub <#185 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG6KRLTJI6HTBPU5OY6AIIDZWBJSDAVCNFSM6AAAAABN655GU2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNBTHA4DOMJZG4> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[ReHoss]

Just for information:
So last year, I discussed in #94 about an

explosion of the dynamics when applying random control forcing when starting with a null vector field after 40 time steps

Then you answered a fix on the control term has been applied in #132 .

I will come back with a minimal working example (if there is indeed an issue from my side).

[ReHoss]

This code should be self-contained. I changed the signature of compute_steady_state a bit with respect to the example of #186.

You can observe a large amount of energy at the beginning of the run (note that I am a beginner in CFD, and this may be normal to an initiate; however, for data-driven control, the associated value is an anomaly). This could cause issues because the user needs to find stable initial checkpoints for all flows or let the flow evolve for some amount of time before interacting with it. #146 is thus important, I think.

This behaviour is observed for any mesh.

PS: Sorry the code is a bit long, but there are a lot of boilerplate and quite simple functions.

import pathlib
from typing import Type

# import tempfile
import firedrake as fd  # pyright: ignore [reportMissingImports]
import hydrogym.firedrake as hgym

# from tensorboard.compat.proto.types_pb2 import DT_STRING

LIST_MESHES = ["coarse", "medium", "fine"]
LIST_STR_FLOWS = ["cylinder", "pinball", "cavity", "backward_facing_step"]
ACTION_MAGNITUDE = 0.0
N_STEPS = 10
DT = 0.0001


# noinspection DuplicatedCode
def integrate_no_control(gym_env: hgym.FlowEnv, n_steps: int | float) -> None:
    assert (float(n_steps)).is_integer(), "n_steps must be an integer"
    n_steps = int(n_steps)
    array_action = gym_env.action_space.sample() * ACTION_MAGNITUDE
    for i in range(n_steps):
        # gym_env.step(array_action)
        obs, reward, done, info = gym_env.step(array_action)
        print(f"Step {i} - Obs {obs} - Reward {reward} - Done {done} - Info {info}")


def get_hydrogym_flow(name_flow: str) -> Type[hgym.FlowConfig]:
    assert name_flow in LIST_STR_FLOWS, "Invalid flow name"

    if name_flow == "cylinder":
        return hgym.Cylinder
    elif name_flow == "pinball":
        return hgym.Pinball
    elif name_flow == "cavity":
        return hgym.Cavity
    elif name_flow == "backward_facing_step":
        return hgym.Step
    else:
        raise ValueError("Invalid flow name")


def reynolds_curriculum(
        name_flow: str,
        reynolds_number: int,
) -> list[int]:
    assert name_flow in LIST_STR_FLOWS, "Invalid flow name"
    assert reynolds_number > 0, "Reynolds number must be positive"

    list_reynolds_cavity = [500, 1000, 2000, 4000, 7500]
    list_reynolds_pinball = [60, 80, 90, 100, 120, 130]
    list_reynolds_cylinder = [100, 120, 130]

    if name_flow == "cavity" and reynolds_number > min(list_reynolds_cavity):
        # Extract sub-list of Reynolds numbers
        list_curriculum = [
            reynolds for reynolds in list_reynolds_cavity if reynolds < reynolds_number
        ]
        list_curriculum.append(reynolds_number)
    elif name_flow == "pinball" and reynolds_number > min(list_reynolds_pinball):
        # Extract sub-list of Reynolds numbers
        list_curriculum = [
            reynolds for reynolds in list_reynolds_pinball if reynolds < reynolds_number
        ]
        list_curriculum.append(reynolds_number)
    elif name_flow == "cylinder" and reynolds_number > min(list_reynolds_cylinder):
        # Extract sub-list of Reynolds numbers
        list_curriculum = [
            reynolds
            for reynolds in list_reynolds_cylinder
            if reynolds < reynolds_number
        ]
        list_curriculum.append(reynolds_number)
    else:
        list_curriculum = [reynolds_number]

    return list_curriculum


def compute_steady_state(flow: hgym.FlowConfig,
                         reynolds_number: int,
                         name_mesh_resolution: str,
                         stabilization: str = "none", solver_parameters=None):
    if solver_parameters is None:
        solver_parameters = {}
    assert stabilization in [
        "none",
        "gls",
    ], "Invalid stabilization method"  # TODO: Augment stabilization methods
    assert name_mesh_resolution in LIST_MESHES, "Invalid mesh resolution"
    assert reynolds_number > 0, "Reynolds number must be positive"

    # Generate
    solver = hgym.NewtonSolver(
        flow=flow,
        stabilization=stabilization,  # pyright: ignore [reportCallIssue]
        solver_parameters=solver_parameters,
    )

    # Degree of freedom
    dof: int = flow.mixed_space.dim()
    hgym.print(f"Total dof: {dof} --- dof/rank: {int(dof / fd.COMM_WORLD.size)}")

    # Get the Reynolds curriculum
    list_reynolds_curriculum = reynolds_curriculum(
        # Extract flow name from the class name
        name_flow=flow.__class__.__name__.lower(),
        reynolds_number=reynolds_number
    )

    for reynolds_number in list_reynolds_curriculum:
        flow.Re.assign(reynolds_number)
        hgym.print(f"Steady solve at Re={reynolds_number}")
        solver.solve()

    # pressure = flow.p
    # velocity = flow.u
    # vorticity = flow.vorticity()

    # if path_output_data is not None create a temporary directory
    # if path_output_data is None:
    #     path_output_data = tempfile.TemporaryDirectory().name

    # if path_output_data is not None:
    #     # noinspection PyTypeChecker
    #     # Create output directory with Pathlib
    #     pathlib.Path(path_output_data).mkdir(parents=True, exist_ok=True)
    # 
    # with tempfile.TemporaryDirectory() as tmpfile:
    #     if path_output_data is None:
    #         path_output_data = tmpfile
    #     # Save the solution
    #     flow.save_checkpoint(f"{path_output_data}/{reynolds_number}_steady.h5")
    #     # noinspection PyUnresolvedReferences
    #     pvd = fd.output.VTKFile(f"{path_output_data}/{reynolds_number}_steady.pvd")
    #     pvd.write(velocity, pressure, vorticity)
    return flow


if __name__ == "__main__":

    def main():
        name_flow = "pinball"
        reynolds_number = 10
        name_mesh_resolution = "fine"
        stabilization = "none"
        solver_parameters = {"snes_monitor": None}

        # Output directory
        # name_directory = "_".join(
        #     [name_flow, str(reynolds_number), name_mesh_resolution,
        #      stabilization])
        # path_output_data = (f"{PATH_PROJECT_ROOT}/data/steady_state"
        #                     f"/{name_flow}/{name_directory}")

        path_output_data = None

        # Create output directory with Pathlib
        if path_output_data is not None:
            # noinspection PyTypeChecker
            pathlib.Path(path_output_data).mkdir(parents=True, exist_ok=True)

        n_steps = N_STEPS
        dt = DT
        # path_initial_vectorfield =  
        # Integrate

        dict_hydrogym_config = {
            "flow": get_hydrogym_flow(name_flow=name_flow),
            "flow_config": {
                # "restart": path_initial_vectorfield,
                "mesh": name_mesh_resolution,
            },
            "solver": hgym.SemiImplicitBDF, # pyright: ignore [reportAttributeAccessIssue]
            "solver_config": {
                "dt": dt,
            },
            # "callbacks": None,
        }
        gym_env = hgym.FlowEnv(env_config=dict_hydrogym_config)
        
        
        # DEBUG: This should set the steady state solution as the initial condition !
        compute_steady_state(flow=gym_env.flow,  # pyright: ignore
                             reynolds_number=reynolds_number,
                             name_mesh_resolution=name_mesh_resolution,
                             stabilization=stabilization,
                             solver_parameters=solver_parameters)
        # 
        # # Create hydrogym environment

        integrate_no_control(gym_env=gym_env, n_steps=n_steps)


    main()

Ouput

Remark the large values at the beginning.

Total dof: 38106 --- dof/rank: 38106
Steady solve at Re=10
0 SNES Function norm 3.589688336694e+00
1 SNES Function norm 1.176074927565e+00
2 SNES Function norm 1.157146281863e-01
3 SNES Function norm 7.309637971174e-03
4 SNES Function norm 5.920949188727e-05
5 SNES Function norm 1.737109035698e-09
Step 0 - Obs [-19.71199202726886, -4141.936487309615, 4138.994569165966, 43156.93187654883, 52274.125307610004, 52214.78581158276] - Reward -14.764584299574159 - Done False - Info {}
Step 1 - Obs [25.986989800584595, 9.377672598959188, 22.33844326259098, 252.15517078253066, 297.92852592642237, 310.92994517247496] - Reward -0.0861013641881428 - Done False - Info {}
Step 2 - Obs [-51.77645279040027, -35.724117486610275, -27.169124829131523, -382.6068655870437, -472.2187654926844, -492.4637475663498] - Reward 0.1347289378646078 - Done False - Info {}
Step 3 - Obs [39.25704738086866, 16.7749695303573, 30.759032906040208, 357.89431362035606, 426.57085443712697, 445.20574919403657] - Reward -0.12296709172515197 - Done False - Info {}
Step 4 - Obs [-8.356666184689445, -10.509596489031885, 0.24122777655627434, -30.56284778266375, -44.037023144241175, -45.73571306513429] - Reward 0.012033558399203923 - Done False - Info {}
Step 5 - Obs [0.23638784538091678, -5.572398397666319, 5.755475547434987, 39.0979001966516, 40.654540587241456, 42.5604375562592] - Reward -0.012231287834015227 - Done False - Info {}
Step 6 - Obs [0.2849260437916428, -5.461298095436318, 5.65625975803994, 38.74140364581311, 40.66041094320934, 42.58572921936703] - Reward -0.01219875438083895 - Done False - Info {}
Step 7 - Obs [0.2829257303934749, -5.412848974229338, 5.601033561219623, 38.228175515553644, 40.38570436294117, 42.28292935215816] - Reward -0.012089680923065298 - Done False - Info {}
Step 8 - Obs [0.2857196344707256, -5.35366016549313, 5.5324282718672455, 37.7351503080007, 40.12504382037394, 42.00212166142268] - Reward -0.011986231578979732 - Done False - Info {}
Step 9 - Obs [0.28882952932360434, -5.298212809177979, 5.469371118606588, 37.268072530671944, 39.87817411544902, 41.734629113009355] - Reward -0.011888087575913035 - Done False - Info {}

[ReHoss]

Regarding this,

Just for information: So last year, I discussed in #94 about an

explosion of the dynamics when applying random control forcing when starting with a null vector field after 40 time-steps

Then you answered a fix on the control term has been applied in #132 .

I do not observe any explosion with actions sampled uniformly up to twice the maximum magnitude of the action space anymore.

[jcallaham]

Okay, I think what's happening in the case of your code is that the BDF solver is created when the flow field is uniformly zero, so all the historical fields u_prev are initialized to zero. Then you're taking the flow and solving for the steady state, but those old fields are still zero. So when the BDF solver starts, the estimates of the time derivative are large rather than basically zero, and you still get transient behavior.

I guess the assumption in the FlowEnv initialization is that the field passed in as the restart is the one that should actually start integration, and so in your case that would have to be manually passed to the solver. Or you could do something slightly hacky like creating a fresh solver before calling integrate_no_control:

# ... steady-state solve
# Create hydrogym environment
gym_env.solver = hgym.SemiImplicitBDF(gym_env.flow, dt=dt)  # Create a new solver initialized with the steady-state fields
integrate_no_control(gym_env=gym_env, n_steps=n_steps)

Otherwise I would recommend pre-computing the steady state and storing it as a restart and then initializing the FlowEnv with that. Alternatively, if you don't necessarily need the "Env" interface and you want to do some computation up front, what I have usually done is use the hgym.integrate interface, which basically avoids this problem by creating the transient solver after steady-state solving or whatever else is already done.

However, I do think that this should work if you were to call gym_env.solver.reset(), which currently won't do much of anything, so I'll leave the ticket open as a reminder to fix that by setting the u_prev fields to the current value of flow.q[0]

[ReHoss]

Thank you very much.

Otherwise I would recommend pre-computing the steady state and storing it as a restart and then initializing the FlowEnv with that.

This is what I do in practice.

Alternatively, if you don't necessarily need the "Env" interface and you want to do some computation up front, what I have usually done is use the hgym.integrate interface, which basically avoids this problem by creating the transient solver after steady-state solving or whatever else is already done.

Does integrate avoid this problem because the creation of the solver is necessarily done after the flow creation; thus, the correct values are initialised? (I guess it is what you mean, but the steady-state solving made me a bit confused^^)

However, I do think that this should work if you were to call gym_env.solver.reset(), which currently won't do much of anything, so I'll leave the ticket open as a reminder to fix that by setting the u_prev fields to the current value of flow.q[0]

If I get it correctly, the zero initialisation of the previous values creates very large derivatives for the first steps, which causes the observed behaviour?

[jcallaham]

Does integrate avoid this problem because the creation of the solver is necessarily done after the flow creation; thus, the correct values are initialised? (I guess it is what you mean, but the steady-state solving made me a bit confused^^)

Yes exactly. If you create a NewtonSolver and then call .solve(), the steady-state solution will be stored in flow.q. Then when you call integrate, it will initialize all the u_prev fields of the BDF solver to whatever flow.q[0] is, so the time derivative should be approximately zero.

If I get it correctly, the zero initialisation of the previous values creates very large derivatives for the first steps, which causes the observed behaviour?

Yes, the time derivative term in the NS residual is estimated like this for BDF:

    k_idx = k - 1
    u_BDF = sum(beta * u_n for beta, u_n in zip(_beta_BDF[k_idx], self.u_prev))
    alpha_k = _alpha_BDF[k_idx]
    u_t = (alpha_k * u - u_BDF) / h  # BDF estimate of time derivative

So if the u_prev fields are zero, then the time derivative estimate will be way off, so the residual minimization will give a totally different field for the next time step. If the u_prev fields are the steady-state velocity, then u_t ~= 0 and the rest of the weak form residual should also be near zero, so the solution at the next time step should be about the same as the steady-state.

[ReHoss]

Hi @jcallaham ,

I have some questions regarding the behaviour of the flow under control.
At low Reynolds, when injecting a white noise (totally non-continuous) $u(t) \sim Uniform(|\epsilon|)$ for any $t$. Please note that the $\epsilon$ is very low (~1e-2).
I observe that the noise is totally communicated to the dynamics, could you explain me why.
In layman's terms, the noise is not dumped.

Thank you very much!!:)

Illustrations:
Obs
Control
Reward

Cylinder Flow

$Re=120$

The observation $y_0$ is not impacted the obs $y_1$ is fully impacted.

Cavity Flow

$Re=500$


The observation is directly impacted by the noise!

$Re=7500$

The observation seems not to be impacted by the noise!

[jcallaham]

I would say this is likely a result of what you might call "feedthrough". In a linear system

ẋ = Ax + Bu
y = Cx + Du

if the matrix D is not zero then the output depends directly on the input. I've only done this linear analysis for the cylinder, but at least in that case I can confirm that the feedthrough matrix is nonzero. With that in mind, here's how I would explain your results:

• For the cylinder, the observations are (lift, drag). So what you seem to see is that there is a feedthrough term for drag, but not lift.
• For the low-Reynolds cavity there appears to be a weak feedthrough effect. Note that this is a linearly stable flow at this Reynolds number.
• For the higher-Reynolds cavity the feedthrough term appears to be gone. But note that the y-axis of this plot is on the order of 100x magnified compared to the low-Reynolds flow, so there may be a similar feedthrough signal that just doesn't appear in the plot.

You might confirm this by looking for coherence between the input and output at frequencies well above the natural oscillation frequency of the flow. If there is coherence, it would tend to support the idea that there is a feedthrough effect.

Unless you suspect this is a bug, I might go ahead and close this issue - we can continue in the discussion section if you want.

[ReHoss]

@jcallaham Thank you very much. I will read a bit about "Feedthrough" and related concepts soon. Indeed, you can close this issue or, if possible, transform it into a discussion now.

Don't you think that there are possible issues with integration (on the boundary) since the signal $u$ is extremely irregular (nowhere continuous!)?

Best,