Code: sandbox

Page with all code information, that will be split later into 3 pages

Under construction

Conventions in the code

Before every method name, the _ prefix is used whenever the method is not intended to be used outside of the body of the class.
U, P (capital) refer to the full field $U(x,t)$, while u, p (small) refer to the perturbation field $u'(x,t)$ (see Numerical details). For boundary conditions, BC, bc follow the same convention.

Code: basic use

Basic use

1️⃣ Choose a use-case: inherit the `FlowSolver` abstract class

The simulation revolves around the abstract class FlowSolver that implements core features such as loading the mesh, defining the function spaces & trial/test functions, variational formulations, numerical schemes and solvers, handling the time-stepping and exporting fields and timeseries. The class is abstract as it does not implement a simulation case per se, but only provides utility for doing so. It features two abstract methods, that are redefined for each use-case:

_make_boundaries provides a definition and naming of the boundaries of the mesh in a pandas DataFrame.

@abstractmethod
def _make_boundaries(self) -> pd.DataFrame:
    pass

The expected DataFrame has the following simple structure:

boundaries_as_df = pandas.DataFrame(
    index=boundaries_names_as_list: list[str], 
    data={"subdomain": subdomains_as_list: list[dolfin.SubDomain]}
)

⚠️ For each use-case, the user is expected to provide a mesh in xdmf format (see Third-party tools), that is compatible with their definition of boundaries.

_make_bcs provides a description of the boundary conditions on the boundaries defined above, in a dedicated class BoundaryConditions containing two lists.

@abstractmethod
def _make_bcs(self) -> BoundaryConditions:
    pass

BoundaryConditions is a utility class that contains two list fields: bcu (velocity boundary conditions for the perturbation field) and bcp (pressure boundary conditions for the perturbation field). See below:

@dataclass
class BoundaryConditions:
    bcu: list[dolfin.DirichletBC]
    bcp: list[dolfin.DirichletBC]

We give two examples with the code (the flow past a cylinder, and the flow over an open cavity) that inherit from FlowSolver: they are respectively CylinderFlowSolver and CavityFlowSolver.

2️⃣ Attach sensors and actuators to an instance of a `FlowSolver` subclass

In order to perform sensing and actuation (with the objective to close the loop), two dedicated abstract classes are proposed: Sensor and Actuator. Both these classes implement behaviors common to all sensors or actuators. They are not aimed at being instantiated, but rather inherited.

The sensors and actuators are attached to a FlowSolver object as a list, embedded in a ParamControl object. By attaching several sensors or actuators, it is possible to generate Multiple-Input, Multiple-Output configurations for control. The call to Sensors and Actuators is made automatically by FlowSolver.

For the cylinder case, we give an example below. We create two actuators forcing boundary conditions (on the top and bottom poles of the cylinder, respectively), and three point probes at different locations in the wake. They are gathered in a ParamControl object, which is passed as an argument to initialize a CylinderFlowSolver.

# Actuators
actuator_bc_1 = ActuatorBCParabolicV(angular_size_deg=10)
actuator_bc_2 = ActuatorBCParabolicV(angular_size_deg=10)
# Sensors
sensor_feedback = SensorPoint(sensor_type=SENSOR_TYPE.V, position=np.array([3, 0]))
sensor_perf_1 = SensorPoint(sensor_type=SENSOR_TYPE.V, position=np.array([3.1, 1]))
sensor_perf_2 = SensorPoint(sensor_type=SENSOR_TYPE.V, position=np.array([3.1, -1]))
# Gather actuators and sensors in ParamControl object
params_control = flowsolverparameters.ParamControl(
    sensor_list=[sensor_feedback, sensor_perf_1, sensor_perf_2],
    actuator_list=[actuator_bc_1, actuator_bc_2],
)

3️⃣ Run a (closed-loop) simulation

Once a use-case is defined by implementing the corresponding class inheriting FlowSolver, the basic feedback syntax has the following philosophy:

The FlowSolver subclass is instantiated with user-defined parameters
The base flow (stationary solution) is computed first
The object is prepared for time-stepping (e.g. we define operators, solvers, numerical schemes)
(Optional) A Controller is synthesized or read from a file
Time loop: iterate the FlowSolver.step(u) method, providing the 1D vector input u (open-loop or closed-loop using the Controller output)

A draft is given below. See the folder examples for more descriptive code.

# Instantiate and initialize FlowSolver object
fs = CylinderFlowSolver(...)
fs.compute_steady_state(...)
fs.initialize_time_stepping(...)

# Instantiate Controller (e.g. load from .mat file)
Kss = Controller.from_file(...)

# Time loop
y_meas = fs.y_meas
for _ in range(fs.params_time.num_steps):
    u_ctrl = Kss.step(y=-y_meas[0], dt=fs.params_time.dt)
    y_meas = fs.step(u_ctrl=u_ctrl)

