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.

@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.

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).

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 {}

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.

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]

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?

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.