Skip to content

Code: sandbox

Jump to bottom
William Jussiau edited this page Feb 24, 2025 · 73 revisions

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

Under construction

Conventions

  • 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 be able 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:

  1. The FlowSolver subclass is instantiated with user-defined parameters
  2. The base flow (stationary solution) is computed first
  3. The object is prepared for time-stepping (e.g. we define operators, solvers, numerical schemes)
  4. (Optional) A Controller is synthesized or read from a file
  5. 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).

See examples for a more detailed description.

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

Defining 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 Expression. Its value should be 0.0 when the method load_expression(self, flowsolver) exits. Its value may be set to something else in the body of the method (see ActuatorForceGaussianV),
  • Assign self.expression = expression at the end of def load_expression(self, flowsolver).

Sensors

Principle

Sensor is an abstract class that provides a method eval(self, up: dolfin.Function) -> float. Classes SensorPoint (point probe) and SensorIntegral (integration on a subdomain) are examples of subclasses that implement the eval method.

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

Defining new sensors

User is responsible for parallel execution when implementing new sensors.

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)
  • See other bullet points in Home
## Additional uses of the toolbox
The toolbox provides additional utility related to flow control:
* Compute dynamic operators A, B, C, D and mass matrix E,
* Restart a simulation from a saved file,
* Define an arbitrary number of actuators and sensors (e.g. feedback and performance sensors),
* Export time series (measurements from sensors, perturbation kinetic energy...) and fields for visualization,
* Modify the equations, the numerical schemes and the solvers used for the time simulation,
* Leverage parallel execution native to FEniCS,
* 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.](https://ieeexplore.ieee.org/abstract/document/10884641/)).

Code: utility

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

Powered by GitHub

Clone this wiki locally