The simulation should run seamlessly while providing information on the computed fields and potentially exporting information (as xdmf and csv).

More details

Initializations in general

Mesh

Mesh as xdmf. User is responsible for mesh coherence with respect to boundary conditions.

Initialization of a FlowSolver

Make all parameter structures from flowsolverparameters.py Timeline of function calls at __init__ for example.

All fields of a FlowSolver

Initialize time-stepping

initialize_time_stepping -> functions with two different ways, whether Tstart=0 (start from ic) or not (restart)

Perturbation of initial state:

_default_initial_perturbation(self, xloc: float = 0.0, yloc: float = 0.0, radius: float = 1.0) -> dolfin.Function

Base flow computation

Good practice is to do Picard iterations, then Newton

Initial guess for Picard only:

_default_steady_state_initial_guess(self) -> dolfin.UserExpression

Inlet flow profile

By default uniform u=1, can be changed in make_BCs. The perturbation bc is still 0 on the inlet.

FlowField, FlowFieldCollection, self.fields

Utility class to gather fields (velocity, pressure and mixed)

Actuators

Principle Actuator is an abstract class that encapsulates a dolfin.Expression and other parameters. An actuator are passed as a parameter to a FlowSolver for instantiation, through an actuator_list in the ParamControl object.

Actuator have an assigned type, defined as an integer enumeration: ACTUATOR_TYPE(IntEnum). It may be one of the following:

ACTUATOR_TYPE.FORCE: the actuator provides a volumic forcing. Its expression is automatically included in the momentum equation (in variational form).
ACTUATOR_TYPE.BC: the actuator modifies the boundary conditions dynamically. It should be reflected by the user when overriding _make_boundaries(), _make_bcs(). An example can be found in examples/cylinder/cylinderflowsolver.py:

def _make_bcs(self):
    ...
    bcu_actuation_up = dolfin.DirichletBC(
        self.W.sub(0),
        self.params_control.actuator_list[0].expression,
        self.get_subdomain["actuator_up"],
        )
    bcu_actuation_lo = dolfin.DirichletBC(
        self.W.sub(0),
        self.params_control.actuator_list[1].expression,
        self.get_subdomain["actuator_lo"],
        )
    ...
    return BoundaryConditions(bcu=bcu, bcp=[])

The expression of each actuator needs to be loaded after the FlowSolver is instantiated (the analytic dolfin.Expression is projected onto the FEM function spaces), which is handled automatically by the code.

Examples of actuators

ActuatorBCParabolicV: boundary condition actuator, 2nd component on velocity has parabolic profile

Mathematical expression:

$${v_{act}}({x}, t) = - \dfrac{(x_1-l)(x_1+l)}{l^2} u(t)$$, with $l = \frac{1}{2} D \sin \left( \frac{\delta}{2} \right)$ and $\delta$ is the tunable actuator opening in degrees.

FEniCS syntax:

def load_expression(self, flowsolver):
    L = (
        1
        / 2
        * flowsolver.params_flow.user_data["D"]
        * np.sin(1 / 2 * self.angular_size_deg * dolfin.pi / 180)
    )
    expression = dolfin.Expression(
        [
            "0",
            "(x[0]>=L || x[0] <=-L) ? 0 : u_ctrl * -1*(x[0]+L)*(x[0]-L) / (L*L)",
        ],
        element=flowsolver.V.ufl_element(),
        L=L,
        u_ctrl=0.0,
    )

    self.expression = expression

ActuatorForceGaussianV: force actuator, gaussian-shaped on the 2nd component of velocity

Mathematical expression:

$$B({x})u(t)=\left[ 0, \eta \exp\left( \frac{\left(x_1 - x_1^0\right)^2 + \left(x_2 - x_2^0\right)^2}{2\sigma_0^2} \right)\right]^T u(t)$$ with $\eta$ such that $\int_\Omega B({x})^T B({x}) d\Omega = 1$.

FEniCS syntax:

def load_expression(self, flowsolver):
    expression = dolfin.Expression(
        [
            "0",
            "u_ctrl * eta*exp(-0.5*((x[0]-x10)*(x[0]-x10)+(x[1]-x20)*(x[1]-x20))/(sig*sig))",
        ],
        element=flowsolver.V.ufl_element(),
        eta=1,
        sig=self.sigma,
        x10=self.position[0],
        x20=self.position[1],
        u_ctrl=1.0,
    )

    BtB = dolfin.norm(expression, mesh=flowsolver.mesh)
    expression.eta = 1 / BtB
    expression.u_ctrl = 0.0
    self.expression = expression

Define new actuators One can readily define a new actuator by inheriting the base class Actuator and providing a dedicated expression through the load_expression(self, flowsolver) method.

