diff --git a/README.md b/README.md index c7f8211..4ef5ec4 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black) ## Gym-preCICE Tutorials -This repository contains ready-to-run tutorial cases serve as for the Gym-preCICE adapter. These tutorials can be used as a foundation for creating similar control cases. +This repository contains ready-to-run tutorial cases that use the Gym-preCICE adapter. These tutorials can be used as a foundation for creating similar control cases. To define new control cases, you need to follow a simple file structure with three key components: ``` diff --git a/quickstart/baseline_uncontrolled.py b/quickstart/baseline_uncontrolled.py index cf1b64f..a86c109 100644 --- a/quickstart/baseline_uncontrolled.py +++ b/quickstart/baseline_uncontrolled.py @@ -37,8 +37,7 @@ def act(self): # step through the environment and control it for one complete episode (8 seconds, 320 steps) while not terminated: - with torch.no_grad(): - action = controller.act() + action = controller.act() _, _, terminated, _, _ = env.step(action) print("The baseline case is done.") diff --git a/quickstart/closed_loop_ppo_controller.py b/quickstart/closed_loop_ppo_controller.py index 6f2f837..43af173 100644 --- a/quickstart/closed_loop_ppo_controller.py +++ b/quickstart/closed_loop_ppo_controller.py @@ -52,7 +52,7 @@ def act(self, obs): obs, info = env.reset() print("\n...") - print("The control case is running!") + print("The PPO control case is running!") print( "This task is expected to be completed in about 5 minutes on a system with two cores @ 2.10GHz." ) diff --git a/quickstart/environment.py b/quickstart/environment.py index 2fb478e..50c98eb 100644 --- a/quickstart/environment.py +++ b/quickstart/environment.py @@ -1,3 +1,4 @@ +"""AFC environment for jet cylinder control.""" import logging import math from os.path import join @@ -18,13 +19,50 @@ class JetCylinder2DEnv(Adapter): - def __init__(self, options, idx=0) -> None: + r"""An AFC environment for jet cylinder control. + + ## Description + This environment corresponds to the DRL-based control case described by Rabault, Jean, et al. in + ["Artificial neural networks trained through deep reinforcement learning discover control strategies for active flow control"](https://doi.org/10.1017/jfm.2019.62). + We control flow rate of a synthetic jet attached to the poles of a rigid cylinder immersed in a two-dimensional incompressible channel flow. + The goal is to reduce drag forces acting on the cylinder. + + ## Action Space + The action is a `ndarray` with shape `(1,)` which can take values in the range of `[-2.5e-4, 2.5e-4]` with the value corresponding to + the control flow rate of the synthetic jet. + **Note**: Each control action (flow rate of the synthetic jet) is smoothly and linearly distributed over `self.action_interval` simulation time steps. + + ## Observation Space + The observation is a `ndarray` with shape `(151,)` which can take values in the range of `[-Inf, Inf]` with the values corresponding to + pressure sensor probes allocated within the flow filed. + + ## Rewards + Since the goal is to keep the drag forces acting on the cylinder as minimum as possible, a reward is defined based on the following relation: + reward = - 0.2 * ||, where and are respectively drag and lift coefficients of the cylinder averaged over the latest 0.335 s simulation time. + including the termination step, is allotted. The threshold for rewards is 500 for v1 and 200 for v0. + + ## Episode End + The episode ends after two seconds of flow field simulation. + + Args: + Adapter: gymprecice adapter super-class + """ + + def __init__(self, options: dict = None, idx: int = 0) -> None: + """Environment constructor. + + Args: + options: a dictionary containing the information within gymprecice-config.json. It is a return of `gymprecice.utils.fileutils.make_result_dir` method called within the controller algorithm. + idx: environment index. + """ super().__init__(options, idx) self.cylinder_origin = np.array([0, 0, 0.0223788]) self.cylinder_radius = 0.05 - self.jet_angle = [90, 270] - self.jet_width = [10, 10] + jet_angle = [90, 270] + jet_width = [10, 10] + self.theta0 = [np.radians(x) for x in jet_angle] + self.w = [np.radians(x) for x in jet_width] self._min_jet_rate = -2.5e-4 self._max_jet_rate = 2.5e-4 @@ -73,35 +111,58 @@ def __init__(self, options, idx=0) -> None: for solver_name in self._solver_list: if solver_name.rpartition("-")[-1].lower() == "openfoam": openfoam_case_name = solver_name - break - self._openfoam_solver_path = join(self._env_path, openfoam_case_name) - interface_patches = get_interface_patches( + + openfoam_interface_patches = get_interface_patches( join(openfoam_case_name, "system", "preciceDict") ) - actuators = [] - for patch in interface_patches: - if patch in self._actuator_list: - actuators.append(patch) - self.actuator_geometric_data = get_patch_geometry(openfoam_case_name, actuators) - actuator_coords = [] + action_patch = [] + self.action_patch_geometric_data = {} + for interface in self._controller_config["write_to"]: + if interface in openfoam_interface_patches: + action_patch.append(interface) + self.action_patch_geometric_data = get_patch_geometry( + openfoam_case_name, action_patch + ) + action_patch_coords = {} + for patch_name in self.action_patch_geometric_data.keys(): + action_patch_coords[patch_name] = [ + np.delete(coord, 2) + for coord in self.action_patch_geometric_data[patch_name]["face_centre"] + ] - for patch_name in self.actuator_geometric_data.keys(): - actuator_coords.append( - [ - np.delete(coord, 2) - for coord in self.actuator_geometric_data[patch_name]["face_centre"] + observation_patch = [] + for interface in self._controller_config["read_from"]: + if interface in openfoam_interface_patches: + observation_patch.append(interface) + self.observation_patch_geometric_data = get_patch_geometry( + openfoam_case_name, observation_patch + ) + observation_patch_coords = {} + for patch_name in self.observation_patch_geometric_data.keys(): + observation_patch_coords[patch_name] = [ + np.delete(coord, 2) + for coord in self.observation_patch_geometric_data[patch_name][ + "face_centre" ] - ) - self._set_precice_vectices(actuator_coords) + ] + + patch_coords = { + "read_from": observation_patch_coords, + "write_to": action_patch_coords, + } + + self._set_precice_vectices(patch_coords) @property def n_probes(self): + """Get the number of pressure probes.""" return self._n_probes @n_probes.setter def n_probes(self, value): + """Set the number of pressure probes.""" self._n_probes = value self._observation_info["n_probes"] = value self.observation_space = gym.spaces.Box( @@ -110,20 +171,24 @@ def n_probes(self, value): @property def n_forces(self): + """Get the number of data columns in forces file (excluding the time column).""" return self._n_forces @n_forces.setter def n_forces(self, value): + """Set the number of data columns in forces file (excluding the time column).""" assert value >= 3, "Number of forceCoeff columns must be greater than 2" self._n_forces = value self._reward_info["n_forces"] = value @property def min_jet_rate(self): + """Get the minimum flow rate of the synthetic jet.""" return self._min_jet_rate @min_jet_rate.setter def min_jet_rate(self, value): + """Set the minimum flow rate of the synthetic jet.""" self._min_jet_rate = value self.action_space = gym.spaces.Box( low=value, high=self._max_jet_rate, shape=(1,), dtype=np.float32 @@ -131,10 +196,12 @@ def min_jet_rate(self, value): @property def max_jet_rate(self): + """Get the maximum flow rate of the synthetic jet.""" return self._max_jet_rate @max_jet_rate.setter def max_jet_rate(self, value): + """Set the maximum flow rate of the synthetic jet.""" self._max_jet_rate = value self.action_space = gym.spaces.Box( low=self._min_jet_rate, high=value, shape=(1,), dtype=np.float32 @@ -142,10 +209,12 @@ def max_jet_rate(self, value): @property def latest_available_sim_time(self): + """Get the starting time of the flow field simulation.""" return self._latest_available_sim_time @latest_available_sim_time.setter def latest_available_sim_time(self, value): + """Set the starting time of the flow field simulation.""" if value == 0.0: value = int(value) self._latest_available_sim_time = value @@ -156,13 +225,18 @@ def latest_available_sim_time(self, value): self._prerun_data_required = value > 0.0 def step(self, action): + r"""Distribute the control action smoothly and linearly over `self.action_interval` simulation time steps.""" return self._smooth_step(action) def _get_action(self, action, write_var_list): - velocity = self._action_to_patch_field(action) - return {write_var_list[0]: velocity} + acuation_interface_field = self._action_to_patch_field(action) + write_data = { + var: acuation_interface_field[var.rpartition("-")[-1]] + for var in write_var_list + } + return write_data - def _get_observation(self): + def _get_observation(self, read_data, read_var_list): return self._probes_to_observation() def _get_reward(self): @@ -182,22 +256,24 @@ def _close_external_resources(self): raise err def _action_to_patch_field(self, action): - theta0 = [np.radians(x) for x in self.jet_angle] - w = [np.radians(x) for x in self.jet_width] origin = self.cylinder_origin radius = self.cylinder_radius patch_flow_rate = [-action, action] # velocity field of the actuation patches - U = [] - for idx, patch_name in enumerate(self.actuator_geometric_data.keys()): + U_profile = {} + for idx, patch_name in enumerate(self.action_patch_geometric_data.keys()): patch_ctr = np.array( - [radius * np.cos(theta0[idx]), radius * np.sin(theta0[idx]), origin[2]] + [ + radius * np.cos(self.theta0[idx]), + radius * np.sin(self.theta0[idx]), + origin[2], + ] ) - magSf = self.actuator_geometric_data[patch_name]["face_area_mag"] - Sf = self.actuator_geometric_data[patch_name]["face_area_vector"] - nf = self.actuator_geometric_data[patch_name]["face_normal"] - w_patch = w[idx] + magSf = self.action_patch_geometric_data[patch_name]["face_area_mag"] + Sf = self.action_patch_geometric_data[patch_name]["face_area_vector"] + nf = self.action_patch_geometric_data[patch_name]["face_normal"] + w_patch = self.w[idx] # convert volumetric flow rate to a sinusoidal profile on the interface avg_U = (patch_flow_rate[idx] / np.sum(magSf)).item() @@ -219,11 +295,12 @@ def _action_to_patch_field(self, action): Q_final = sum([u.dot(s) for u, s in zip(U_patch, Sf)]) if np.isclose(Q_final, patch_flow_rate[idx]): - U.append(U_patch) + U_profile[patch_name] = np.array( + [np.delete(item, 2) for item in U_patch] + ) else: raise Exception("Not a synthetic jet: Q_jet1 + Q_jet2 is not zero") - U_profile = np.array([np.delete(item, 2) for sublist in U for item in sublist]) return U_profile def _probes_to_observation(self): diff --git a/quickstart/open_loop_sinusoidal_controller.py b/quickstart/open_loop_sinusoidal_controller.py index 1168700..b3e0cf3 100644 --- a/quickstart/open_loop_sinusoidal_controller.py +++ b/quickstart/open_loop_sinusoidal_controller.py @@ -37,7 +37,7 @@ def act(self): _, _ = env.reset() print("\n...") - print("The control case is running!") + print("The sinusoidal control case is running!") print( "This task is expected to be completed in about 5 minutes on a system with two cores @ 2.10GHz." ) diff --git a/quickstart/physics-simulation-engine/gymprecice-config.json b/quickstart/physics-simulation-engine/gymprecice-config.json index 47b61e8..ece1909 100644 --- a/quickstart/physics-simulation-engine/gymprecice-config.json +++ b/quickstart/physics-simulation-engine/gymprecice-config.json @@ -1,14 +1,18 @@ { "environment": { - "name": "jet_cylinder_2d" + "name": "jet_cylinder" }, - "solvers": { - "name": ["fluid-openfoam"], + "physics_simulation_engine": { + "solvers": ["fluid-openfoam"], "prerun_script": "prerun.sh", "reset_script": "reset.sh", "run_script": "run.sh" }, - "actuators": { - "name": ["jet1", "jet2"] + "controller": { + "read_from": {}, + "write_to": { + "jet1": "Velocity", + "jet2": "Velocity" + } } } diff --git a/quickstart/quickstart.ipynb b/quickstart/quickstart.ipynb index f2bd586..e895f60 100644 --- a/quickstart/quickstart.ipynb +++ b/quickstart/quickstart.ipynb @@ -101,13 +101,15 @@ "metadata": {}, "outputs": [], "source": [ - "from os import path\n", + "from os import path, getcwd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import sys\n", + "sys.path.append(getcwd())\n", "\n", - "baseline_case = \"./fluid-openfoam\"\n", - "ppo_case = \"./fluid-openfoam\"\n", - "sinusoidal_case = \"./fluid-openfoam\"\n", + "baseline_case = \"gymprecice-run/jet_cylinder_baseline/env_0/fluid-openfoam\"\n", + "ppo_case = \"gymprecice-run/jet_cylinder_ppo/env_0/fluid-openfoam\"\n", + "sinusoidal_case = \"gymprecice-run/jet_cylinder_sinusoidal/env_0/fluid-openfoam\"\n", "\n", "# control signal information for the baseline (nothing but zero values, read as ckeck!)\n", "baseline_jet1_file = path.join(baseline_case, \"postProcessing/flowRateJet1/0/surfaceFieldValue.dat\")\n", @@ -172,9 +174,9 @@ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", - "baseline_case = \"./gymprecice-run/jet_cylinder_2d_baseline/env_0/fluid-openfoam\"\n", - "ppo_case = \"./gymprecice-run/jet_cylinder_2d_ppo/env_0/fluid-openfoam\"\n", - "sinusoidal_case = \"./gymprecice-run/jet_cylinder_2d_sinusoidal/env_0/fluid-openfoam\"\n", + "baseline_case = \"./gymprecice-run/jet_cylinder_baseline/env_0/fluid-openfoam\"\n", + "ppo_case = \"./gymprecice-run/jet_cylinder_ppo/env_0/fluid-openfoam\"\n", + "sinusoidal_case = \"./gymprecice-run/jet_cylinder_sinusoidal/env_0/fluid-openfoam\"\n", "\n", "# control response information for the baseline (nothing but zero values, read as ckeck!)\n", "baseline_forces_file = path.join(baseline_case, \"postProcessing/forceCoeffs/0/coefficient.dat\")\n",