⚠️ Do not forget

Include a u_ctrl field in the dolfin.Expression. Its value may be changed in the body of the method (see ActuatorForceGaussianV), but it should be 0.0 when the method load_expression(self, flowsolver) exits.
Assign self.expression = expression at the end of def load_expression(self, flowsolver).

Sensors

Principle Sensor is an abstract class that gathers a behavior common to all sensors: it exhibits an abstract method eval(self, up: dolfin.Function) -> float to evaluate the measurement on a mixed-field (u,p).

The Sensor abstract class is expecte to be inherited by specific kinds of sensors. For example, the classes SensorPoint (point probe) and SensorHorizontalWallShear (integration on a subdomain) are subclasses that implement the Sensor.eval() abstract method.

The evaluation of sensors is handled automatically by the FlowSolver in the step() method.

Some sensors, for example those inheriting from SensorIntegral (e.g. SensorHorizontalWallShear) need to be loaded in some way (e.g. to define a subdomain of integration). They implement a load() method that is called by FlowSolver if the boolean Sensor.require_loading is set to True.

Just like actuators, sensors hold a SENSOR_TYPE(IntEnum), but it serves a different purpose. The SENSOR_TYPE.U, SENSOR_TYPE.V, SENSOR_TYPE.P, SENSOR_TYPE.OTHER is merely a shortcut to evaluate point probes on the right component of the field.

Examples of sensors

A simple SensorPoint(Sensor) has a straightforward definition of its eval() method: it evaluates the field at the given position (self.position) and on the given component (self.sensor_type):

Mathematical expression:

$y(t) = u(x_s, t)$ or $y(t) = v(x_s, t)$ or $y(t) = p(x_s, t)$ where $x_s$ is the sensor location.

FEniCS syntax:

def eval(self, up):
    return up(self.position[0], self.position[1])[self.sensor_type]

For a SensorHorizontalWallShear(SensorIntegral) (where SensorIntegral inherits directly from Sensor), the definition of the eval() method is more complex. First of all, the load() function defines a subdomain of integration with a given index: int and an associated ds: dolfin.Measure. The eval() method integrates (assemble) on the sensor subdomain (self.ds(int(self.sensor_index))) the quantity $\frac{\partial u}{\partial x_2}$.

Mathematical expression:

$y(t) = \int_{x \in S} \frac{\partial u}{\partial x_2} dx$

FEniCS syntax:

def eval(self, up):
    return dolfin.assemble(up.dx(1)[0] * self.ds(int(self.sensor_index)))

Define new sensors

New sensors can be defined by inheriting existing classes.

⚠️ The user is responsible for the compatibility of their Sensor.eval()code with parallel execution (MPI) of the code.

Saving files

save every

Controller

What is this class

control.StateSpace encapsulating a current state x and potentially a file

What is `step`

Wrapper around control.StateSpace.forced_response

Code: advanced use

Restart (show graph)
Modify the equations, the numerical schemes and the solvers used for the time simulation,
Use it as backend in an optimization tool (as in Jussiau, W., Demourant, F., Leclercq, C., & Apkarian, P. (2025). Control of a Class of High-Dimensional Nonlinear Oscillators: Application to Flow Stabilization. IEEE Transactions on Control Systems Technology.).

Code: utility

Operator computation
Frequency response computation
Export utils (spy, save)
Debug utils (export subdomains)
Mpi utils (= shortcuts)

Powered by GitHub

Navigation Menu

Code: sandbox

Conventions in the code

Code: basic use

Basic use

1️⃣ Choose a use-case: inherit the `FlowSolver` abstract class

2️⃣ Attach sensors and actuators to an instance of a `FlowSolver` subclass

3️⃣ Run a (closed-loop) simulation

More details

Initializations in general

Mesh

Initialization of a FlowSolver

Initialize time-stepping

Base flow computation

Inlet flow profile

FlowField, FlowFieldCollection, self.fields

Actuators

Sensors

Saving files

save every

Controller

What is this class

What is `step`

Code: advanced use

Code: utility

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally

Code: sandbox

Conventions in the code

Code: basic use

Basic use

1️⃣ Choose a use-case: inherit the FlowSolver abstract class

2️⃣ Attach sensors and actuators to an instance of a FlowSolver subclass

3️⃣ Run a (closed-loop) simulation

More details

Initializations in general

Mesh

Initialization of a FlowSolver

Initialize time-stepping

Base flow computation

Inlet flow profile

FlowField, FlowFieldCollection, self.fields

Actuators

Sensors

Saving files

save every

Controller

What is this class

What is step

Code: advanced use

Code: utility

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally

1️⃣ Choose a use-case: inherit the `FlowSolver` abstract class

2️⃣ Attach sensors and actuators to an instance of a `FlowSolver` subclass

What is `step`